Reduced Two-Electron Interactions in Anharmonic Molecular Vibrational Calculations Involving Localized Normal Coordinates

: Spatially localized vibrational normal mode coordinates are shown to reduce the importance of calculating the full set of two-electron terms in the molecular electronic Schrödinger equation. Electron correlation and dispersion interactions become less significant in ( E,E )-1,3,5,7-octatetraene vibrational self-consistent field calculations when displacing remote atoms along multiple coordinates. Electron correlation interactions between spatially remote modes are also found to be less important compared to their corresponding uncorrelated interaction terms. Attenuation of the Coulomb operator indicates that the two-electron terms between remote electrons become less important for accurately describing the strongly contributing mode-coupling terms between sets of localized vibrational modes.


Introduction
Experimental observations are often linked to a proposed quantum mechanical system by matching the observed spectroscopic results with the state transitions that are predicted for that quantum system.This approach relies on the availability of theories that can accurately predict these spectra on human accessible time scales.
The complexity of the two-electron terms in the molecular electronic Schrödinger equation describing the system's quantum mechanical behavior rapidly increases for larger many-body systems.3][14] However, two-electron integral terms still make up the majority of the calculation, and the electron correlation energy missing from Hartree-Fock theory is known to play an important role in accurately calculating molecular vibrations.6][17][18][19][20][21][22] Increasing molecular sizes can lead to exponential increases in the dimension of both the electronic structure calculations and the nuclear PES, creating computational bottlenecks that limit accurate anharmonic calculations to systems with smaller scales.
One method for addressing scaling issues in the electronic part of the calculation is to use Kohn-Sham density functional theory (DFT) with approximate exchangecorrelation functionals. 46] Benchmarking studies using DFT-D methods have shown improvements in the accuracy of vibrational frequency calculations compared to uncorrected hybrid density functionals, even when considering small covalently bound semirigid molecules. 27In addition, dispersion interactions are known to be particularly important when calculating accurate vibrational properties in systems that contain weakly bound molecules, with Barone and co-workers having used B3LYP-D3 and dispersion corrected pseudopotential methods to calculate anharmonic frequencies of several nucleobase dimers 28 and Pitsevich et al. having looked at the anharmonic frequencies of the water dimer using second-order vibrational perturbation theory (VPT2). 29ve function based electronic correlation theories can also be made more tractable by neglecting some of the less important aspects of the electronic wave function when solving the Schrödinger equation in a basis of molecular orbitals that have been localized in physical space.Similar developments have been made to address the scaling of the PES in vibrational calculations.Unitary transformation of the normal coordinates that result from assuming that the PES is a quadratic function around a minimum energy nuclear geometry can optimize these coordinates by maximizing various measures of their spatial localization.1][32] Calculations using coordinates that are localized in this manner have been shown to remain physically meaningful while reducing the number of terms in the anharmonic PES that are needed to calculate vibrational frequencies with a given level accuracy.n both the electronic and nuclear cases, the choice of coordinates can make the mathematics describing the system more or less complex, and methods for choosing appropriate coordinates to describe a particular system are therefore important tools for increasing the computational tractability of quantum chemical calculations in general.
To date, these approaches have mainly focused on addressing computational bottlenecks in the electronic and nuclear Schrödinger equations separately; however, the electronic structure and nuclear PES calculations are not entirely disconnected.Localized vibrational coordinates often involve significantly displacing relatively few atoms compared to normal coordinates, and Cheng and Steele have previously noted that the motion of a small number of atoms within a localized coordinate describing a sufficiently large system could be expected to cause a relatively insignificant perturbation to the overall electronic structure. 331] The results of these studies have indicated that accurately describing PES electronic interactions becomes less important for both multicoordinate terms involving the displacement of remote atoms and when describing the interactions between remote atoms within coordinates. 60e focus of the current study is therefore to investigate if similar schemes can be used to reduce the computational complexity of specific two-electron terms in the electronic structure calculations that are used to calculate anharmonic vibrational frequencies without introducing significant errors into these calculations.

