Quantum-chemical calculation of two-dimensional infrared spectra using localized-mode VSCF/VCI

Computational protocols for the simulation of two-dimensional infrared (2D IR) spectroscopy usually rely on vibrational exciton models, which require an empirical parametrization. Here, we present an efficient quantum-chemical protocol for predicting static 2D IR spectra that does not require any empirical parameters. For the calculation of anharmonic vibrational energy levels and transition dipole moments, we employ the localized-mode vibrational self-consistent field (L-VSCF) / vibrational configuration interaction (L-VCI) approach previous established for (linear) anharmonic theoretical vibrational spectroscopy [Panek and Jacob, ChemPhysChem 15 , 3365–3377 (2014)]. We demonstrate that with an efficient expansion of the potential energy surface using anharmonic one-mode potentials and harmonic two-mode potentials, 2D IR spectra of metal carbonyl complexes and of dipeptides can be predicted reliably. We further show how the close connection between L-VCI and vibrational exciton models can be exploited to extract the parameters of such models from those calculations. This provides a novel route to the fully quantum-chemical parametrization of vibrational exciton model for predicting 2D IR spectra.


Introduction
2][3][4][5] Since its inception by Hamm, Lim, and Hochstrasser 25 years ago, 6 it has flourished into a broadly applicable experimental method. 7Applications in chemical biology include the elucidation of the structure 8,9 and dynamics of proteins [10][11][12][13][14][15][16] in solution.Recent experimental advances allow for studying biomolecules at interfaces, 17,18 enable transient spectroscopy for following non-equilibrium protein dynamics, 19 and implement spectro-electrochemistry to study redox-active proteins. 20,21 oreover, 2D IR spectroscopy is employed to investigate ground-state and excited-state dynamics in model systems such as metal carbonyl complexes [22][23][24][25] and can be applied in material science to shine light on the ion dynamics in lithium ion batteries. 26[29] Only by comparison to computational predictions, detailed information on the structure and dynamics can be extracted from experimental 2D IR spectra.Simply speaking, in 2D IR spectroscopy one performs a pump-probe experiment, in which a pump pulse excites a molecule from the vibrational ground state to a singly-excited vibrational state, followed by a probe pulse that results in stimulated emission or in an excitation to a doubly excited vibrational state. 4,30 his way, 2D IR spectroscopy reveals the anharmonicities of the vibrational modes as well as couplings between them.By spreading out the spectroscopic information in two dimensions, 2D IR spectroscopy can reveal information that is hidden in conventional vibrational spectra.Furthermore, by varying the delay between the pump and the probe pulses, dynamical information can be assessed.Because it specifically probes anharmonicity, for a computational treatment of 2D IR spectroscopy one has to go beyond the harmonic approximation 31 that is commonly used in theoretical vibrational spectroscopy. 2,4,32,33 Coon protocols for the simulation of 2D IR spectroscopy rely on empirical frequency maps in combination with vibrational exciton models (for an exhaustive review, see Ref. 34).Such models parametrize the frequencies and transition dipole moments of local vibrational modes as well as the harmonic couplings between them as functions of the molecular structure and the local electric field.Furthermore, the anharmonicity of the local modes is included as an empirical shift.This allows for the efficient calculation of the vibrational energy levels, which can be combined with molecular dynamics simulations and used as input for the simulation of 2D IR spectra.Such a protocol is efficient and has been successfully applied in numerous studies over the past 25 years.2][43] However, the construction of frequency maps for new systems can be tedious and its parametrization requires a significant effort.Therefore, the reliance on such parametrizations has slowed down the progress in 2D IR spectroscopy.
For instance, the interpretation of amide I 2D IR spectra of peptides containing proline required a dedicated parametrization. 44For promising 2D IR experiments studying proteins beyond the amide I vibrations 45 or investigating nucleic acid bases or DNA [46][47][48] suitable parametrizations became available only years after the initial experiments [49][50][51][52][53] and the interpretation of such experiments partly remains elusive.
Quantum-chemical methods of theoretical vibrational spectroscopy allow for the calculation of vibrational energy levels without the need for system-specific parametrizations and can thus close this gap.The calculation of vibrational spectra within the harmonic approximation is well-established and routinely applied for molecules with hundreds of atoms (see, e.g., Refs.[73][74] Here, we extend one such approach, the localized-mode vibrational self-consistent field (L-VSCF) method in combination the with localized-mode vibrational configuration interaction (L-VCI) method, 75,76 for the calculation of 2D-IR spectra.When applied together with an anharmonic one-mode potential energy surface with harmonic two-mode couplings between localized modes, this approach has been shown to be able to predict anharmonic vibrational spectra accurately and efficiently. 73,77 oreover, it allows for a direct connection with vibrational exciton models, 78 without the need for any empirical parametrization.
Previously, some attempts have been made to obtain 2D IR spectra from quantumchemical calculations.Mukamel and coworkers employed a fourth-order Taylor expansion of the potential-energy expansion as basis of an exciton Hamiltonian. 79Subsequently, second-order vibrational perturbation theory has been applied for modeling 2D IR spectra of metal carbonyls 23 and of deuterated glycolaldehyde. 80Sibert and co-workers constructed a model of harmonically-coupled Morse oscillators and applied it to model the 2D IR spectra of carbonyl complexes. 81More recently, Besley and coworkers 82 employed calculations in the harmonic approximation in combination with a localization of normal modes 83 to extract the parameters of an exciton model (local-mode frequencies and harmonic couplings).They then applied an empirical anharmonicity on top of this model to simulate the amide I 2D IR spectra of small peptides.Bowman, Tokmakoff, and coworkers employed high-level VSCF/VCI calculations in combination with a local monomer approximation for protonated water clusters in order to interpret 2D IR spectra. 84re, our goal is to devise a computational protocol that is efficient and parameter-free and that allows for both the direct prediction of 2D IR spectra for arbitrary molecular systems and can serve as reference for the parametrization of vibrational maps as needed for empirical exciton models.This work is organized as follows.In Section 2, we review the methodology applied for the quantum-chemical calculation of anharmonic vibrational energy levels and transition dipole moments with L-VSCF and L-VCI (Sections 2.1 and 2.2), discuss the connection to vibrational exciton models (Section 2.3), and recall the theory for the calculation of 2D IR spectra (Section 2.4).In Section 3 we give the details of our computational protocol.
In Section 4, we apply this protocol for selected test cases, specifically for two carbonyl complexes and two dipeptides.We compare the different approximations for the potential energy surface in L-VSCF/L-VCI and exciton models derived from the quantum-chemical calculations.Finally, concluding remarks and an outlook are presented in Section 5.

Harmonic approximation and localized modes
In the harmonic approximation, the starting point for the calculation of vibrational spectra is the calculation of the Hessian matrix H (m) , which contains the second derivatives of the potential energy surface with respect to mass-weighted Cartesian coordinates. 31,85 y diagonalizing this Hessian, one obtains the harmonic normal modes Q as eigenvectors with the eigenvalues H (q) ij , which are related to the vibrational frequencies ω i = H (q) ij , and give rise to the vibrational energy levels, where M is the number of vibrational degrees of freedom.In the double-harmonic approximation, the transition dipole moments for the fundamental transitions are given by the derivatives of the molecular dipole moment with respect to the normal-mode coordinates, i.e., µ i = dµ/dq i .
Using mass-weighted normal-mode coordinates q i , the harmonic approximation to the molecular potential energy surface is given by 7]86 To this end, all vibrational modes that contribute to a specific band or are of interest are collected into the subset Q sub .For 2D IR spectroscopy, this will usually be only a small number of vibrational modes, such as the carbonyl stretching vibrations in metal carbonyls or the amide I vibrations in peptides.For this subset of modes, a unitary transformation U is determined such that the resulting modes, are maximally localized, as measured by a suitable criterion. 83Here and in the following, a tilde denotes quantities that refer to localized modes.The Hessian with respect to localized-mode coordinates, is no longer diagonal, but the localization re-introduces non-diagonal elements that describe the harmonic coupling between localized modes, i.e., the harmonic approximation to the potential energy surface now becomes where ωi = Hii are the localized-mode frequencies and Hij are the off-diagonal elements of the Hessian with respect to localized modes.

