Vibrational embedding theory

Many vibrational phenomena of interest in complex (bio-) molecular systems are of local nature, for example, the vibrational spectral signatures of water wires in bacteriorhodopsin.¹ For the computation of these spectra, subspace approaches are often applied,^2–7 accounting only for the vibrational degrees of freedom within a certain set of atoms of the system. The underlying potential energy surfaces may be harmonic^2–5 or anharmonic^6,7 and can be obtained in cluster models^2,4 or accounting for environmental atoms in a multi-scale manner.^2,3,6,7 Still, the description of the vibrational structure treats all considered vibrational degrees of freedom equally. In most cases, however, we are not interested in all vibrational degrees of freedom to the same extent. For instance, in the case of bacteriorhodopsin, the spectral region dominated by N–H and O–H stretches is of particular interest. Similarly, also in vibrational spectra with non-local contributions, typically not all vibrational degrees of freedom are of equal interest. For example, in the vibrational spectra of proteins, the so-called amide-I bands are often used as marker bands. In the computation of those bands by vibrational-structure theory, subspace approaches that only consider the vibrational modes of interest are adopted.^8–12 

These subspace approaches may be considered analogs to the cluster approaches in electronic-structure theory. Multi-level methods, instead, focus the effort on the modes important for the spectral region of interest while considering the remaining degrees of freedom on less rigorous footings. These approaches promise an improved accuracy–cost ratio, similar to their successful analogs in electronic structure theory. Such approaches are, however, sparse in vibrational-structure theory,^13–19 as have recently been reviewed in Ref. 20. Most of these methods have been developed in view of highly accurate calculations on small systems, e.g., a single water molecule,^13,19 although the selection of the mode–mode correlation toward specific modes of interest is also reported.^16,17

In this work, we suggest a vibrational embedding scheme, which allows the embedding of an anharmonic vibrational wave function in a mean-field description of the other vibrational degrees of freedom. The embedding is achieved by the addition of an embedding potential to the Hamiltonian, which may be obtained from a mean-field or correlated description of the environmental vibrational space. The resulting embedded Hamiltonian has the same form as the subspace analog. Hence, any vibrational wave function may be employed with this Hamiltonian, including local-correlated^21–27 and response treatments.^28–31 We show this for vibrational coupled-cluster response treatments, which have, to date, not been applied within multi-level vibrational structure approaches.

We first outline the theoretical background and methodology (Sec. II) and then present our implementation in MidasCpp (Sec. III). After the computational details (Sec. IV), we present numerical examples of small molecules and water clusters, as well as a water wire in bacteriorhodopsin (Sec. V). Finally, we summarize, conclude on the results, and present an outlook for further developments (Sec. VI).

As mentioned above, in realistic calculations of all but the smallest systems, the configurational space needs to be restricted. The most common approach here is to limit the number of modes in non-reference occupations. This number is given in square brackets behind the method acronym. For instance, VCI[2] refers to a VCI ansatz, in which only configurations with up to two non-reference occupations are considered. Analogously, the cluster operator for VCC[2] only includes excitation operators that substitute modals in up to two modes. Higher models, such as VCC[3] and VCC[4], are defined analogously.

For the calculation of vibrational spectra, response formalisms^28–31 offer a valuable alternative to the optimization of each individual vibrational excited state. The vibrational coordinates, i.e., the basis, the vibrational space is spanned in, is decisive, when dividing the vibrational space in multi-level approaches. Many different vibrational coordinates have been suggested in recent years (see the Ref. 20 for an overview of rectilinear vibrational coordinates). In the present work, we apply potential energy surfaces in the following types of rectilinear vibrational coordinates:

• Normal modes: Eigenvectors of the multi-dimensional harmonic problem.

• Optimized coordinates: Coordinates that minimize the VSCF energy for a given potential energy surface.^37–40 

• Hybrid optimized and localized coordinates (HOLCs): Coordinates optimized with a mixed optimization criterion consisting of a localization⁴¹ and an (in this case harmonic) energy criterion.⁴² 

• Group-localized coordinates (or local normal modes⁴³): These coordinates are obtained from the partial Hessian approach,^44–47 where the Hessian is diagonalized in blocks belonging to disjoint subsets of atoms. Within group-localized coordinates, vibrations moving the groups relative to each other are not well defined.