Theory
The most common approach for calculating nuclear vibrational frequencies is to use the harmonic approximation. 1 In this case the nuclear potential energy surface is assumed to be a quadratic function around an optimized geometry, reducing the PES evaluation to second derivatives of the energy with respect to nuclear displacements.The nuclear vibrational Schrödinger equation then factorizes into 3N-6 (or 3N-5 for a linear molecule) one-dimensional quantum harmonic oscillators, where N is the number of atoms in the system and where each of the resulting harmonic oscillators is contained along a different orthogonal normal coordinate.
Harmonic vibrational frequencies, wave functions, and normal coordinates are computed from the eigenvalues and eigenvectors of the mass-weighted Hessian matrix, (H), corresponding to the nuclear Cartesian displacements for each atom, where x i represents a Cartesian displacement coordinate of an atom with mass m i .
In anharmonic frequency calculations, the PES is often approximated as a truncated Taylor series expansion in normal mode displacement coordinates around an equilibrium molecular geometry.Derivatives of the energy with respect to nuclear displacement are then often evaluated up to the fourth order.This is referred to as the quartic force field (QFF). 62[64] The anharmonic nuclear wave function can also be approximated in a number of different ways, including using vibrational self-consistent field theory (VSCF), [65][66][67][68][69][70] which describes the nuclear wave function as one-mode wavefunctions that see an averaged potential from the other modes.VSCF mode wave functions can then be correlated using vibrational Møller-Plesset perturbation theory (VMP2), 63,[71][72][73] vibrational configuration interaction (VCI), [74][75][76][77][78] or vibrational coupled cluster (VCC) 79- 80 theories, and when VSCF is performed using localized modes, it is termed L-VSCF. 33,35 he localized modes used in L-VSCF procedures,  , have to date been those originally proposed by Jacob and Reiher 30 and are calculated using a unitary transformation, , of normal modes, , which maximizes an objective function (ξ),  = .

2
The objective function used here is analogous to the molecular orbital localization scheme of Boys and maximizes the sum of the squares of the atomic contributions to the localized modes, 31 where μ and i are the localized modes and nuclei, respectively, N is the total number of nuclei, and where R i is the position vector of each nucleus with respect to the molecular origin.
The resulting localized modes can then be used to transform the mass-weighted where  m,s   is the effective mean-field potential of mode m in state s, the subscripts new and old denote the current and previous VSCF cycles, respectively, and b is a real number that defines the extent of the mixing between the two potentials.The value of b was set to 0.9 during these calculations unless the VSCF procedure failed to converge with this mixing parameter (see the Supporting Information).
To quantify the impact of removing or reducing the two-electron terms when describing anharmonic mode coupling, the position of each mode p is defined as and the distance between two modes p and q is given by Coupling terms between modes that are centered beyond a certain cutoff distance can then be identified and modified as needed.Modifying aspects of the PES in this manner is conceptually similar to the multi-level methods used in multi level quasi-Monte Carlo 82 and recently for anharmonic molecular vibrational frequencies calculated through the single-to-all (STA)-in-VSCF and nonrelaxed (nr)-VCI-in-VSCF approaches. 83However, the distance decay of interactions between localized coordinates provides a natural basis for partitioning a multilevel approaches and further motivates the investigation into these relationships.
The minimum energy geometries predicted by two levels of theory can also differ, leading to differences in the spatial location of otherwise equivalently assigned terms in the Taylor series PES expansion.These differences are often considered to be negligible when combining anharmonic force fields that are calculated using ab initio levels of theory as these geometries tend to remain relatively convergent. 17,22,60 Hver, when extending the methodology to include methods that map from quadrature-based representations of the PES, it is worth noting that these differences may be exaggerated when using lower-resolution quadrature grids, which could in turn increase the computational complexity necessary for a given level of accuracy.