L-VSCF and L-VCI
To account for anharmonic effects in theoretical vibrational spectroscopy, one needs to go beyond the harmonic approximation to the potential energy surface given by Eq. (3)   or Eq. ( 6).Commonly, this is done by using a truncated n-mode expansion, which we perform in localized-mode coordinates.When including up to two-mode contributions, this corresponds to 87, 88 Here, the one-mode potentials Ṽ (1) i (q i ) and the two-mode potentials Ṽ (2) ij (q i , qj ) are evaluated numerically on suitable grids of displacements along localized-mode coordinates. 75 a further approximation, the full two-mode potentials can be replaced by their harmonic counterparts, i.e., Note that for this approximation, the number of required single-point calculations scales linearly with the number of considered localized modes, as the harmonic two-mode couplings can be extracted from the Hessian with respect to localized modes.
Using these approximations, the vibrational Schrödinger equation can be solved using L-VSCF / L-VCI. 75,77 -VSCF uses a product ansatz for the vibrational wavefunction, with the modals (i.e., one-mode wavefunctions) ϕ n i i (q i ) for the i-th localized mode qi , with the corresponding vibrational quantum number n i .The modals can be obtained from the self-consistent solution of one-mode equations, 89 − 1 2 in which the effective potential V eff,0 i (q i ) depends on the one-mode wavefunctions.
In L-VCI, one then uses a configuration interaction expansion of the vibrational wavefunction as linear combination of the L-VSCF wavefunctions Ψ L-VSCF n 1 ,n 2 ,... , i.e., 90 Generally, one has to truncate the L-VCI expansion.Here, we use an L-VCISD expansion, which includes the ground state as well as singly and doubly excited states, which is the minimum necessary for modeling 2D IR spectroscopy.This approximation is usually sufficient when considering weakly coupled localized modes. 73,77,78 Biagonalizing the L-VCI matrix, one can obtain the expansion coefficients c (n) I as well as the corresponding energy eigenvalues E n .The anharmonic excitation frequencies are then given by ωn = E n − E 0 , where the bar is used to distinguish the anharmonic frequencies from the harmonic fundamental frequencies ω i .
Using the L-VCI expansion coefficients, the transition dipole moments between the corresponding vibrational states.
can be calculated with the help of a one-mode expansion of the dipole moment surfaces.
For a more detailed discussion of L-VSCF/L-VCI and well as benchmark calculations assessing its accuracy, we refer to Refs.73, 75, 77.