• FALCON coordinates: Most FALCON coordinates⁴⁸ are strictly local to a well-defined set of atoms in the system. Delocalized coordinates are only needed to describe the relative motion between the sets of atoms.

By choosing the MCR, the accuracy of the vibrational wave function can, hence, be stirred to be more accurate for selected degrees of freedom that reflect certain vibrational features. Similarly, Refs. 16 and 17 restrict a perturbative treatment (not correlated wave function) to focus on O–H stretching vibrations. In this work, all remaining modes are correlated with that mode of interest, but no correlation between environmental modes is considered.

Also correlated or non-correlated response calculations may be performed in an embedded subspace. These embedded-subspace response calculations are, however, not equivalent to the response calculations in full space with the corresponding restriction in the mode-combination range: In the full-space formulation, the environmental modals contribute to the vibrational excitations, which is not possible in the embedded-subspace formulation.

t_A, and $t̄A$ collect the cluster amplitudes and Lagrangian multipliers, respectively, in subset A. $ΛB$⁠, t_B, and $t̄B$ are defined analogously for subset B.

Hence, we can also address the VCC environment by means of an effective embedding potential, which can be used analogously to the VSCF-embedding potential in the reference state as well as response calculations.

Since the embedding potential of subset A $(vAemb,VCC)$ depends on the VSCF modals and VCC amplitudes and Lagrangian multipliers of subsystem B, the embedded VCC equations for the two subsystems are coupled. In the embedded-subset calculation of subset A, however, the embedding potential $vAemb,VCC$ is kept constant, i.e., neither the amplitudes nor the Lagrangian multipliers of the subsystem B are altered in this calculation. Analogously, the $vBemb,VCC$ is kept constant in the embedded-subset calculation of subset B. Still, the coupled equations can be solved in an iterative manner by using the results of the embedded subset calculation for subset A to update the embedding potential $vBemb,VCC$ for subset B and vice versa until self-consistency. In the resulting VCCinVCC scheme, both VSCF modals and VCC amplitudes and Lagrangian multipliers are mutually optimized.

We have, hence, defined different flavors of vibrational embedding theory as follows:

• VCCinVSCF embeds a correlated VCC treatment in a VSCF mean-field description of the environmental modes using an additive potential. It also allows for local response calculations. This approach shares some characteristics with the partially separable vibrational self-consistent field approach,¹³ which can be considered a full VCI (FVCI) in VSCF approach, as well as the vibrational mean-field configuration interaction (VMFCI) approach.^14,15

• VCCinVCC is an extension to VCCinVSCF to multiple subsets treated by means of correlated VCC wave function, being mutually optimized in an iterative manner. This approach is related to the vibrational active space VSCF (VASSCF) approach by Mizukami and Tew.¹⁹ Since VCCinVCC also consists of subset problems, response calculation leads also here to vibrationally excited states that are local to this subset.

• GroupVCC uses the MCR to restrict the correlation of the reference vibrational wave function within certain groups of modes as outlined in Sec. II B. In this case, a reference-state wave function is obtained, which includes all vibrational degrees of freedom, and subsequent response calculations can account for resonances between the vibrational excitations of different subsets. Note that for multiple subsets of modes, the reference-state GroupVCC approach is similar but not equivalent to VCCinVCC, as in VCCinVCC also the VSCF reference is updated in the iterative process, but not in GroupVCC.

Hence, the generation of the VCC-embedding potentials requires only the one-mode integrals $hqb,qbb,t$ in natural modals as well as the respective occupation numbers. Both have been extracted from vibrational wave function calculations in a locally modified version of MidasCpp.⁵² With these, the embedding contribution was directly added to the subset potential. The respective developments in MidasCpp are included in the MidasCpp 2023.04.0 release.⁵³ The embedded Hamiltonian was then employed in the embedded vibrational wave function calculation within MidasCpp. The described workflow has been implemented on a Python framework denoted EmbVib.