Computational Details
The geometry optimization and vibrational frequency calculations shown here were performed using the Q-Chem quantum chemical software package 84 using Hartree-Fock and second-order Møller-Plesset perturbation theory (MP2) in combination with the cc-pVDZ electronic basis set. 857] DFT calculations were performed using the global hybrid B3LYP exchange-correlation functional [88][89] and also including various DFT-D empirical dispersion corrections.The DFT-D energy corrections examined in this study include the B3LYP-D2 26 and B3LYP-D3(CSO) 66,67 methods, as these theories have been shown to give the largest absolute frequency shift and the lowest overall deviation from experiment, respectively, in a recent benchmarking study. 27DFT-D3(CSO) has been used with and without the optional addition of a three-body correction to the energy based on the Axilrod-Teller-Muto triple dipole term proposed by Grimme, 90 and the three-body corrected variation of D3(CSO) is denoted here with an asterisk as D3(CSO)*.All DFT calculations were performed using the 6-31G(d,p) basis set.Finally, the Coulomb-attenuated Schrödinger equation (CASE) approximation to HF theory, also known as attenuated HF theory, was used to test the effect of reducing the two-electron interactions within the localized and canonical normal coordinate terms. 91en calculating vibrational frequencies, harmonic calculations were initially performed at the optimized geometries, and unitary transformations were used to localize the resulting normal mode vibrational coordinates by repeated Jacobi sweeps until iteration no longer increased the objective function by more than 1×10 - 6 . 30 The normal modes were localized according to the Boys-like distance-based localization objective function shown in Eq. 3, with modes vibrating at harmonic frequencies below 1000 cm -1 , between 1000 and 3000 cm -1 , and above 3000 cm -1 localized separately.The three frequency localization windows were taken to contain modes 1-21, 22-38, and 39-48, respectively, when numbered in order of ascending vibrational energy.Anharmonic third-and fourth-order derivatives of the energy with respect to nuclear displacement were used within a Taylor series expansion of the QFF PES and were calculated using numerical differentiation of the Hessian matrices for the HF and DFT methods.Numerical energy differences were used for the RI-MP2 and CASE calculations, and when comparing with the CASE approximation the HF reference calculation has also been calculated using finite energy differences.A finite difference step size of 0.05291 Å was used throughout when calculating the anharmonic terms.Anharmonic frequencies were calculated using a 2MR representation of the QFF unless otherwise stated, and (L-)VSCF transition frequencies were calculated using a nuclear basis set consisting of the 20 lowest harmonic oscillator wave functions along each of the coordinates, with a convergence threshold set to 1×10 -4 cm -1 unless otherwise stated.The fundamental transition frequencies were taken from the eigenvalues of the ground-state onemode wave functions.Mode-coupling terms involving modes beyond a specified cutoff distance have either been removed entirely or are calculated using the alternative levels of theory stated in the text (vide infra).
Computing the 2MR derivatives formally requires evaluating the energy at a number of nuclear geometries that scales with the second power of the number of modes.
Energy, gradient, and Hessian finite difference schemes, and their combinations, have different scaling relationships that are discussed elsewhere. 64The numbers of localized mode pairs in the Hamiltonian that have distances between their mode centers that fall beyond specific cutoff distances have been shown here for selected methods to provide a measure of the increased computational efficiency that each method reduction is expected to afford.Mean absolute deviations (MADs) have been quoted at several points to aid in overall comparison between methods, and the underlying frequency data sets for all of the methods used are available in the Supporting Information.