Vibrational exciton models
Vibrational exciton models, 4,6,34 commonly expressed in second quantization, approximate the vibrational Hamiltonian as, Here, b † i and bi are bosonic creation and annihilation operators with respect to a local-mode basis, such as the L-VSCF modals (i.e., b † i creates a vibrational quantum for the i-th L-VSCF modal).The exciton Hamiltonian contains as parameters the local-mode harmonic frequencies, ωi , the harmonic coupling constants J ij , and the local-mode anharmonic shifts For the case of two modes (and setting ∆ 1 = ∆ 2 = ∆), the matrix representation of the above exciton Hamiltonian in the basis of the ground state as well as singly and doubly excited states becomes, Diagonalizing the matrix representation of the exciton Hamiltonian H exciton yields the vibrational energy levels E n and the corresponding eigenvectors C. Note that the manifolds of the ground state, of the singly-excited states, and of the doubly-excited states are not coupled in the exciton Hamiltonian.
Within the vibrational exciton model, the transition dipole moment operator is defined as μ = i μi , where μi is the transition dipole moment of the fundamental excitation of the i-th local mode.This can be used to set up a matrix representation of the dipole moment operator.For the case of only two modes [in analogy to (15)], one obtains, Transforming this matrix using the matrix C of eigenvectors of the exciton Hamiltonian, the transition dipole moments between the vibrational states can be obtained. 4,6 ferent strategies can be employed for obtaining the parameters of the above vibrational exciton Hamiltonian.First and most commonly, vibrational maps can be used, 34 which provide an empirical parametrization of the local-mode frequencies ωi , the harmonic coupling constants J ij , and the local-mode transition dipoles μi as functions of the molecular structure and the electrostatic environment of the local oscillators.This is commonly combined with a fixed value for the local-mode anharmonic shift ∆.
A second strategy was proposed by Hanson-Heine et al. 82 which is based on quantumchemical calculations.From a harmonic frequency calculation followed by a localization of the relevant normal modes (see Section 2.1), the localized-mode frequencies ωi and localmode transition dipoles μi can be extracted directly.The harmonic coupling constants J ij can be obtained as elements of the coupling matrix J = ( Hsub ) 1/2 , i.e., the matrix square root of the Hessian with respect to localized modes. 83These parameters are then used in combination with an empirical anharmonic shift ∆ to construct the vibrational exciton Hamiltonian.We will refer to this parametrization as "∆-exciton model".
Here, we will explore a third strategy, which exploits the connection between the L-VCISD matrix and the vibrational exciton Hamiltonian. 78To this end, we use the frequencies of the L-VSCF eigenstates (i.e., the diagonal elements of the L-VCISD matrix) on the diagonal of the vibrational exciton Hamiltonian ("L-VSCF exciton model").This way, the anharmonicity is taken into account directly and no empirical anharmonic shift is required.
The off-diagonal elements of the vibrational exciton Hamiltonian and the transition dipole moments are handled as above (i.e., they are obtained from the harmonic localized modes).