All vibrational wave function calculations have been performed with a locally modified version of MidasCpp 2021.4.1.⁵² We employed a bspline basis set⁵⁴ with the highest order ten and a primary basis set density of 0.8. Eight VSCF modals per mode were considered in the correlated vibrational wave function calculations. The VCC calculations were performed with the so-called V3 transformer. All reported vibrational excitation energies are obtained with the response module in MidasCpp using targeting of fundamental excitations.⁵⁵ Within the VCCinVCC scheme, we applied five cycles, which proved to be sufficient to converge the vibrational excitation energies. For the methylfurfural, dicyclopropyl ketone, and maleimide examples, we employed Hamiltonian mode–mode coupling measures⁴⁹ to divide the vibrational spaces. For this analysis, we considered eight modals per mode.

All employed potential energy surfaces (PESs) are reported in earlier studies: For the methylfurfural and dicyclopropyl ketone examples, we employed the reference potential in FALCON coordinates reported in Ref. 56, which contains up to two-mode couplings (RI-MP2/cc-pVTZ). The calculations for maleimide were obtained for the potentials reported in Ref. 57 for normal modes, HOLCs, and FALCON coordinates. In all three cases, we here apply the potential up to four-mode couplings, with screened terms and polynomial fit (see Ref. 57 for details). The quartic potentials for the water dimer and trimer in normal modes and optimized coordinates contain only vibrational degrees of freedom dominated by the motions within the individual water molecules and have been generated by Yagi et al.³⁹ The potential energy surface for bacteriorhodopsin was also provided by Yagi and Sugita (Ref. 7), generated for the protein system only.

Often, only certain vibrational features are of interest rather than the full spectrum. C=O stretches are typical marker bands. We investigate the systematic increase of the vibrational subspace treated by means of VCC response in a VCCinVSCF scheme for the C=O stretching vibrations of methylfurfural and dicyclopropyl ketone. For these calculations, we employ the PESs in FALCON coordinates first reported in Ref. 56. In the applied vibrational coordinates, three groups of atoms are defined in which the internal vibrational degrees of freedom are strictly local within these groups of atoms. The initial groups are depicted in Figs. 1(a) and 1(b).

For the methylfurfural PES applied here, these groups are the CHO group, the furan ring as well as the methyl group. Additionally, the vibrations moving these groups relative to each other are defined in the FALCON scheme. In the present PES, this concerns two sets of inter-group vibrational coordinates: One set of six vibrational coordinates moving the CHO and furan groups relative to each other while the position of the atoms in the methyl group is fixed, and an additional set of six vibrational coordinates moving the CHO–furan entity with respect to the methyl group.

Assuming that the distance to the C=O group correlates with the coupling to the C=O stretching frequency, we set up the following hierarchy of vibrational spaces with an increasing number of degrees of freedom (DOF):

• S1: Only the C=O stretch (1 DOF)

• S2: Only the vibrations in the CHO group (3 DOFs)

• S3: S2 plus the relative vibrations to the furan ring (9 DOFs)

• S4: S3 plus the vibrations in the furan ring (23 DOFs)

• full: All vibrational degrees of freedom (36 DOFs)

Figure 2(a) depicts the VCC[2] C=O vibrational frequency for these vibrational spaces neglecting all other degrees of freedom (subspace approach) as well as within VCC[2]inVSCF embedding. Including only the C=O stretching vibration (S1), the subspace approach yields a vibrational excitation energy closer to the full system VCC[2] reference than the VCC[2]inVSCF embedded approach. Increasing the vibrational space by subsequently including the remaining vibrations in the CHO group and the relative motion to the furan group (S1 and S3) leads to a drop in vibrational frequency and thereby to a larger disagreement with the VCC[2] reference for both schemes. Only the inclusion of the internal degrees of freedom of the furan group (S4) leads to good agreement with the reference for the VCC[2] inVSCF case, whereas the subspace VCC[2] treatment for S4 overestimates the vibrational frequency.

This shows that the inclusion of internal vibrations of the furan ring is required. It is, however, fair to assume that not all vibrations of the furan ring contribute equally. We, hence, suggest alternative measures to estimate the coupling prior to the VCC calculation. We, here, apply the Hamiltonian measures for mode–mode couplings reported in Ref. 49. Notably, a clear correlation of the Hamiltonian measure to the FALCON-type hierarchy is obtained (see Table S-I supplementary material). The resulting fundamental frequencies of the C=O stretch in methylfurfural for different thresholds of that Hamiltonian measure are shown in Fig. 2(b). The measure-based expansion of the vibrational space exhibits an overall similar behavior as the FALCON-based expansion of the vibrational space. The smaller steps possible in the measure-based expansion, however, allow convergence already at a lower number of degrees of freedom included in the active subspace.