A. Electron Correlation
4] Vibrational state energies calculated with HF/6-31G(d,p) were found to be within 0.7 cm -1 of the reference when excluding coupling terms between vibrational modes localized using the distance-based objective function and centered more than 7 a 0 apart. 33In a previously study, the current author also found similar results at the B97-1/6-31G(d,p) level of theory, where a cutoff of 4.4 a 0 maintained wavenumber accuracy. 60Larger electronic basis sets were found to become less important for remote mode-couplings, and deviations of less than a wavenumber from the 6-31G(d,p) basis set result were reported when calculating couplings between modes centered more than 2.8, 1.6, and 1.2 a 0 apart using the STO-3G, 3-21G, and 6-31G basis sets, respectively.Basis set reductions combined with modecoupling removal led to improvements over truncation alone. 60These calculations showed that a smaller basis set was needed to describe the electronic wave function when calculating couplings between more spatially remote localized modes while also retaining a high degree of overall accuracy.This study follows a similar approach to examine if reductions in the calculation of electron correlation energy and two-electron terms between remote mode-coupling terms can be used to generate similar computational improvements.
Small differences are expected between anharmonic frequency calculations using the HF and correlated RI-MP2 theories on the basis that differences in the PES described by these theories will also cause differences in the geometry and normal coordinates that they calculate.However, the frequencies that are calculated when removing the cubic and quartic terms from these calculations, shown in Figure 1, suggest that the correlated RI-MP2 calculation requires marginally more remote mode coupling in order to maintain a given level of accuracy when compared to the uncorrelated HF calculation.Calculations where the correlation energy is removed from the RI-MP2 mode-coupling terms between remote modes, and HF is used to model them instead, also show that specific levels of accuracy can be achieved while removing the correlation energy from a larger number of remote modes than can be removed altogether (see Figure 2).Subwavenumber accuracy in the zeropoint energy (ZPE) is achieved at a cutoff of 1.2 a 0 , due to error cancellation, and then from cutoff distances of 6.2 a 0 and greater.Errors of less than 10 cm -1 are achieved with cutoff distances of 0.2 a 0 , corresponding to the removal of the RI-MP2 correlation energy from ca. 92% of the mode pairs, compared with a cutoff distance of 3.8 a 0 and ca.29% of the mode pairs when removing those coupling terms entirely.Given that the electronic correlation energy plays an important role in anharmonic vibrational frequencies and that the computational complexity involved in accurately calculating the correlation energy is relatively high, this trend shows promise for reducing both the absolute computational complexity and the scaling of the complexity during these calculations.

B. Empirical Dispersion
The vibrational frequencies calculated using D3(CSO) dispersion corrections to the B3LYP DFT energy surface remain close to the frequencies calculated for the 2MR force field without dispersion.Only vibrational modes 24/25, 32, and 35/36 differ by more than ca. 2 cm -1 .The relative contributions from mode coupling are similar between all of the DFT methods tested, with subwavenumber MAD and ZPE differences from the 2MR result recovered when removing coupling terms between modes centered more than 4.4 and 6.8 a 0 apart, respectively.Larger absolute shifts of ca. 3 cm -1 per mode are seen using D2 dispersion corrections, with shifts of 5 cm -1 in the carbon-hydrogen stretching modes.The results of distance-dependent removal of the mode-coupling terms in the B3LYP-D2 potential energy surface are shown in Figure 3. Removing only the dispersion energy from the mode-coupling terms involving remote modes, shown in Figure 4, leads to a relatively small change in the ZPE of ca. 10 cm -1 at a 0.2 a 0 cutoff distance, and remains below 1 cm -1 for the larger cutoff distances tested after a cutoff distance of 4.8 a 0 .A total MAD in the fundamental transition energies of less than 1 cm -1 is also seen at a cutoff of just 0.6 a 0 .These results for the empirical dispersion corrections show a similar trend to those seen when removing the RI-MP2 electron correlation energy terms, although with a much smaller magnitude owing in part to the smaller impact that the dispersion corrections have on the L-VSCF anharmonicity of (E,E)-1,3,5,7octatetraene.