2D IR spectra
The anharmonic vibrational excitation energies ωn and transition dipole moments µ nm calculated either quantum-chemically or with the help of a vibrational exciton model can be used as input for the calculation of 2D IR spectra. 4Here, we will only consider the static case, i.e., the 2D IR spectra are calculated for a fixed molecular structure.For additional theoretical background as well as the discussion of the dynamical case, we refer to Refs.33, 91, 92.2D IR spectroscopy probes the third-order nonlinear response.In general, three IR laser pulses are used, which lead to different Feynman pathways (see Fig. 1), rephasing and non-rephasing, depending on the wave vector k = ∓k 1 ± k 2 + k 3 .The third-order response functions R n (t 1 , t 2 , t 3 ) are thus described by six double-sided Feynman diagrams that generate 2D IR spectra for any number of coupled eigenstates.The rephasing and non-rephasing response functions for any given number of eigenstates are calculated as and Here, ωn are the (anharmonic) frequencies of the vibrational excitations, with the indices i and j referring to singly excited states and the index k referring to doubly excited states.
The variables t i are the delay times between laser pulses and T 2 is the dephasing time where Ên are polarizations of the laser pulses (which are determined by the experimental setup) and μα are the transition dipole moments.The angles θ ab (with latin letters) are angles between the laser pulses, whereas θ αβ (with greek letters) are angles between transition dipole moments.The specific combination of transition dipole moments relevant for each of the six response functions is given by the subscripts of ⟨C α,β,γ,δ ⟩ in Eqs.(17)   to (22), where single indices i, j, k refer to transitions from the ground state to the specified excited state, and double indices ik and jk refer to the transition from the singly-excited state i or j to the doubly excited state k.
Since the third-order response is commonly visualized in the frequency domain, the response functions are double Fourier transformed with respect to the t 1 and t 3 coherence times, 4,22,79 S n (ω with the variable parameter t 2 determined by the delay time between the second and third laser pulse in the experiment.Now, the real parts of the response functions are summed up and double Fourier-transformed in order to get purely absorptive spectra Other evaluation methods, e.g., plotting the imaginary part or the photon echo signals that are only generated in the k I direction, are easily accessed by altering the function in Eq. ( 25), but will not be considered here.

Computational Details
All quantum-chemical calculations were performed using density-functional theory (DFT) with the Turbomole V7.5.1 program package. 94We employed the B3LYP exchangecorrelation functional 95 and the def2-TZVP basis set. 96Solvation effects were included using the COSMO model. 97,98 or the metal carbonyl complexes, we used default COSMO parameters for n-hexane, whereas for the dipeptides we used default parameters for water.
For the dipeptides, the N -deuterated structures were considered in all vibrational calculations.The quantum-chemical calculations were automated using the PyADF scripting framework. 99,100 st, the molecular structures were optimized using Turbomole and harmonic vibrational frequencies and normal modes were calculated using the SNF program. 56,85 or the localization of the relevant normal modes (consisting of the CO stretch modes for the metal carbonyls and of the amide I vibrations for the dipeptides), we employed our LocVib Python package. 56,83,101 Tis localization provides localized-mode coordinates as well as localized-mode harmonic frequencies, transition dipole moments, and harmonic coupling constants. 83e anharmonic potential energy surfaces in terms of localized modes were constructed using our Vibrations code for anharmonic theoretical vibrational spectroscopy 75,77,102 in combination with PyADF and Turbomole.Here, we used two different approximations of the potential energy surface, either including both anharmonic one-mode and twomode potentials [2-mode, see Eq. ( 7)], or including only anharmonic one-mode potentials in combination with harmonic two-mode potentials [1-mode+2h, see Eq. ( 8)].In both cases, we used grids of 16 points along each local mode coordinate as described in Ref. 75.
Subsequently, L-VSCF/L-VCI were performed with Vibrations.Here, the singly and doubly excited states were included in the VCI expansion, corresponding to an L-VCISD approximation.For further details, we refer to Refs.75, 77.
These L-VSCF/L-VCI calculations provided anharmonic vibrational frequencies and transition dipole moments for singly and doubly excited states, which could then be used as input for the calculation of the 2D IR spectra.Alternatively, parameters for the construction of a vibrational exciton model as described in Section 2.3 were extracted from the harmonic calculation and the localization as well as the L-VSCF calculation.We considered two different exciton models, termed L-VSCF exciton model and ∆-exciton model.
By diagonalizing the matrix representation of the vibrational exciton model, the same input data for the calculation of 2D IR spectra (i.e., anharmonic vibrational frequencies and transition dipole moments) can be obtained.
Using the anharmonic vibrational frequencies and transition dipole moments for singly and doubly excited states, the 2D IR spectra are calculated as described in Section 2.4.
For the dephasing time, which influences the peak broadening, we chose T 2 = 2 ps for all calculations as this is a time scale often found in amide I vibrations. 4For the waiting time t 2 we used no time delay (i.e., S(ω 3 , t 2 = 0, ω 1 )).Our code for simulating 2D IR spectra is written using Python 3 and makes use of Numpy. 103,104 ll spectra are plotted using Matplotlib. 105,106 ata set containing coordinates of all considered molecular structure, output data from the harmonic vibrational frequency calculations, the localization of normal modes and the anharmonic L-VSCF/L-VCI calculation as well as Python code and Jupyter notebooks that can be used for the calculation of 2D-IR spectra from this data and for generating the plots presented here is available at Ref. 107.

Results and Discussion
For our quantum-chemical calculations we consider four test cases that correspond to typical applications of 2D IR spectroscopy in experiment.First, we investigate the CO stretching vibrations in two metal-carbonyl complexes, dicarbonyl(acetylacetonato)rho-  The first test case is the well-studied DAR molecule, which shows a prototypical 2D IR spectrum. 22,30,108 I features two strongly-coupled anharmonic CO stretch vibrations, with an asymmetric stretch mode at ωa and a symmetric stretch mode at ωs , both of which carry strong IR intensity.These two fundamental vibrations give rise to three overtone and combination bands.
The calculated 2D IR spectra of DAR are shown in Fig. 3.The fundamentals show up as negative peaks on the diagonal, which are each accompanied by a negative peak that is due to the excitation to the corresponding overtone.The coupling between the modes results in cross-peaks, where the anharmonically shifted cross-peaks correspond to the excitation of the combination band. 30Different polarization conditions vary the intensities of the cross-peaks, which are comparably weak in the ⟨ZZZZ⟩ polarization condition, but can be strongly enhanced by using the ⟨ZZXX⟩ polarization condition. 22 our calculations of the 2D IR spectrum of DAR, we compare different computational approximations.Our most sophisticated approach uses a two-mode approximation potentialenergy surface [see Eq. ( 7)] in combination with L-VSCF/L-VCISD calculations.The corresponding results are shown in Fig. 3a and b for the ⟨ZZZZ⟩ and ⟨ZZXX⟩ polarization conditions, respectively.These spectra serve as reference for the following, more approximate models.As a first simplification, we only calculate the anharmonic one-mode potentials and approximate the two-mode potential harmonically [see Eq. ( 8)].This drastically cuts down the computational costs for the construction of the potential energy surface.The resulting 2D IR spectra are shown in Fig. 3c and d and are virtually indistinguishable from those obtained with the fully anharmonic two-mode potentials, which confirms that this approximation is valid.
As another approximation, we construct a vibrational exciton model using the anharmonic L-VSCF energy levels for the local-mode vibrations, the harmonic local-mode transition dipole moments and the harmonic coupling constants (L-VSCF exciton model).The corresponding 2D IR spectra are shown in Fig. 3e and f, which are almost identical to those from the L-VSCF/L-VCI calculations.This agreement demonstrates that an exciton model with suitable parameters is able to reproduce the results of fully anharmonic L-VSCF/L-VCI calculations, i.e., the approximations that are made in the construction of exciton models 78 are justified here.
The local-mode frequencies and coupling constants of this L-VSCF exciton models are summarized in Table I.The anharmonic L-VSCF local-mode frequencies in this model can be approximated by introducing an anharmonic shift of ∆ = 21 cm −1 on top of the harmonic local-mode frequencies, resulting in a simplified ∆-exciton model [see Eq. ( 15)].
Again, this ∆-exciton model results in 2D IR spectra (see Fig. 3g and h) that are indis- tinguishable from the other approximations.
Comparing our results to the experimental spectra measured by Khalil et al. 22 shows a good agreement with the overall line shapes and relative peak positions, but we find that our spectra are shifted to higher wavenumbers.For a more detailed comparison, For all calculations, the positions of the fundamentals are shifted compared to experiment, which is due to the shortcoming of density-functional theory as employed for the vibrational frequency calculations.However, our calculations provide a rather good agreement for the splitting between the two fundamentals, which is calculated as 61 cm −1 , while it In addition, Table II compares three different anharmonicities: the anharmonicity of the doubly-excited antisymmetric CO stretch mode ∆ 2a,a (i.e., the difference between the two peaks in the lower left corner of the 2D IR spectra), the anharmonicity of the doublyexcited symmetric CO stretch mode ∆ 2s,s (i.e., the difference between the two peaks in the upper right corner), and the anharmonicity of the the combination band ∆ as,a/s (i.e., the difference between the two peaks in the upper left and in the lower right corner).All calculations are in good agreement with the experimental anharmonicities, even though they underestimate them by a few wavenumbers.
As second example, we consider the 2D IR spectrum for the CO stretch vibrations in iron(0)-pentacarbonyl (IPC).The five local CO stretch vibrations give rise to five normal modes, of which three are IR active.Two of them are the degenerate equatorial vibrational modes ωeq , the other one is the axial vibrational mode ωax .The remaining two modes are not IR active.
Figure 4 shows the calculated 2D IR spectra for the two considered polarization condi-   A comparison of the 2D IR spectra calculated using L-VSCF/L-VCISD with full anharmonic two-mode potentials (see Fig. 4a and b) and with harmonically approximated two-mode potentials (see Fig. 4c and d) shows only minor differences, in particular for the Fermi resonance of the overtone of the axial CO stretch vibrations.This confirms that the computationally efficient approximation of using L-VSCF/L-VCI with anharmonic one-mode potentials and harmonic two-mode potentials (1-mode+2h) is applicable for the CO stretch vibrations in carbonyl complexes.Based on this approximation, it is possible to set up the L-VSCF exciton model as described above.The frequency parameters used in this model are listed in Table IV.The resulting 2D IR spectra (see Fig. 4e and f) are in very good agreement with the L-VSCF/L-VCISD results, even though the peak positions are shifted by a few wavenumbers (see also Table III).These differences are due to the approximations inherent to the vibrational exciton model compared to a full L-VCISD treatment. 78The L-VSCF exciton model can be closely approximated by using the harmonic local-mode frequencies in combination with an anharmonic shift of ∆ = 21 cm −1 (see Table IV).The corresponding 2D IR spectra for this ∆-exciton model (see Fig. 4g and h) are almost indistinguishable from those based on the L-VSCF exciton model.
Comparing our calculated 2D IR spectra to the experimental spectra recorded in solution 109 and in the gas-phase 111 shows an overall good agreement (see also Table III).Note that in experiment the peaks are broadened because of additional dynamical effects that are not included in our static calculations.As for DAR, the positions of the fundamentals are shifted, but the splitting between the two vibrational modes are reproduced quite well (28 cm −1 vs. 21 cm −1 in the experimental spectrum).The splitting of the anhamonically shifted peaks in not resolved in experiment, but the anharmonic shifts are in very good agreement.The 2D IR spectra of DKP and NAGMA calculated with L-VSCF/L-VCISD with anharmonic one-mode and harmonic two-mode potentials are shown in Fig. 5a, b and in Fig. 5c, d, respectively.The parameters of the vibrational exciton models based on our harmonic and L-VSCF calculations are listed in Table V.The 2D IR spectra calculated using L-VSCF/L-VCISD with full anharmonic two-mode potentials as well as those based on the L-VSCF and ∆-exciton models are given in Figs.S1 and S2 in the Supporting Information.As for the previous test cases, the 2D IR spectra calculated with all four approaches are almost indistinguishable.Table S2 in the Supporting Information further lists the fundamental transition frequencies as well as anharmonic shifts extracted from our calculations.
For DKP, the 2D IR spectrum shows only one negative peak corresponding to ωa on the diagonal, together with an anharmonically shifted positive peak.The calculated anharmonic shift amounts to ∆ a,2a = 4 cm −1 .For NAGMA, the calculated 2D IR spectrum shows the two fundamentals on the diagonal, along with the corresponding overtones, with an anharmonic shift of ∆ a,2a = 7 cm −1 for the more intense antisymmetric mode at lower wavenumbers, and of ∆ s,2s = 4 cm −1 for the less intense symmetric mode.On the off-diagonal, peaks corresponding to the combination band appear, with an anharmonic shift of ∆ as,a/s = 2 cm −1 .
A direct comparison to experimental spectra is not possible for DKP and NAGMA, but we can compare to the spectra calculated by Hanson-Heine et al. using a vibrational exciton model based on harmonic localized-mode frequencies and coupling constants. 82Overall, the calculated 2D IR spectra show a very good agreement.This is not surprising, since the parameters of the vibrational exciton models corresponding to our calculations (see  6 and used in all empirical amide I models.Here, we do not rely on an empirical parameter, but are able to calculate this local-mode anharmonic shift directly.From our calculations, we obtain ∆ = 4 cm −1 and ∆ = 7 cm −1 for DKP and NAGMA, respectively.These local-mode anharmonic shifts, which is significantly lower than the commonly used empirical value, and also indicate that this anharmonic shift might show significant variations.

Conclusions and Outlook
As 2D IR spectroscopy is commonly applied to probe anharmonic vibrations of coupled local oscillators, methods of anharmonic theoretical vibrational spectroscopy that are based on localized vibrational modes are ideally suited for computational 2D IR spec-troscopy.Here, we have demonstrated that the L-VSCF/L-VCI method 75,77 can be used for predicting 2D IR spectra.
We have restricted the treatment to the relevant modes (e.g., CO stretch vibrations or amide I modes) and have shown that a potential energy surface constructed using anharmonic one-mode potentials and harmonic two-mode potentials is sufficiently accurate.
This provides a computationally efficient quantum-chemical route to the relevant anharmonic energy levels.In addition to a conventional, harmonic frequency calculation (which requires one calculation of the Hessian matrix), only a fixed number of additional singlepoint energy calculations per considered mode is required.This is significantly more efficient than other quantum-chemical methods for the calculation of anharmonic energy levels, such as second-order vibrational perturbation theory, which requires at least two additional calculations of the Hessian per considered normal mode.
The approach presented here is solely based on quantum-chemical calculations and thus does not rely on any empirical parameters.If necessary, all approximations within the L-VSCF/L-VCI approach can be alleviated by (a) using a more reliable methods for the quantum-chemical calculations in the construction of the potential energy surface, (b) using a more complete representation of the potential-energy surface by including additional modes and/or by extending the number of n-mode potentials that are calculated explicitly, and (c) by extending the size of the expansion space in the L-VCI calculations, albeit each of these will increase the required computational effort. 77We believe that the efficient protocol applied here provides a reasonable starting point for the prediction of 2D IR spectra in the most common cases, as demonstrated by the examples considered here.
Because the L-VSCF/L-VCI approach is parameter-free, it will be particularly useful for the interpretation of 2D IR spectra in cases in which no reliable empirical exciton models are available, i.e., if vibrations other than the amide I region in polypeptides or the stretching vibrations in water are examined.Examples include 2D IR spectroscopy of nucleic acid bases or DNA [46][47][48] and the 2D IR spectroscopy of polypeptides beyond the amide I region, 45 which we plan to address in the future.
In this work, we only considered static 2D IR spectra that were calculated for a single molecular structure.When extending the presented approach to the dynamical case, it will become necessary to consider a large number of additional structures.As any quantum-chemical approach will quickly become infeasible for larger systems and when considering thousands of structural snapshots from molecular dynamics simulations, this usually requires the use of empirically parametrized vibrational exciton models.However, as we have demonstrated L-VSCF/L-VCI also provides access to all the parameters required for the construction of such models. 78The local-mode anharmonic energy levels can be extracted from the L-VSCF calculation, whereas the harmonic coupling constants as well as the local-mode transition dipole moments are available from the harmonic frequency calculation in combination with the localization of the normal modes. 82,83 his can be used as basis for the parametrization of vibrational exciton models, which we will also explore in our future work in order to enable the quantum-chemical prediction of 2D IR spectra of dynamical systems.

Figure 1 :
Figure 1: Double-sided Feynman diagrams for any number of coupled eigenstates.R 1 , R 2 , and R 3 are rephasing diagrams describing the third-order nonlinear response in the phase-matching k I = −k 1 + k 2 + k 3 direction, whereas R 4 , R 5 , R 6 describe the nonrephasing diagrams in the k II = +k 1 − k 2 + k 3 direction.

Figure 3 :
Figure 3: Calculated 2D IR spectra of DAR for the ⟨ZZZZ⟩ (left) and ⟨ZZXX⟩ (right) polarization conditions.(a) and (b) L-VSCF/L-VCISD calculations with anharmonic oneand two-mode potentials, (c) and (d) L-VSCF/L-VCISD calculations with anharmonic one-mode potentials and harmonic two-mode potentials, (e) and (f) anharmonic L-VSCF exciton model, (g) and (h) harmonic ∆-exciton model with an anharmonic shift of ∆ = 21 cm −1 .All plots are normalized to the intensity of the highest peak.Horizontal and verticals lines indicate the frequencies of the fundamentals.

amounts to 69 cm − 1
in experiment.The calculations of Hamm and Zanni in Ref. 4 yield a similar splitting of 59 cm −1 , while the model of Moran et al. in Ref. 79 predicts a splitting of only 27 cm −1 .

Figure 4 :
Figure 4: Calculated 2D IR spectra of IPC for the ⟨ZZZZ⟩ (left) and ⟨ZZXX⟩ (right) polarization conditions.(a) and (b) L-VSCF/L-VCISD calculations with anharmonic oneand two-mode potentials, (c) and (d) L-VSCF/L-VCISD calculations with anharmonic one-mode potentials and harmonic two-mode potentials, (e) and (f) anharmonic L-VSCF exciton model, (g) and (h) harmonic ∆-exciton model with an anharmonic shift of ∆ = 21 cm −1 .All plots are normalized to the intensity of the highest peak.Horizontal and verticals lines indicate the frequencies of the fundamentals.
. 110 for the assignment of the vibrational modes in IPC).For each of these two vibrations, the calculated 2D IR spectra show two anharmonically-shifted positive signals.An inspection of the relevant transitions in the L-VSCF/L-VCISD calculations allows us to provide an assignment of these features: For ωeq (see lower left corner of the 2D IR spectra) this is because the two degenerate vibrations of the equatorial CO stretch modes give rise to two overtones as well as one combination band, which both carry IR intensity and appear at different frequencies.For ωeq (see upper right corner of the 2D IR spectra), the overtones of the axial CO stretch vibration are close in energy with the overtone of one of the totally-symmetric CO stretch vibrations, which gives rise to a Fermi resonance.The 2D IR spectra show no visible cross-peaks.

Figure 5 :
Figure 5: Calculated 2D IR spectra for the ⟨ZZZZ⟩ (left) and ⟨ZZXX⟩ (right) polarization conditions from.(a) and (b) are for the DKP molecule, (c) and (d) are for the NAGMA molecule.
4,33,91Conceptually, each response function represents a different process.The R 1 and R 4 functions represent stimulated emission, R 2 and R 5 represent ground-state bleaching, and R 3 and R 6 represent excited-state absorption.

Table II
lists our calculated fundamental transition frequencies and anharmonic shifts and compares them to the values extracted from the experimental 2D IR spectra in Ref.22as well as to those from previous calculations in Refs.79 and 4.

Table II :
Comparison of fundamental transition frequencies for the antisymmetric and symmetric stretch modes, ωa and ωs , and of the anharmonic shifts ∆ 2a,a , ∆ as,a/s , and ∆ 2s,s for the DAR molecule.Included are the values calculated with our different approximations as well as previous experimental and computational results.All values are given in cm −1 .An extended version of this table is provided in the Supporting Information (TableS1).

Table III :
Comparison of fundamental transition frequencies for the equatorial and axial CO stretch modes, ωeq and ωax , and of the corresponding anharmonic shifts ∆ eq and ∆ ax for the IPC molecule.Included are the values calculated with our different approximations as well as previous experimental results.All values are given in cm −1 .Table III lists the fundamental frequencies as well as the anharmonic shifts found in these spectra.The degenerate modes ωeq coincide and result in peaks that are twice as intense as the ωax peaks (see also Ref

Table IV :
Local-mode frequencies and coupling constants used in the vibrational exciton model for IPC.Values obtained within the harmonic approximation are given for comparison.For the L-VSCF exciton model, the local-mode frequencies are extracted from the anharmonic L-VSCF calculation.For the ∆-exciton model, an anharmonic shift of ∆ = 21 cm −1 is applied.All values are given in cm −1 .

Table V :
Local-mode frequencies and coupling constants used in the vibrational exciton models for DKP and NAGMA.Values obtained within the harmonic approximation are given for comparison.For the L-VSCF exciton models, the local-mode frequencies are extracted from the anharmonic L-VSCF calculations.For the ∆-exciton models, an anharmonic shift of ∆ = 4 cm −1 and of ∆ = 7 cm −1 is applied for DKP and NAGMA, respectively.All values are given in cm −1 .

Table V
) are very similar to those used in Ref. 82.The main difference is in the anharmonic shifts.Hanson-Heine et al. used an empirical anharmonic shift of ∆ = 16 cm −1 for the local modes, as originally proposed by Hamm et al.