An analogous analysis for dicyclopropyl ketone is shown in Fig. 3. The atomic subsets in the applied FALCON coordinates are the C=O group as well as the individual cyclopropene rings [cf. Fig. 1(b)]. Due to the symmetry of the molecule, only three different partitionings appear sensible, these are as follows:

• S1: Only the C=O stretch (1 DOF)

• S2: C=O stretch and relative motions of the C=O group with respect to both cyclopropyl groups (12 DOFs)

• full: all 48 DOFs

The corresponding convergence of the fundamental vibrational excitation energy seems better for the subset approach than for the VCC[2]inVSCF embedding [cf. Fig. 3(a)]. Since the embedding approach incorporates more effects than the subspace approach, we assign this behavior to fortunate error cancellation. This picture is different for the measure-based expansion of the vibrational space [Fig. 3(b)], applying the same Hamiltonian-based measures as above. Although the one-mode subset result lies closer to the full-space VCC[2] reference, the VCC[2]inVSCF results are converged at a lower number of modes included in the VCC space than the subspace results.

As an example of targeting coupled C=O stretching vibrations, we have chosen the maleimide molecule, which contains two symmetrically equivalent C=O groups, whose stretching vibrations are highly coupled and form a symmetric and an asymmetric linear combination. For this example, we further explore different vibrational coordinates to span the vibrational space. These are normal coordinates, FALCON coordinates, and HOLC coordinates. The subgroups of modes for the FALCON treatment are the two C=O groups, individually, the NH group, and the HC=CH group [cf. Fig. 1(c)].⁵⁷ This means, in FALCON coordinates, the C=O stretching vibrations exclusively involve movements of C and O atoms, and none of the other atoms in the molecule move. Also, in the HOLC coordinates for maleimide, the C=O stretches are localized on the individual C=O groups, but all other atoms in that molecule can have (small) contributions. Since these coordinates span the vibrational space differently, the results for the different coordinate systems are only equivalent in a full-space treatment, which is clearly out of reach for the maleimide molecule. In Ref. 57, it is shown that converging the results for different coordinate systems for maleimide to the same limit is challenging. However, all types of coordinates support similar agreement with the experiment.⁵⁷ 

We note that for all investigated vibrational coordinates, we observe a high degree of mixing in the fundamental excitations, complicating the assignment of the modes. This is why we start the analysis with VCC[1] for the active subset. A consequence of this is that the reference results differ for the different coordinate types. It still allows us to assess the impact of the environmental modes. Below, we will also briefly discuss the respective VCC[2] results.

Figure 4 shows the convergence of the symmetric and asymmetric C=O stretching frequencies with increasing VCC[1] space in the reduced space and VCC[1]inVSCF response treatment, respectively. The subspace has been increased based on the Hamiltonian coupling parameters to either of the modes dominated by the C=O stretching vibrations. The tick “i” labels computations with one-mode subsets. Since local or localized C=O stretches in FALCON coordinates and HOLCs are symmetry equivalent, the same value is obtained for the symmetric and asymmetric vibrations, respectively. “CO” then denotes subspaces in which only these two vibrational degrees of freedom are considered together, and “full” denotes the space including all 24 vibrational degrees of freedom for maleimide. All other ticks correspond to the threshold of the Hamiltonian measure applied. The scale of the x-axis in Fig. 4 represents the number of modes considered in the subspace of interest.

For normal modes, we observe that the VCC[1]inVSCF treatment for one mode in the active space is already in quite close agreement with the full result. This holds particularly for the symmetric vibration, whereas for the asymmetric vibration, a somewhat alternating behavior is obtained with increasing subspace, which, however, is already converged for a threshold of 5 × 10⁻⁴ (corresponding to seven degrees of freedom). In the subspace approach in normal modes, the convergence of both vibrational energies is significantly slower than in the embedded case.