C. Coulomb Attenuation
The reduced basis set dependence noted for anharmonic vibrational frequency calculations when computed using localized rather than canonical normal coordinates suggests that the importance of accurately describing remote electronic interactions diminishes when computing certain PES terms in the localized coordinate representation. 60In this case, the contribution of two-electron interactions between remote electrons is also examined.The CASE approximation and associated attenuated HF theory provide a suitable method for testing the importance of distant two-electron terms by approximating the unattenuated theory by splitting the Coulomb operator in the electronic Schrödinger equation into a longrange and short-range piece, shown in Eq. 8, and then only preserving the shortrange two-electron interactions 91 The use of the error function in the Coulomb operator preserves the short-range twoelectron terms while smoothly and rapidly switching off the two-electron terms around the separation distance that is determined by the attenuation parameter (ω).
Applying the CASE approximation to the diagonal cubic and quartic terms in the potential energy surface causes the localized diagonal (1MR) frequencies to be more strongly affected than the corresponding set of normal mode frequencies (Table 1).
The localized modes show absolute errors of more than a wavenumber on average for an attenuation parameter of 0.25 a 0 -1 and higher, compared to 4.5 a 0 -1 for the set of normal modes.Errors of over a wavenumber are seen for individual modes at attenuation parameter values of 0.20 and 0.25 a 0 -1 , respectively.These trends are somewhat counterintuitive given that the more localized atomic displacements are expected to cause more localized perturbations to the electronic structure.However, the smaller shifts may be attributed to systematic cancellation of error when strongly displacing a greater number of more remote atoms along the more delocalized normal coordinates.
The coupled normal mode frequencies show a more intuitive pattern and are more sensitive to the attenuation of the two-electron terms than the corresponding set of localized mode frequencies (see Table 2).These trends are particularly noticeable in modes 38 and 48 in ascending energy order, which are localized onto the terminal carbon-carbon and carbon-hydrogen bonds in the local coordinate basis, and are therefore the most spatially remote modes from the perspective of other atoms in the octatetraene molecule.The errors in these modes when calculated using Coulomb attenuation, given in Table 3, show that attenuation of the two-electron terms in the set of mode couplings leads to smaller errors in the localized coordinate mode frequencies compared to the corresponding high-energy normal modes.Normal mode errors of more than a wavenumber are seen at an attenuation parameter of 0.25 a 0 -1 compared to 0.30 and 0.35 a 0 -1 for the two localized modes, respectively.Combining reductions in both the two-electron and two-mode interactions, by taking the localized mode-coupling terms in the 2MR force field between modes centered beyond a given cutoff distance and replacing them with equivalent mode-coupling terms calculated using the CASE approximation, introduces errors with respect to the full 2MR result that are shown in Figure 5 for the ZPE.There is a significant amount of variation in the accuracy of the calculated ZPEs as progressively more mode-coupling terms are attenuated to different degrees, particularly when smaller cutoff distance thresholds are used.However, the lower attenuation parameters introduce smaller errors into the overall calculation at a number of different cutoff distances and allow for more modes to be removed while retaining specific levels of accuracy.Absolute errors in the ZPEs of below 10 cm -1 are consistently achieved when using mode-coupling cutoff distances of at least 3.40, 3.00, and 2.80 a 0 for the Coulomb attenuation parameters of 0.20, 0.25, and 0.30 a 0 -1 , respectively.These cutoff distances correspond to the attenuation of ca.39, 44, and 45% of the mode pairs in octatetraene.
The different rates at which errors are introduced into the vibrational frequencies when additional mode-coupling terms are described using greater degrees of twoelectron Coulomb attenuation suggest that similar distance-based reductions in the electron correlation may be possible when using "local-correlation" methods.Several of these methods use molecular orbital localization procedures during the calculation of the correlation energy, 92 and further increases in the computational efficiency may therefore be possible through the use of properly tuned distance-based cutoffs and attenuation parameters in both the electronic and nuclear calculations.Even greater reductions in the computational complexity may also be possible through the development of colocalization procedures that localize both the normal mode and molecular orbital spaces at the same time in order to achieve an optimal balance of coordinates.A recent study has suggested that a similar combination of orbital and normal mode localizations could show potential for improving the efficiency of calculating vibrational frequencies for electronically excited states when using the restricted virtual space approximation to time-dependent DFT (RVS-TDDFT). 93wever, this is considered to be outside the scope of the current investigation.
Similarly, given the potential that two-electron attenuation has to improve efficiency in the calculation of vibrational properties, future investigations may also be warranted into the distance decay behavior of the correlated two-electron terms in attenuated MP2 theory.
Figure 5. Absolute deviations in VSCF ZPE from the 2MR HF result when calculating coupling terms between modes that lie beyond a given cutoff distance using the CASE with the various attenuation parameters (ω) given in units of a 0 -1 (top), and the number of mode pairs that are located beyond each cutoff distance (bottom).