For HOLCs, we observe a similar behavior as for normal modes as soon as both vibrational modes of interest are correlated (CO). Also here, the frequency of the asymmetric stretch is roughly converged with a threshold of 5 × 10⁻⁴. In HOLCs, however, this corresponds to 13 vibrational modes in the VCC[1] subspace. Again, the convergence is quicker for the VCC[1]inVSCF approach than in the subspace approach. In FALCON coordinates, the convergence of the asymmetric stretch seems less smooth.

The respective subspace VCC[2] and VCC[2]inVSCF results are shown in Fig. 5. Here, a large degree of ambiguity occurs for all types of vibrational coordinates. This poses challenges in the assignment and comparison of the different vibrational coordinates. This is why several vibrational excitation energies are reported with the weights of the C=O stretching vibrations in Fig. 5. Note that the weights are not directly comparable for the different coordinate schemes, due to the different ways to span the vibrational space. A direct comparison, however, is not in the scope of the present work. We note that for VCI wave functions, a one-to-one comparison for different vibrational coordinates has been reported.⁵⁸ 

For all coordinate systems, the VCC[2] result is faster approached by the VSCF embedding than by the subspace approach. The differences in the full VCC[2] reference results for the different coordinates stem from non-converged results with respect to the mode–mode coupling expansion.⁴² 

As a first numerical example for multiple subsets, we have chosen inner vibrational degrees of freedom for water in a water dimer [cf. Fig. 1(d)] with the potential energy surfaces from Ref. 39. Table I shows the VCC[3] fundamental vibrational excitation energies for all internal degrees of freedom of the water dimer in normal modes for different embedding schemes with energy-based partitioning. This means we group the low-energy bending modes in one subset and the high-energy stretching modes in the other. The VCC[3] subset calculation leads to a deviation of 30–60 cm⁻¹. This deviation is reduced to −3 to 7 cm⁻¹ in VCC[3]inVSCF embedding. Additional mutual optimization in a VCC[3]inVCC[3] treatment leads to only minor improvement. Additionally, the results for a supersystem with VCC mode–mode correlations only within the given subsets (GroupVCC[3]) are given. This coupling does not improve significantly on the results in this case either.

In the water dimer, the normal modes can clearly be assigned to the individual water molecules. We have, hence, performed local partitioning, grouping all normal modes assigned to one water molecule. By this, we obtain a better performance of the subset approach and a slightly larger deviation of the embedding and GroupVCC approaches [see yellow graph in Fig. 6(a)]. Additionally, we conducted both partitionings also with optimized coordinates, for which the local character is more pronounced than for the normal modes. For energy-based partitioning [green graph in Fig. 6(a)], the two types of coordinates perform very similarly. In optimized coordinates with local partitioning [red graph in Fig. 6(a)], already the subset approach yields very good results.

For the water trimer [cf. Fig. 1(e)], the results, however, look rather different [see Fig. 6(b)]. Again, the energy-based partitioning for normal modes and optimized coordinates yields almost identical results. However, the deviations obtained for the local partitioning with all embedding schemes are with mean average deviations (MADs) above 15 cm⁻¹ rather high. The clearly best performance of the local partitioning is obtained for the GroupVCC[3] scheme. Here the local partitioning leads to better agreement with the reference than the energy-based partitioning. Figures 7(a)–7(c) show the deviations for bending vibrations, asymmetric stretches, and symmetric stretches, respectively. This shows that the performance of the local partitioning is clearly worst for the symmetric stretches, whereas the other vibrational transitions are better described. This result is expected due to the high symmetry of the water trimer: The symmetry of the fundamental vibrational excitations cannot be represented in the local embedding approaches. The GroupVCC[3] scheme, however, can recover the symmetric delocalized vibrational excitations.

As a second numerical example of multiple subsets, we have chosen 16 fundamental stretching and bending vibrations of the pentagonal hydrogen-bonded network in bacteriorhodopsin (Fig. 8), following the labeling of the modes in Ref. 7.

The chosen vibrational modes were previously selected by Yagi and Sugita in Ref. 7, where the potential energy surface was generated in group-localized coordinates. Employing this potential, we applied different vibrational embedding schemes (SubVCC[4], VCC[4]inVSCF, VCC[4]inVCC[4], and GroupVCC[4]). We performed a VCC[4] linear response calculation with the full set of the 16 modes focusing on all fundamental transition vibrations as reference. We observe two resonance states for two fundamental frequencies assigned to 18₁ and 29₁ in Ref. 7 with combination bands or overtones of modes within the same local entity (see Table S-X in the supplementary material for details). Here, we adopted the nomenclature from Ref. 7, in which 18₁ is the fundamental (first) vibrational excitation of mode 18, whereas 18₂ is the second vibrational excitation of the same mode. All other contributions are defined analogously. Figure 9 shows the 16 individual fundamental frequencies for the different vibrational embedding schemes and the VCC[4] reference with an energy-based partitioning (left) and local partitioning (right) of the modes. The mixed excitations described above are shown in dashed lines in Fig. 9. All values are also reported in Tables S-XI and S-XII in the supplementary material.

We first focus on the discussion of energy-based partitioning, where the full set of modes is divided into two subsets: one subset consists of all bending modes and the other includes all stretching modes. Within this partitioning of the full set of modes, the three embedding schemes SubVCC[4], VCC[4]inVSCF, and VCC[4]inVCC[4] exhibit only one vibrational eigenstate for the stretching fundamental frequencies corresponding to 18₁ and 29₁. These three embedding schemes are clearly unable to capture the resonances due to the lack of description of the correlation between the subsets and are, therefore, excluded from the MAD values discussed below.

The most significant individual deviations of the SubVCC[4] approach from the VCC[4] reference are obtained for the bending vibrations 46₁, 28₁, and 37₁ of the respective three water molecules in the hydrogen-bonded network (cf. Fig. 8). Adding an embedding potential, either through the VCC[4]inVSCF embedding or the VCC[4]inVCC[4] embedding significantly reduces the deviation by a factor of 10 (e.g., for 46₁ from 34.26 to 2.68 cm⁻¹). This is also obtained for the other fundamental frequencies, except for 59₁: Although with the SubVCC[4] approach, a minor deviation of 0.92 cm⁻¹ is observed, the addition of an embedding potential increases the deviation in the case of VCC[4]inVSCF and VCC[4]inVCC[4] by 6.39 cm⁻¹ in this case.

The MAD (cf. Table II) of the SubVCC[4] approach of all 14 fundamental frequencies exhibits a value of 16.51 cm⁻¹, whereas the MAD obtained from the bending vibrations only is slightly larger with 17.45 cm⁻¹ than the MAD obtained from the stretching vibrations only with 15.26 cm⁻¹. The MADs for bending, stretching, and all modes are significantly reduced with an embedding potential. In the case of the VCC[4]inVSCF embedding scheme, we obtain a MAD of 2.6 cm⁻¹, which is marginally higher than the MAD with the VCC[4]inVCC[4] embedding scheme with 2.3 cm⁻¹. The fundamental frequencies obtained with the GroupVCC[4] scheme exhibit the smallest MAD with 1.8 cm⁻¹ and obtained further the two resonance states as in the VCC[4] reference response calculation of the full set of modes.

Figure 9 and Table II further report the individual vibrational excitation energies and their MADs for local partitioning of the subsets from the unpartitioned VCC[4] reference. Here, the subsets are divided into groups analogously to the grouping of atoms in Ref. 7. The modes within one group are indicated by the rectangles in Fig. 8. The SubVCC[4], VCC[4]inVSCF, and GroupVCC[4] approaches are able to capture the two resonance states of the VCC[4] reference. The respective vibrational excitations show similar resonances; however, the weights of the fundamental transitions 18₁ and 29₁ are more dominant in comparison to the full VCC[4] reference (cf. Table S-X in the supplementary material). Therefore, we included the respective two frequencies of the resonance state with the lower energy in the MAD calculation, whereas the frequencies with higher energies are excluded in the MAD determination and are indicated with dotted lines in Fig. 9.