Conclusions
Expressing the nuclear potential energy surface in terms of localized vibrational displacement coordinates has been shown to significantly reduce the importance of calculating the two-electron terms in the electronic Schrödinger equation compared to using harmonic normal coordinates when studying the (E,E)-1,3,5,7-octatetraene molecule.This discovery can be used to reduce the computationally expensive post-SCF ab initio treatment of electron correlation when describing a large portion of the PES necessary to achieve a high degree of accuracy in the molecular vibrational frequency calculations.
Removing components of the electron correlation energy from mode-coupling terms involving spatially remote displacement coordinates did not introduce significant errors into L-VSCF calculations when using RI-MP2 and DFT-D electronic structure theories.The inclusion of uncorrelated or less correlated electronic energies for remote coupling terms led to lower errors than their outright removal and allowed for the correlation energy to be removed from a larger number of coupling terms.
Improvements in the computational cost of these calculations are therefore expected.Attenuating the electronic Coulomb operator also indicated that twoelectron interactions between spatially remote electrons become less important for achieving a given level of accuracy.These findings are similar to the reduced basis set dependence previously noted for localized coordinate vibrational calculations.
Localized coordinates therefore have significant potential for improving the scaling of computational costs when calculating anharmonic force constants and molecular vibrations for larger systems, and these results indicate that a degree of redundancy is present when combining unmodified canonical molecular orbital and canonical normal mode theories.

4 L
Hessian matrix from atomic Cartesian displacement coordinates into localized mode displacement coordinates,  =  †  † .-VSCF convergence is improved by slowing down the change in the effective potential during the L-VSCF procedure by computing the effective potential at each new VSCF cycle from the weighted sum of the unaltered new potential and the effective potential obtained during the previous iteration, 81  m,s   =   m,s    + 1 −   m,s    , 5

Figure 1 .
Figure 1.Absolute deviations in VSCF ZPE from the 2MR result when removing the coupling terms between modes that lie beyond a given cutoff distance for HF and RI-MP2 PES calculations (top), together with the number of mode pairs that are located beyond each cutoff distance when using HF (bottom).

Figure 2 .
Figure 2. Absolute deviations in VSCF ZPE from the 2MR RI-MP2 result when removing the electronic correlation energy from coupling terms between modes that lie beyond a given cutoff distance (top), together with the number of mode pairs that are located beyond each cutoff distance (bottom).

Figure 3 .
Figure 3. Absolute deviations in VSCF ZPE from the 2MR B3LYP-D2 result when removing the coupling terms between modes that lie beyond a given cutoff distance (top), together with the number of mode pairs that are located beyond each cutoff distance (bottom).

Figure 4 .
Figure 4. Absolute deviations in VSCF ZPE from the 2MR B3LYP-D2 result when removing dispersion from coupling terms between modes that lie beyond a given cutoff distance (top), together with the number of mode pairs that are located beyond each cutoff distance (bottom).

Table 1 .
Mean absolute deviations (in cm -1 ) from the HF result for 1MR (L-)VSCF frequencies when calculated using Coulomb attenuated cubic and quartic terms in normal and localized coordinates.

Table 2 .
Mean absolute deviations (in cm -1 ) from the HF result for 2MR (L-)VSCF frequencies when calculated using Coulomb attenuated cubic and quartic terms in normal and localized coordinates.

Table 3 .
Absolute deviations (in cm -1 ) of modes ν 38 and ν 48 from their HF (L-)VSCF frequencies when calculated using Coulomb attenuated cubic and quartic terms.