For SubVCC[4], the deviations for the fundamental modes range from 0.00 cm⁻¹ (8₁) to 21.57 cm⁻¹ (47₁). Similar deviations are obtained with an additional embedding potential (cf. Table S-XII in supplementary material), and the vibrational embedding approaches VCC[4]inVSCF and VCC[4]inVCC[4] show no significant improvement to the SubVCC[4]. Only the GroupVCC[4] shows a significant improvement with the largest deviation of 1.37 cm⁻¹ (29₁). The SubVCC[4] with local partitioning leads to a MAD including all 16 fundamentals of 2.84 cm⁻¹, which is significantly lower than with the SubVCC[4] approach using energy-based partitioning. This MAD value is not significantly lowered by VSCF or VCC embedding, but again, GroupVCC can improve the results.

For the bacteriorhodopsin, we, hence, find a superior performance in local partitioning compared with the energy-based partitioning similar to the water dimer discussed above. In this context, we recall that the applied potential energy surface considers higher coupling terms within a local group (grid-based two-mode coupled terms plus quartic force field terms up to three-mode coupling) than in between the group (harmonic terms only). This potential, hence, may tend to a higher sensitivity for local correlation than for inter-group correlation.

We further performed vibrational embedding calculations with the two partitioning schemes (energy-based and local) using lower-order VCC expansions. All MADs of the vibrational embedding calculations show similar overall trends for VCC[2] and VCC[3] compared with the VCC[4] vibrational embedding calculations (cf. Sec. S4.C in the supplementary material).

We have presented a framework for the vibrational embedding theory that incorporates the effect of environmental modes by means of a mean-field interaction. It can be applied to correlated as well as response subset treatments. The environmental description may either be obtained from a mean-field VSCF (VCCinVSCF) calculation of the environmental degrees of freedom or a correlated treatment (VCCinVCC). The latter also allows mutual optimization of the different subsets, which can be considered a multi-reference approach.¹⁹ These approaches have been implemented and assessed for local C=O vibrational excitations on methylfurfural and dicyclopropene and coupled C=O vibrations in maleimide, as well as for a water dimer, water trimer, and a water cluster in bacteriorhodopsin. In all cases, we compare with full reference calculations and those, where the grouping is only implemented in the reference vibrational wave function, but not in the response framework (GroupVCC).

We generally find VCCinVSCF to converge clearly quicker than the respective subspace calculations. An additional mutual optimization in a VCCinVCC framework has only a minor impact. This observation is in line with earlier findings for VCIinVCI embedding schemes.¹⁹ The inclusion of coupling in the response framework (GroupVCC), however, generally further improves the agreement with the full reference.

Notably, we observe a strong dependence on the performance of the different embedding and group correlation schemes for several groups on the applied partitioning: For the water dimer, local partitioning is found to be more beneficial than the energy-based analog. In contrast to that, embedding methods (VCCinVSCF and VCCinVCC) with local partitioning do not yield reasonable agreement with the reference for fundamental vibrational excitation in the symmetric water trimer, as the delocalized excitation cannot be recovered in a local embedding scheme. For the larger water cluster in bacteriorhodopsin, however, the local partitioning led to an astonishing close agreement in the subspace approach, with only minor changes by embedding (VCCinVSCF and VCCinVSCF). The GroupVCC approach leads to an improved agreement also in this case.

These findings suggest that a future two-step scheme establishing excitonic-coupling schemes for selected local vibrational excitations (e.g., fundamentals of the O–H stretches in the water trimer) on top of the embedding vibrational wave function is a promising route to also cover delocalization effects.

For the case of maleimide, we showed briefly that our analysis for the different vibrational coordinates is hampered by a high degree of excitation mixing. This complication may be circumvented by a non-resonant response framework^59,60 allowing us to calculate the absorption infrared cross-section at a certain frequency.^61–65 For this, however, also dipole surfaces are required. This means incorporating the environmental effect in these so-called damped response schemes^61–65 requires also embedding schemes for the dipole operator, which will be part of future work in our group.

See the supplementary material for additional data underlying the discussed calculations of methylfurfural (S1) in Sec. V A of the main text, dicyclopropyl ketone (S2) in Sec. V B of the main text, malmeimide (S3) in Sec. V C of the main text, and bacteriorhodopsin (S4) of Sec. V D of the main text. Furthermore, for bacteriorhodopsin results of embedded calculations, with the two partitioning schemes (local and energy-based) with lower order VCC (VCC[2] and VCC[3]) are reported and discussed.