Localized Normal Coordinates in Accurate Vibrational Structure

Jun 6, 2019 - The impact of localized CH-stretching normal coordinates in comparison to canonical normal coordinates on the performance of accurate ...
1 downloads 0 Views 978KB Size
Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. 2019, 15, 4187−4196

Localized Normal Coordinates in Accurate Vibrational Structure Calculations: Benchmarks for Small Molecules Benjamin Ziegler and Guntram Rauhut* Institute for Theoretical Chemistry, University of Stuttgart, Pfaffenwaldring 55, 70569 Stuttgart, Germany

Downloaded via NOTTINGHAM TRENT UNIV on August 13, 2019 at 08:54:08 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: The impact of localized CH-stretching normal coordinates in comparison to canonical normal coordinates on the performance of accurate vibrational structure calculations has been studied for simple molecules of up to eight atoms. Two aspects have been considered in detail, (a) the (pre)screening of coupling terms within an n-mode expansion of the multidimensional potential energy surface and (b) the demands in vibrational configuration interaction calculations (VCI). All calculations have been performed in a realistic setup, and the effect of any approximation has been measured in deviations of the final VCI frequencies, which allows for a direct comparison with experimental data.

1. INTRODUCTION Within the calculation of accurate molecular vibrational spectra beyond the harmonic approximation the choice of a meaningful coordinate system is a crucial step, which controls all subsequent calculations. Consequently, there are as many opinions as coordinate systems about the right choice, and the discussion about this is endless. In a simplified manner, one may distinguish between rectilinear coordinate systems and internal ones, the latter making use of bond lengths, angles, and dihedral angles or any linear combinations of these. It is easy to see that the latter definition of a coordinate system introduces some arbitrariness, which has extensively been discussed within the framework of geometry optimization many years ago.1 Rectilinear coordinate systems are mainly reduced to normal coordinates as arising from standard normal-mode analyses. It is their simplicity, their ease of generation, and their uniqueness which renders them attractive. For semirigid molecules, for which the harmonic approximation usually is a good one, in most cases rectilinear normal coordinates, among many others, are a proper choice as long as the potential energy surface (PES) to be considered is local. As a result vibrational frequencies can be determined to very high accuracy based on these coordinates.2 Within the last years a discussion evolved,3−13 if a unitary transformation of normal coordinates may result in preferable coordinate systems. This of course reintroduces some arbitrariness as an infinite number of unitary transformations can be defined, but the most meaningful is not known in advance. Different strategies have been presented to determine the transformation matrix automatically.4−7,13 This discussion comprises localized normal coordinates (LNC)3,5−7,14,15 including different strategies to determine subsets of coordinates to be included in the localization process,4,13 optimized coordinates (OC)8−11 and hybrid optimized local coordinates (HOLC).12 Of course, the idea of using local modes or any linear combinations of normal coordinates within vibrational structure calculations is © 2019 American Chemical Society

not new and has been pioneered many years ago by Henry, Truhlar, and others.16,17 Two ideas play a major role within these different concepts. The first one, which is the foundation of the optimized coordinates, aims at a minimization of the zero-point vibrational energy (ZPVE) at the level of vibrational selfconsistent field (VSCF) theory.8 Once this energy has been minimized, the number of correlating functions should be smaller than in standard calculations as the correlation energy is much smaller. As a consequence of that, subsequent vibration correlation methods should be much faster. However, this idea is hampered by the fact that the information about the ZPVE at the VSCF level already requires a complete PES. In principle, several routes can be chosen to overcome this problem. These rely either on the recalculation of (parts of) the PES or the rotation of it into the new coordinate system. If we assume that the PES is represented in terms of an n-mode expansion,18,19 which is one of the most common choices, a rotation of the PES typically introduces some errors. These arise from the fact that n-mode representations do not conserve the orders of the expansion, e.g., upon rotation a 3mode coupling term may become a 2-mode coupling term in the new coordinate system and vice versa.20 Once the n-mode representation of the PES is truncated fairly early, this may lead to low-order potentials of lower quality. Likewise, the use of multilevel potentials may lead to unreliable low-order terms of the PES. In summary, the idea of minimizing the ZPVE is very nice but not without obstacles. The second idea is the same one, which is used in local electron correlation methods for very many years:21 in principle, normal coordinates are delocalized over the entire molecule. A localization procedure, as for example that of Pipek and Mezey,22 can be applied to minimize the number of Received: April 19, 2019 Published: June 6, 2019 4187

DOI: 10.1021/acs.jctc.9b00381 J. Chem. Theory Comput. 2019, 15, 4187−4196

Article

Journal of Chemical Theory and Computation contributing atoms within the displacement vectors of the individual coordinates.14 This must lead to a decoupling of the contributions to the potential, which is most obvious for large molecules, e.g., when two localized coordinates include completely different sets of atoms. The advantage of this approach is the fact that the localization can be performed prior to any calculation of the PES and thus it avoids the problems discussed above altogether. However, it is by no means clear, which coordinates to localize. Previous studies have shown that a localization of all normal coordinates usually does not lead to any advantages.13 It is the introduction of subsets, which leads to improvements. However, this again introduces an arbitrariness and the assignment of coordinates to groups is not obvious at first glance. Therefore, several groups have developed techniques to generate such subsets in an automated fashion.4,6,13 Most interestingly, calculations have shown that optimized coordinates and localized coordinates, in which only the CHstretching modes and the CH-bending modes have been localized, are quite similar and due to that the localized coordinates do reduce the ZPVE in subsequent VSCF calculations as well.10,12 Moreover, the largest impact arises from the CH-stretching modes and thus the simplest definition of a set of localized coordinates relies on a localization of the CH-stretching coordinates only, while keeping all other normal coordinates unchanged. For sure, this is by no means the most sophisticated approach, but a well-defined one, which most likely retains most of the benefits. Therefore, in all what follows we will restrict the discussion to this very simple setup and we will not contribute to any strategies, which aim at an even further optimization of coordinates. For this aspect see the work of Steele et al.,13 Hanson-Heine,6 and Jacob and coworkers.4 According to the literature, the use of localized coordinates results in a couple of advantages, but it is also obvious that they introduce a number of severe problems, particularly for small molecules. A nice feature of localized coordinates is the invariance of the Watson Hamiltonian23 to any unitary transformations, i.e. all implemented programs can be used without limitations and there is no need for a specially adapted VSCF procedure etc. Concerning the advantages of localized coordinates Jacob and co-workers4 write “If the transformed coordinates {qĩ } are more localized than the normal-mode coordinates {qi} this will generally result in a faster convergence of the N-mode expansion of the potential energy surfaces, i.e. fewer high-order couplings will be required.... In turn, this can further result in a faster convergence of the VCI expansion with respect to the size of the excitation space.” These sentences nicely summarize the expectations set into localized normal coordinates. Disadvantages arise from the fact that the majority of small molecules shows some symmetry elements, which result in significant speed-ups within the calculation of the PES and vibration correlation methods. Essentially, calculating PESs of small molecules without exploiting symmetry is a waste of CPU-time. Normal coordinates belong to the irreducible representations of the molecular point group of the molecule. Upon localization this symmetry will often be destroyed but may introduce permutational symmetry at the same time. As we have shown quite recently,24 within the calculation of potential energy surfaces this situation can be handled by an algorithm, which exploits molecular point group symmetry and permutational symmetry simultaneously. Our benchmarks24

for small molecules have shown that the use of localized coordinates may lead to speed-ups within the calculation of potential energy surfaces, but this does not hold in general. However, these benchmarks did not include any (pre)screening of the individual terms of the n-mode expansion, i.e., individual terms of the n-mode expansion have not been neglected. This aspect shall be considered in this study here. Moreover, the loss of molecular point group symmetry also destroys the block structure of the Hamiltonian matrix in subsequent VCI calculations. Of course, again permutational symmetry will be introduced, but it is significantly more difficult to exploit this in order to accelerate the calculations. For example, we are not aware of an eigenvalue solver, which accounts for permutational symmetry in the Hamiltonian matrix. This may be a severe drawback of using localized coordinates but has not yet been studied in detail. The motivation for this study here is 2-fold: (1) While our recent study24 about the impact of different coordinate systems within the generation of potential energy surfaces has shown that localized normal coordinates do not necessarily lead to a reduction of CPU time, we study here the performance of different prescreening or screening schemes, which reduce the number of terms within the n-mode expansion within a given order. We are also interested in the question, if the n-mode expansion can be truncated earlier due to the use of localized normal coordinates. The second aspect (2) concerns VCI calculations. Do localized coordinates lead to faster VCI calculations due to fewer important configurations (Hartree products) in the correlation space? In order to study these effects most realistically, we did not use any model potentials or low-quality potentials but tried to reproduce experimental values to high accuracy as one would do in any application. Consequently, all our calculations presented below rely on a multilevel scheme;25,26 that is, the different orders of the nmode expansion of the PES have been evaluated at different levels of electronic structure theory, which leads to significant savings in CPU time.

2. THEORY As mentioned above, it has been shown by several authors that the number of small or even negligible high-order couplings within the n-mode expansion increases when using LNCs.4,13,27 It is thus desirable to neglect these terms prior to their explicit calculation, which requires the definition of meaningful criteria for the importance of a coupling term. Steele and co-workers suggest a distance criterion, which is based on a definition for a position vector of the individual LNCs.13,27 This criterion is in analogy to the distance criterion in local electron correlation methods, but falls short once not all coordinates have been localized.28 Therefore, we prefer a criterion which exploits already existing information about the potential. Following this strategy, a prescreening criterion for a d-dimensional coupling term can only rely on quantities computed from lower orders, i.e., (d − 1), etc. We assume that the individual coupling surfaces are given in terms of grid points. In that case the importance I of a coupling surface is related to its mean norm, i.e. 1 Ii1,..., id = d ||Vi1,..., id|| (1) n where n is the number of grid points in each direction. As the importance of a surface must be considered relatively, it must be divided by the corresponding mean norms of the underlying 4188

DOI: 10.1021/acs.jctc.9b00381 J. Chem. Theory Comput. 2019, 15, 4187−4196

Article

Journal of Chemical Theory and Computation

criteria shown in eqs 5 and 6 for prescreening (P) and screening (S), respectively.

subsurfaces. For example, the criterion, which decides about the fact, if a 3D subsurface shall be kept or not, can be estimated from Ii , j

3. COMPUTATIONAL DETAILS For all molecules, the equilibrium structure - being the expansion point of the n-mode representation of the PES - and the normal coordinates were determined at the level of explicitly correlated coupled-cluster theory in combination with an augmented basis set of triple-ζ quality, i.e. CCSD(T)F12a/cc-pVTZ-F12.31 As needed, the normal coordinates of the CH-stretching modes obtained this way were localized on the basis of Jacobi rotations as described in detail in ref 14. Within this procedure the Euclidean norm of the displacement vectors l ⃗ of a set of coordinates  is maximized according to

Ij , k

Ii , k

1 ≤ f pskip ,3 Ii + Ij Ii + Ik Ij + Ik ωiωjωk

(2)

The weight factor depending on the harmonic frequencies ω pays respect to the different nature of the individual modes. In order to generalize this empirical criterion, let us define an index vector i of length d. Moreover, let ( be the set of all partial index vectors, which arises for a given subsurface. For example, for i = (i1,i2,i3) the set ( is given as ( = {(i1 , i2); (i1 , i3); (i2 , i3)} ∪ {(i1); (i2); (i3)}

(3)

( = (2 ∪ (1

(4)

nuclei

∑∑ i∈

Note, (d is the set of index vectors with d indices. With that the generalized criterion can be written as ij yz Ij jj zz jj ∏ zz ∏ ωα−1 ≤ f pskip ,d d−2 jj zz 1 d−1 ∑ ∑ I d ′ d ′= 1 k∈( k { α∈( k j∈(

∑ |Vi*,...,i ,η| ∏ 1

η

d

α ∈{i1,..., i 2}

L1 ≤ f sskip ,d Lα

(7)

For all molecules, the n-mode expansion of the PES has been truncated after the 4-mode coupling (4D) terms. While the 1D and 2D terms have been computed at the CCSD(T)-F12a/ccpVTZ-F12 level, the explicitly correlated distinguishable clusters approach of Kats and Manby32,33 has been used for the 3D terms, i.e., DCSD-F12a/cc-pVDZ-F12. It was shown by Kats et al.33 and by Martin and co-workers34 that the distinguishable clusters approach yields very accurate geometries and harmonic vibrational frequencies of near-CCSD(T) quality. The benefit of using this coupled-cluster variant arises from the fact, that the calculation of the time-consuming (perturbative) triple-excitations can completely be avoided. Moreover, the use of explicitly correlated ab initio methods allows to employ smaller basis sets and thus the accuracy of calculations at the DCSD-F12a/cc-pVDZ-F12 level is roughly comparable to those of conventional CCSD(T)/cc-pVQZ calculations or slightly worse. Therefore, we consider DCSDF12a/cc-pVDZ-F12 theory an attractive intermediate computational level. As the number of 4D terms can be very large, we have chosen conventional MP2/cc-pVDZ calculations to limit the computational effort. All potential energy surfaces have been generated in a fully automated manner29 and symmetry has been exploited 2-fold, i.e., due to surface symmetry and within the electronic structure calculations. In the case of LNCs, permutational symmetry of the coordinates has been exploited. For the understanding of the work presented here, it is important to notice that (a) we do not use a Taylor expansion of the potential, in which the expansion coefficients can be obtained by numerical differentiation. Even more importantly, (b) the individual terms (subsurfaces) of the n-mode expansion are generated by an iterative procedure rather than a constant number of ab initio calculations. In the iterative scheme the number of ab initio points will be increased until convergence of the subsurface is achieved. This ensures that all terms of the n-mode expansions are well balanced and that the number of ab initio calculations is kept as small as possible. This feature allows for detailed insight into the performance of different coordinate systems. For details see refs 29 and 24. In contrast to that, PES generators with a fixed number of ab initio points per mode coupling surface lead to a large computational overhead and lower accuracy in case that the fixed number of grid points was chosen too small. Once a PES has been generated in terms of grid points, it has been transferred to a polynomial representation making

(5)

Of course, this criterion is purely empirical in nature, and there is no guarantee that it will work in all cases. It simply assumes that, in the case that all underlying subsurfaces are flat, the coupling surface to be investigated will also be flat. There is no proof for this assumption. However, in practice it was found that this criterion works quite well.29 Alternatively, one may screen a surface of interest at a lower computational level and may then decide if this surface should be computed or not, which has been introduced many years ago.30 Of course, this is not prescreening rather than screening and a meaningful criterion needs also to be defined. We used a criterion, which is based on the values of the vertices of a surface, V i*1,..., id , η , multiplied by a factor related to the relative extension of a surface, Lα, along a given coordinate, i.e. 1 2d

K

(|| liK⃗ ||2 )2 → max

(6)

Note, the asterisk in V* denotes the vertices of a subsurface obtained by an electronic structure calculation and L1 belongs to the mode of lowest frequency. The sum runs over the number of vertices for a given dimension. For example, for a 4D surface, this screening criterion requires 16 ab initio calculations in total and thus it might be fairly demanding to do screening instead of prescreening. However, as a qualitative answer rather than a quantitative one is sufficient at this point, we used very fast semiempirical calculations (AM1), to decide about the importance of a surface. Single-point calculations at this level of theory require just a tenth or hundredth of a second, and thus at this level, the screening time is negligible in comparison to the expense of the generation of the PES. Moreover, we found the accuracy of these low-level calculations to be sufficient to use them for estimating the importance of a surface. Our criterion can be related to the criterion of Steele and co-workers,13 who measure the coupling strength of a subsurface by the square root of the expection value of the square of the respective surface. This expectation value is evaluated from d-dimensional harmonic oscillator wave functions. Our benchmark study presented below relies on the 4189

DOI: 10.1021/acs.jctc.9b00381 J. Chem. Theory Comput. 2019, 15, 4187−4196

Article

Journal of Chemical Theory and Computation use of our fitting procedure based on Kronecker products.35 Note that this step essentially needs no time. Subsequently, VSCF modals have been determined from a basis of 20 modedependent distributed Gaussians. Corrections due to vibrational angular momentum terms (VAM) were added to the energy a posteriori.36 In the subsequent vibrational configuration interaction calculations, ground-state based modals were employed throughout. The VCI calculations included up to quintuple excitations. The VCI space was further restricted by single-mode excitations up to the sixth root and a maximal sum of quantum numbers of 12. In case of the 8-atomic molecules of the test set this scheme resulted in about 23 million configurations (Hartree products). From this initial configuration space, significant configurations have been selected in an iterative manner using a perturbational selection criterion.37−40 As the computational cost of these calculations is mainly controlled by the number of selected configurations, the loss of the block structure of the VCI matrix due to localized coordinates is less important in these configurationselective approaches. Likewise it is not so much a question about the order of the expansion of the correlation space, i.e., quadruple or quintuple excitations, rather than the individual number of selected configurations until convergence of a vibrational state is achieved. Consequently, a configurationselective algorithm allows for more insight with respect to the performance of a VCI calculation in dependence on the coordinate system. Eigenpairs have been determined in an iterative fashion using our residuum based algorithm (RACE).41 VAM terms were included in all VCI calculations, but the expansion of the μ-tensor was truncated after the constant term.36 This approximation is known to result in very small errors for semirigid systems of the size considered here.

Table 1. Maximum Correction of Fundamental VCI Frequencies Due to the Inclusion of 4-Mode Coupling Terms within the Potential Function and Deviations with Respect to Experimental Results NC molecule

MAD

CH3F C2H4 C2H3F C2H6 C2H5F

2.6 2.7 3.0 1.5 2.8

a

b

LNC c

MAX

3D→4D

MAD

6.4 5.9 5.9 3.6 8.4

16.5 31.6 15.3 31.5 21.0

1.8 1.8 1.4 1.9 2.3d

a

MAXb

3D→4Dc

3.6 3.3 4.3 6.7 9.0d

3.3 11.9 20.2 41.8 40.5

a

Mean absolute deviation of the fundamental 4D VCI frequencies with respect to experimental results. bMaximum deviation of the fundamental 4D VCI frequencies with respect to experimental results. c Maximum correction of the 3D VCI frequencies due to the inclusion of 4D terms. dThe lowest torsional mode has been neglected, see text.

the observations of other authors. The same appears to be true for maximum deviations, but for C2H5F the lowest torsional mode shows a significantly larger deviation of 16.8 cm−1 once LNCs are employed. When including this outlier into the statistics, the MAD for C2H5F is 3.1 cm−1. The reason for the poor performance of this particular mode is not yet clear, but it is hardly of any relevance for the results of this study. In any case, 4-mode coupling terms cannot be neglected once high accuracy is needed, neither for NCs nor for LNCs. In all cases the maximal correction, 3D → 4D, shows up for the CHstretching modes. For CH3F and C2H4 the LNCs yield smaller maximum corrections than the NCs, but for all other molecules the opposite is true. As the magnitude of the maximal correction is difficult to predict in advance, we conclude here that the neglect of entire orders within the nmode representation will hardly be possible. Therefore, computational benefits within the generation of PESs arising from LNCs need to be due to an increased sparsity of the coupling terms, but not due to an earlier truncation of the nmode expansion. As an example, in Table 2 we have shown the results for the fundamental modes of the ethane molecule in dependence on the order of the n-mode expansion and on the type of the coordinates. Corresponding tables for the other molecules of the test set are provided in the Supporting Information. Clearly, the 4D results for both sets of coordinates are in excellent agreement with the experimental results. Moreover, the CH-stretching modes show the largest deviations at the 3D level, which is a well-known result. It can be seen from this table that our calculations result in symmetry-broken solutions for the degenerate modes. This results from the incompleteness of the n-mode expansion and thus this effect is much smaller for the 4D calculations than for the 3D calculations. Large symmetry breakings indicate that n-mode couplings of even higher order may be mandatory. In Table 3 we have listed the number of ab initio calculations needed for generating the potential energy surfaces without neglecting any terms. In several cases the LNCs perform worse than the NCs. We mainly explain this with an increase of the complexity of the individual subsurfaces. Once the structure of a given subsurface becomes more challenging for the fitting routines, more iterations are needed to obtain converged results within the surface build-up. This has significant impact for high-order coupling terms. Let us assume a particular 4-mode coupling surface without any symmetry requests 3 iteration steps for

4. DISCUSSION AND RESULTS As outlined above, we have restricted the localization to the CH-stretching coordinates only and the focus of this work here is not on a more sophisticated grouping of coordinates or constrained localization techniques to yield even further reductions in CPU-time. As a consequence of that, we have chosen a benchmark set of small molecules with an increasing number of hydrogen atoms and of different molecular point groups, namely methylfluoride (CH3F, C3v), fluoroethene (C2H3F, Cs), ethene (C2H4, D2h), fluoroethane (C2H5F, Cs) and ethane (C2H6, D3d). 4.1. Reference Calculations. The sparsity of the terms within the individual orders of the n-mode expansion of the PES does not provide any information about the importance of these orders. As long as some terms within a given order are large, these terms cannot be neglected. Therefore, prior to any (pre)screening we have determined the importance of the 4mode coupling terms by simply calculating the maximum shift in the fundamental VCI frequencies obtained from 3D and 4D calculations for the molecules of our test set. Besides this we list here the mean absolute (MAD) and maximum absolute deviation (MAX) with respect to experimental results. The values are summarized in Table 1. Note that in all VCI calculations, on which the data in Table 1 do rely, a large correlation space including up to 5-tuple excitations has been employed, so that effects arising from an incomplete basis must be considered to be very small. Table 1 shows that the overall agreement with experimental data is excellent. On average LNCs provide slightly better MADs with respect to experimental data than NCs. This result is in agreement with 4190

DOI: 10.1021/acs.jctc.9b00381 J. Chem. Theory Comput. 2019, 15, 4187−4196

Article

Journal of Chemical Theory and Computation

Table 2. VCI Results for the Fundamental Modes of Ethane in Dependence on the Order of the PES Expansion and the Type of the Coordinates Exp.a

NC 3D

ΔExp.

NC 4D

ΔExp.

LNC 3D

ΔExp.

LNC 4D

ΔExp.

2985.4 2985.4 2968.7 2968.7 2897.8b 2895.7 1472.0 1472.0 1468.1 1468.1 1397 1379.2 1195.3 1195.3 994.1 821.7 821.7 289.3 MAD MAX

2953.3 2957.3 2940.1 2956.1 2890.0 2891.9 1470.3 1470.8 1467.3 1467.3 1396.1 1377.3 1200.2 1200.6 995.7 821.1 820.6 302.8

(32.1) (28.1) (28.6) (12.6) (7.8) (3.8) (1.7) (1.2) (0.8) (0.8) (0.9) (1.9) (4.9) (5.3) (1.5) (0.6) (1.2) (13.5) 8.2 32.1

2984.8 2984.7 2968.8 2967.4 2894.2 2896.0 1468.5 1469.1 1468.9 1466.7 1395.5 1376.8 1195.8 1195.7 993.7 818.3 818.3 289.6

(0.6) (0.7) (0.1) (1.3) (3.6) (0.3) (3.6) (2.9) (0.8) (1.4) (1.5) (2.4) (0.5) (0.4) (0.4) (3.4) (3.4) (0.3) 1.5 3.6

3018.5 3006.9 2985.8 2950.4 2940.7 2908.1 1466.9 1467.2 1459.1 1459.9 1389.7 1371.1 1189.0 1188.5 994.7 811.8 811.0 283.7

(33.1) (21.5) (17.1) (18.3) (42.9) (12.4) (5.1) (4.9) (9.0) (8.2) (7.3) (8.0) (6.3) (6.8) (0.6) (9.9) (10.7) (5.6) 12.7 42.9

2982.5 2981.2 2967.3 2967.7 2898.9 2893.1 1470.6 1471.1 1470.4 1468.0 1396.9 1378.3 1197.8 1197.6 992.4 820.7 820.0 296.0

(2.9) (4.2) (1.4) (1.0) (1.1) (2.6) (1.4) (0.9) (2.3) (0.1) (0.1) (0.9) (2.5) (2.3) (1.7) (1.1) (1.7) (6.7) 1.9 6.7

Experimental data taken from ref 42. bThis band has been assigned to the overtone of the fundamental occurring at 1472 cm−1 in ref 42.

a

on PESs spanned by NCs or LNCs, respectively. In these calculations, the number of neglected coupling terms within the n-mode expansion of the PESs was steadily increased according to our (pre)screening criteria, i.e., the criteria have been computed for each subsurfaces of a given order and have been sorted by magnitude. Subsequently, a user-defined number of subsurfaces with the lowest impact (the lowest values according to our criteria) have been omitted. Consequently, the number of neglected coupling terms may i 3N − 6 y vary between 0 and jjjj d zzzz. In the latter case, the result must k { coincide with the result for the corresponding (d−1) calculation. We have performed these calculations independently for the 3D and 4D terms of the potential. In summary, for each fundamental mode of a molecule we have performed about 200 VCI calculations to analyze the dependence of the VCI fundamental frequencies on (a) the (pre)screening criterion, (b) the expansion of the PES, and (c) the use of LNCs vs NCs. In Figure 1 we have shown the VCI error for the individual fundamentals of C2H4 relative to a VCI calculation, which was based on a full potential, i.e., without neglecting individual surfaces. In this figure each fundamental is represented by an individually colored line. CH-stretching modes are colored in blue. As must be expected, the error increases with a growing number of neglected surfaces. Let us accept a maximum error of 1 cm−1 in the final VCI results. It can then be seen from Figure 1 that about 20 3D terms can be neglected by prescreeing for the NC potential and about 22 for the LNC potential. In general, the errors for the LNC PESs look less erratic than for the NC PESs. Note that, the maximum errors usually arise from the CH-stretching modes, which, as discussed above, show the largest correction arising from high-order couplings. Screening rather than prescreening leads to more systematic results and many more 3D coupling terms can be neglected for the LNC potential in comparison to the NC potential. Switching to the 4D potentials essentially

Table 3. Number of Terms and Energy Single-Point Calculations (Rounded to Thousands) for the Molecules of the Test Set in Dependence on the Order of the Expansion and the Type of the Coordinates terms

3D Points

4D Points

molecule

3D

4D

NC

LNC

NC

LNC

CH3F C2H4 C2H3F C2H6 C2H5F

84 220 220 816 816

126 495 495 3060 3060

17 000 11 000 56 000 70 000 114 000

24 000 10 000 64 000 175 000 192 000

87 000 55 000 285 000 673 000 956 000

54 000 47 000 162 000 793 000 908 000

NCs, but 4 for LNCs. In that case, 2800 more ab initio calculation would be required, which is a large number. It thus appears that the number of electronic structure calculations is mainly dominated by a large number of iterations for a moderate number of subsurfaces, rather than the large number of coupling surfaces. Of course, this effect cannot be seen, when a constant number of ab initio calculations per subsurface is used. However, this should be avoided, because the iterative build-up of the subsurfaces leads to tremendous savings in CPU time. On average, there appears to be a vague trend that LNC calculations require less 4D points, but more 3D points. Once multilevel schemes are used, in which the 3D single point calculations are more expensive than the 4D points, this may result in more demanding calculations once LNCs are used. According to Table 3 the average number of ab initio calculations per dimension within the 3D terms scatters between 3.5 and 6.5, while 3.1 up to 5.1 points are used within the 4-mode coupling terms. These numbers show that the effect due to the neglect of individual subsurfaces can be easily overcompensated by an additional step within the iterative build-up of the respective coupling terms. 4.2. (Pre)screening of Coupling Terms. For the smaller molecules of the test set, in particular for CH3F, C2H3F, and C2H4, we have performed systematic VCI calculations relying 4191

DOI: 10.1021/acs.jctc.9b00381 J. Chem. Theory Comput. 2019, 15, 4187−4196

Article

Journal of Chemical Theory and Computation

Figure 1. Deviation of the VCI fundamental modes of ethene due to the neglect of an increasing number of coupling terms in dependence on the type of (pre)screening and the coordinate system.

terms. The percentage of negligible 4D terms scatters between 13.2% and 49.2%, which indicates that (pre)screening starts to be efficient. However, the number of neglected ab initio single point calculations deviates significantly from the number of neglected coupling terms. This result is more pronounced for NCs than for LNCs. Clearly, those coupling terms will be neglected, which can be described by a rather low number of ab initio calculations. These are fairly smooth surfaces of low complexity. This is important to notice, because the counting of negligible subsurfaces as a measure for computational savings apparently is misleading once a subsurface can be represented by a variable number of ab initio points, which we consider to be most meaningful. An average ratio (neglected terms vs neglected points in %) of 2.3 can be observed for NCs, while a lower value of 1.4 can be seen for LNCs. This result again supports our argument of an increased complexity of the coupling terms based on LNCs, i.e. due to an pronounced complexity the number of ab initio points per surface increases and thus the neglect of such surfaces leads to the neglect of more ab initio points than in the case on NCs. As a consequence of that, the percentage of neglected ab initio calculations is lower and scatters between 5.7% and 35.9% only. A neglect of about one-third of all 4D ab initio calculations (LNCs) looks rather promising and suggests large computational savings as the number of 4D ab initio calculations is very large. And indeed this argument holds true once a single electronic structure level has been employed for all n-mode coupling terms. Yet, the neglect of one-third of all ab initio calculations in contrast to about 10% in case of NCs does not take into account, that the number of electronic structure

shows very similar results but allows for the neglect of many more surfaces than within the 3D terms. Figure 1 allows for several conclusions: (1) screening rather than prescreening leads to a more reliable selection of negligible coupling terms, (2) many more 4D terms can be neglected than 3D terms, and (3) LNCs perform slightly better than NCs and show a less erratic behavior. This analysis allowed to retrieve thresholds for the criteria of eq 5 and 6, which are listed in Table 4. We have transferred Table 4. Prescreening and Screening Thresholds for Neglecting d-Dimensional Coupling Terms in an n-Mode Expansion of the PESa f skip p,d skip fs,d

d=3

d=4

10−17 10−2

10−26 10−1

a

All values in a.u.

these thresholds to all molecules of the small test set and then performed VCI calculations. The mean and absolute deviations of these calculations with respect to the reference calculations without the neglect of coupling terms are shown in Tables 5 and 6 for C2H6 and C2H5F, respectively. In general, these tables support the findings obtained from Figure 1 (see above). Moreover, the percentage of 3D terms or 3D ab initio points, which can be neglected, is below 10% in all cases and thus it is worthwhile to ask, if these rather limited savings are large enough in order to justify the uncertainty, which has been introduced by (pre)screening. Note that these deviations were found to be as large as 3.5 cm−1. The situation alters for the 4D 4192

DOI: 10.1021/acs.jctc.9b00381 J. Chem. Theory Comput. 2019, 15, 4187−4196

Article

Journal of Chemical Theory and Computation

Table 5. Percentage of Neglected 3D and/or 4D Terms/Points within the PES Due to Screening (S) or Prescreening (P) in Dependence on the Coordinate Systems for C2H6 terms a

order NC

3D 4D 3D+4D

LNC

3D 4D 3D+4D

crit.

3D/%

P S P S P S P S P S P S

2.1 7.0

MADb

points 4D/%

3D/%

4D/%

0.6 2.0 19.8 23.0 23.1 26.3

2.1 7.0 6.0 9.6

38.3 45.8 42.0 49.2

6.0 9.6

8.3 9.7 9.5 11.2

0.6 2.0 1.4 1.8

28.4 33.6 30.9 35.9

1.4 1.8

−1

cm

0.2 0.2 0.4 0.4 0.6 0.5 0.6 0.6 0.8 0.7 0.8 0.8

MAXc cm−1 0.7 0.8 1.7 2.6 2.3 2.6 3.5 1.3 2.7 3.2 2.5 3.3

a Order of the PES expansion, which has been used for (pre)screening, nD refers to an PES expansion up to nD terms and (pre)screening has been applied to the highest order only. 3D+4D refers to a 4D calculation, in which (pre)screening was applied to both, 3D and 4D terms. bMean absolute deviation (in cm−1) of the fundamental VCI frequencies with respect to reference calculations without neglecting individual subsurfaces. c Maximum absolute deviation (in cm−1) of the fundamental VCI frequencies with respect to reference calculations without neglecting individual subsurfaces.

Table 6. Percentage of Neglected 3D and/or 4D Terms/Points within the PES Due to Screening (S) or Prescreening (P) in Dependence on the Coordinate Systems for C2H5F terms ordera NC

3D 4D 3D+4D

LNC

3D 4D 3D+4D

crit.

3D/%

P S P S P S P S P S P S

0.7 2.8

0.7 2.8 5.0 6.3

5.0 6.3

points 4D/%

3D/%

4D/%

0.2 0.6 13.2 25.2 14.4 27.7

30.6 46.3 33.3 48.6

0.2 0.6 1.2 1.3

1.2 1.3

5.7 11.2 6.2 12.2

23.5 33.0 25.6 34.7

MADb

MAXc

cm−1

cm−1

0.1 0.1 0.1 0.5 0.1 0.4 0.1 0.1 0.4 0.4 0.3 0.5

0.9 0.9 0.6 1.7 0.7 1.3 0.8 0.9 1.5 1.6 1.1 2.4

a

Order of the PES expansion, which has been used for (pre)screening, nD refers to an PES expansion up to nD terms and (pre)screening has been applied to the highest order only. 3D+4D refers to a 4D calculation, in which (pre)screening was applied to both, 3D and 4D terms. bMean absolute deviation (in cm−1) of the fundamental VCI frequencies with respect to reference calculations without neglecting individual subsurfaces. c Maximum absolute deviation (in cm−1) of the fundamental VCI frequencies with respect to reference calculation without neglecting individual subsurfaces.

computationally cheaper than the 3D terms. Screening diminishes the CPU times for the 4D terms to 194 min. (NC) and 154 min. (LNC). Putting these 4D savings in relation to the CPU time for the 3D terms, this corresponds to 4% for NCs and 8% for LNCs. These again are very small numbers and hardly justify the uncertainty introduced into the final VCI results. In other words: the neglect of high-order terms due to (pre)screening affects the cheapest electronic structure calculations and thus very high percentages of negligibe coupling-terms are necessary in order to obtain considerable CPU time savings once multilevel schemes are employed, which we strongly recommend and which corresponds to modern local electron correlation methods. In summary, with respect to CPU time savings we do not see any advantages due to (pre)screening, even when LNCs are used. This statement of course is valid only for small systems with probably less than 10 atoms. However, another aspect

calculations with high symmetry will be larger in NC calculations than in LNC calculations. As a consequence of that the advantage of LNC calculations with respect to CPU time diminishes. This aspect has been discussed in detail in ref 24. Moreover, most accurate calculations for larger systems use multilevel schemes, which introduce only very modest errors with respect to one-level calculations once the different levels have been chosen carefully. As outlined above, we have chosen DCSD-F12a/cc-pVDZ-F12 theory for the 3D terms and MP2/ cc-pVDZ theory for the 4D terms and the final results are in excellent agreement with experimental results as shown in Tables 1 and 2. According to this the CPU time needed for the 3D terms of C2H6 is 470 min. (NC) and 1183 min. (LNC), but 212 min. (NC) and 249 min. (LNC) for the 4D terms only. These numbers refer to calculations on 64 cores and do not use (pre)screening at all. We note here that it is a typical result of multilevel calculations that the 4D terms are 4193

DOI: 10.1021/acs.jctc.9b00381 J. Chem. Theory Comput. 2019, 15, 4187−4196

Article

Journal of Chemical Theory and Computation

Table 7. Average (MEAN) and Maximum (MAX) Number of Configurations (Hartree Products, Rounded to Hundreds) Needed in the VCI Calculations of the Fundamental Modes for the Molecules of the Test Set NC 3D

LNC 3D

NC 4D

LNC 4D

molecule

MEAN

MAX

MEAN

MAX

MEAN

MAX

MEAN

MAX

CH3F C2H4 C2H3F C2H6 C2H5F

2000 2500 8600 6000 13 400

3200 5000 20 100 12 900 37 300

2900 7400 10 600 23 300 26 400

4700 12 600 25 400 61 400 70 000

2800 4600 14 600 12 700 25 900

4500 9200 38 500 29 400 78 300

3200 8900 12 700 22 800 26 200

5200 16 600 33 200 51 500 67 200

is based on a perturbational criterion and is performed in an iterative manner. This leads to results very close to the full-VCI limit once the initial space has been chosen properly. The initial set chosen here comprises up to quintuple excitations, which results in 2.3 × 107 configurations for C2H6. The average and maximum numbers of configurations determined for the fundamental modes of our test molecules are listed in Table 7. Usually, many more configurations are needed for a proper description of CH-stretching modes than for all others. Consequently, the maximum numbers presented in Table 7 refer to CH-stretchings and the average numbers depend very strongly on the number of hydrogen atoms within the individual molecules. It can be seen from Table 7 that molecules with high orders of molecular symmetry request fewer configurations than those with low orders, irrespective if NCs or LNCs are used. Moreover, in almost all cases (an exception is C2H3F (4D)) LNC calculations require more configurations than NCs. This appears to be more pronounced for highly symmetric molecules. As the number of selected configurations very much controls the efficiency of the VCI calculations, we conclude that for molecules with high symmetry LNCs will lead to more demanding VCI calculations, even for algorithms, which have been optimized for permutational symmetry. Note that the conclusions drawn above do not allow for any statements concerning the order of the excitation level within the expansion of the VCI space. All calculations presented above have been obtained from a correlation space up to quintuple excitations. The question arises, if LNCs allow to reduce this fairly large space to quadruple excitations without significant loss in accuracy. Therefore, we have performed respective calculations for all systems of the test set. For the small molecules, i.e., CH3F, C2H4, and C2H3F, the changes are indeed very small and their is no need for quintuple excitations once LNC have been used. However, the same holds true for the NC calculations, which simply shows that the correlation space is converged in both cases. However, for the larger systems, the situation changes. For example, VCI calculations including up to quadruple excitations for C2H6 led to mean abolute deviations of 9.2 cm−1 (NC) and 3.6 cm−1 (LNC) and to maximum deviations with respect to experimental results of 20.9 cm−1 (NC) and 22.8 cm−1 (LNC). This example shows that the overall performance of the LNCs indeed tends to be better, but the occurrence of outliers, which lead to large maximum deviations, prohibits the restriction of the correlation space. As mentioned above, this aspect is of limited impact for configuration selective approaches, but of significance for conventional calculations. The outliers observed in these calculations are resonating CH stretching modes with small leading VCI coefficients. In this case, the high order excitations may have significant impact on the sensitive balance on the respective resonances and can thus not

concerns memory savings. Considering again C2H6 the storage of all 4D surfaces in memory may require 3.6 GB of memory once the surfaces are stored on a fine grid of 20 points per degree of freedom, which can be obtained from interpolation. Neglecting now up to 50% of all 4D terms (rather than ab initio points), leads to memory savings of up to 1.8 GB, which is much more significant than the CPU time savings. Therefore, we conclude that the use of LNCs might be preferable when multicore architectures with limited amounts of memory are used for computing the PES. Within the discussion of advantages or disadvantages the question arises, if the benefits of localization are largest, if the symmetry of the molecule is very high. If this would be true, one would see pronounced differences for small molecules, while for larger systems the effects should be smaller. This prompted us to include pairs of related systems into the test set, i.e., C2H4/C2H3F and C2 H6/C2 H5F, which differ significantly in the order of the point group. However, none of the data provided above support this statement, and thus we conclude that such effects must be small if they exist at all. 4.3. Analysis of VCI Calculations. The localization of normal coordinates inevitably destroys the block structure of the VCI matrix and thus LNCs must lead to significantly longer computation times in conventional VCI calculations and - perhaps more importantly - to larger memory demands. E.g. for ethene the destruction of the block structure may lead to an increase of the memory requirements by about a factor 64 once sparsity is not exploited. Of course, conventional VCI calculations are hardly performed nowadays and more sophisticated approaches diminish the bottlenecks of such calculations, e.g., configuration-selective and/or tensor-decomposition based approaches. However, the exploitation of permutational symmetry, which has possibly been introduced by LNCs, is much more difficult to implement into these efficiency-driven methods. This makes a direct comparison of CPU times rather difficult, as one would compare an optimized algorithm for molecular symmetries with a preliminary one for permutational symmetry. Therefore, we have switched off symmetry in our VCI benchmark calculations, being well aware of the fact that symmetry will accelerate the NC calculations. For C2H4 and C2H6 the LNC calculations for all fundamental transitions took 4 times as long as the NC calculations, i.e., 18.7 min (C2H4, NC) and 75.0 min (C2H4, LNC) and 46.5 h (C2H6, NC) and 185.4 h (C2H6, LNC). However, for C2H3F we observed a moderate advantage for the LNC calculation, i.e. 4.7 h (NC) and 3.5 h (LNC). In order to understand why the VCI calculations in most cases took much longer for LNC potentials than for NC potentials, we counted the configurations, which have been selected by our algorithm to be important. Our selection algorithm has been presented in detail elsewhere,37,38 but the principle idea is to select significant configurations from a large initial set. The selection 4194

DOI: 10.1021/acs.jctc.9b00381 J. Chem. Theory Comput. 2019, 15, 4187−4196

Journal of Chemical Theory and Computation



be neglected. Similar effects were found for C2H5F. Therefore, we conclude that the restriction of the order of the VCI correlation space in case of LNCs might be possible for wellbehaved modes, but certainly is not feasible for multiconfigurational modes, which become more prevalent in larger molecules and in high-frequency regions.

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.9b00381. Tables containing VCI results for the fundamental modes of methylfluoride, fluoroethene, and ethene (PDF)

5. SUMMARY AND CONCLUSIONS For a small test set of simple molecules with up to eight atoms, systematic benchmark calculations have been performed to study the impact of a localization of the CH-stretching modes on the performance within the construction of an accurate potential energy surface and subsequent vibrational configuration interaction calculations. The focus of this work here is on the (pre)screening of n-mode coupling terms. As must be expected, screening was found to be more reliable than prescreening and both techniques are more efficient for LNCs than for NCs. However, this effect is diminished, when counting the number of individual ab initio calculations in a realistic setup with a variable number of ab initio grid points within the different terms of the PES expansion. Moreover, the use of multilevel schemes, which is a widely used and very efficient approximation within the construction of highdimensional potential energy surfaces, essentially eliminates any advantages due to (pre)screening once high accuracy is requested. Besides that, an earlier truncation of the n-mode expansion in LNCs without loss of accuracy is not supported by this study. Configuration-selective VCI calculations relying on PESs spanned by LNCs suffer from two aspects, (a) the loss of the block structure of the Hamiltonian matrix and (b) a larger number of significant configurations in comparison to VCI calculations relying on canonical normal coordinates. As VCI calculations close to the full-VCI limit can be very timeconsuming for larger systems, this effect may not be underestimated. Therefore, configuration selective VCI calculations based on NCs tend to outperform those based on LNC potentials. Moreover, sensitive resonances in individual modes prohibit a restriction of the correlation space to lower orders of its expansion. For small molecules, as discussed here, the aspect of memory demands may be more important than a consideration of CPU time. Apparently, the improved neglect of potential terms in the case of LNCs leads to a decrease of memory requirements within the construction of the potential. Contrarily, the loss of the block structure in VCI calculations leads to higher memory demands and thus the individual memory requirements depend very much on the employed coordinate system. In summary we conclude that for small molecules, which show at least some symmetry elements, the use of localized CH-stretching modes does not bear significant advantages, but this result of course depends on the definition of the unitary transformation matrix. We address this outcome to the efficiency of multilevel schemes within the construction of PESs and mainly to the loss of the block structure in the VCI matrix, accompanied by an increased number of important configurations in case of LNCs. For sure the situation will alter for larger systems, for which high-level electronic structure calculations and thus a steep hierarchy in the multilevel approach are prohibitive and which usually do not show any symmetry elements. However, the crossover point is not yet known, but most likely will not show up for systems of less than 10−12 atoms.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +49 (0) 711 685-64405. Fax: +49 (0)711 685-64442. ORCID

Guntram Rauhut: 0000-0003-0623-3254 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS Support was provided by the COST Action CM1405 ”Molecules in Motion (MOLIM)”. This research was also supported in part by the bwHPC initiative and the bwHPC-C5 project provided through associated compute services of the JUSTUS HPC facility at the University of Ulm. bwHPC and bwHPC-C5 are funded by the Ministry of Science, Research and the Arts Baden-Württemberg (MWK). Further financial support was provided by the Deutsche Forschungsgemeinschaft (DFG), project RA 656/20-1, and the Studienstiftung des Deutschen Volkes.



REFERENCES

(1) Eckert, F.; Pulay, P.; Werner, H.-J. Ab initio geometry optimization for large molecules. J. Comput. Chem. 1997, 18, 1473− 1483. (2) Ziegler, B.; Rauhut, G. Accurate vibrational configuration interaction calculations on diborane and its isotopologues. J. Phys. Chem. A 2019, 123, 3367−3373. (3) Panek, P. T.; Jacob, C. R. On the benefits of localized modes in anharmonic vibrational calculations for small molecules. J. Chem. Phys. 2016, 144, 164111. (4) Panek, P. T.; Hoeske, A. A.; Jacob, C. R. On the choice of coordinates in anharmonic theoretical vibrational spectroscopy: Harmonic vs. anharmonic coupling in vibrational configuration interaction. J. Chem. Phys. 2019, 150, 054107. (5) Hanson-Heine, M. W. D. Examining the impact of harmonic correlation on vibrational frequencies calculated in localized coordinates. J. Chem. Phys. 2015, 143, 164104. (6) Hanson-Heine, M. W. D. Intermediate vibrational coordinate localization with harmonic coupling constraints. J. Chem. Phys. 2016, 144, 204116. (7) Hanson-Heine, M. W. D. Reduced basis set dependence in anharmonic frequency calculations involving localized coordinates. J. Chem. Theory Comput. 2018, 14, 1277−1285. (8) Yagi, K.; Keceli, M.; Hirata, S. Optimized coordinates for anharmonic vibrational structure theories. J. Chem. Phys. 2012, 137, 204118. (9) Yagi, K.; Otaki, H. Vibrational quasi-degenerate perturbation theory with optimized coordinates: Applications to ethylene and trans-1,3-butadiene. J. Chem. Phys. 2014, 140, 084113. (10) Thomsen, B.; Yagi, K.; Christiansen, O. Optimized coordinates in vibrational coupled cluster calculations. J. Chem. Phys. 2014, 140, 154102. (11) Thomsen, B.; Yagi, K.; Christiansen, O. A simple state-average procedure determining optimal coordinates for anharmonic vibrational calculations. Chem. Phys. Lett. 2014, 610−611, 288−297.

4195

DOI: 10.1021/acs.jctc.9b00381 J. Chem. Theory Comput. 2019, 15, 4187−4196

Article

Journal of Chemical Theory and Computation (12) Klinting, E. L.; König, C.; Christiansen, O. Hybrid Optimized and Localized Vibrational Coordinates. J. Phys. Chem. A 2015, 119, 11007−11021. (13) Cheng, X.; Talbot, J. J.; Steele, R. P. Tuning vibrational mode localization with frequency windowing. J. Chem. Phys. 2016, 145, 124112. (14) Rauhut, G. Configuration selection as a route towards efficient vibrational configuration interaction calculations. J. Chem. Phys. 2007, 127, 184109. (15) Jacob, C. R.; Luber, S.; Reiher, M. Analysis of Secondary Structure Effects on the IR and Raman Spectra of Polypeptides in Terms of Localized Vibrations. J. Phys. Chem. B 2009, 113, 6558− 6573. (16) Henry, B. R. Use of local modes in the description of highly vibrationally excited molecules. Acc. Chem. Res. 1977, 10, 207−213. (17) Thompson, T. C.; Truhlar, D. G. Optimization of vibrational coordinates, with an application to the water molecule. J. Chem. Phys. 1982, 77, 3031−3035. (18) Carter, S.; Culik, S. J.; Bowman, J. M. Vibrational selfconsistent field method for many-mode systems: A new approach and application to the vibrations of CO adsorbed on Cu(100). J. Chem. Phys. 1997, 107, 10458−10469. (19) Bowman, J. M.; Carrington, T., Jr.; Meyer, H. D. Variational quantum approaches for computing vibrational energies of polyatomic molecules. Mol. Phys. 2008, 106, 2145−2182. (20) Meier, P.; Oschetzki, D.; Berger, R.; Rauhut, G. Transformation of potential energy surfaces for estimating isotopic shifts in anharmonic vibrational frequency calculations. J. Chem. Phys. 2014, 140, 184111. (21) Saebo, S.; Pulay, P. Local Treatment of Electron Correlation. Annu. Rev. Phys. Chem. 1993, 44, 213−236. (22) Pipek, J.; Mezey, P. A Fast Intrinsic Localization Procedure Applicable for Ab Initio and Semiempirical Linear Combination of Orbital Wave-Functions. J. Chem. Phys. 1989, 90, 4916−4926. (23) Watson, J. K. G. Simplification of molecular vibration-rotation Hamiltonian. Mol. Phys. 1968, 15, 479−490. (24) Ziegler, B.; Rauhut, G. Rigorous Use of Symmetry within the Construction of Multidimensional Potential Energy Surfaces. J. Chem. Phys. 2018, 149, 164110. (25) Pflüger, K.; Paulus, M.; Jagiella, S.; Burkert, T.; Rauhut, G. Multi-Level Vibrational SCF Calculations and FTIR Measurements on Furazan. Theor. Chem. Acc. 2005, 114, 327−332. (26) Yagi, K.; Hirata, S.; Hirao, K. Multiresolution potential energy surfaces for vibrational state calculations. Theor. Chem. Acc. 2007, 118, 681−691. (27) Cheng, X.; Steele, R. P. Efficient anharmonic vibrational spectroscopy for large molecules using local-mode coordinates. J. Chem. Phys. 2014, 141, 104105. (28) Sæbø, S.; Pulay, P. Local treatment of electron correlation. Annu. Rev. Phys. Chem. 1993, 44, 213−236. (29) Rauhut, G. Efficient calculation of potential energy surfaces for the generation of vibrational wave functions. J. Chem. Phys. 2004, 121, 9313−9322. (30) Benoit, D. M. Fast vibrational self-consistent field calculations through a reduced mode-mode coupling scheme. J. Chem. Phys. 2004, 120, 562−573. (31) Adler, T. B.; Knizia, G.; Werner, H.-J. A simple and efficient CCSD(T)-F12 approximation. J. Chem. Phys. 2007, 127, 221106. (32) Kats, D.; Manby, F. R. The distinguishable cluster approximation. J. Chem. Phys. 2013, 139, 021102. (33) Kats, D.; Kreplin, D.; Werner, H.-J.; Manby, F. R. Accurate thermochemistry from explicitly correlated distinguishable cluster approximation. J. Chem. Phys. 2015, 142, 064111. (34) Kesharwani, M. K.; Sylvetsky, N.; Martin, J. M. L. Surprising Performance for Vibrational Frequencies of the Distinguishable Clusters with Singles and Doubles (DCSD) and MP2.5 Approximations. AIP Conf. Proc. 2017, 1906, 030005.

(35) Ziegler, B.; Rauhut, G. Efficient generation of sum-of-products representations of high-dimensional potential energy surfaces based on multimode expansions. J. Chem. Phys. 2016, 144, 114114. (36) Neff, M.; Hrenar, T.; Oschetzki, D.; Rauhut, G. Convergence of vibrational angular momentum terms within the Watson Hamiltonian. J. Chem. Phys. 2011, 134, 064105. (37) Neff, M.; Rauhut, G. Toward large scale vibrational configuration interaction calculations. J. Chem. Phys. 2009, 131, 124129. (38) Neff, M.; Rauhut, G. Erratum: Toward large scale vibrational configuration interaction calculations [J. Chem. Phys. 131, 124129 (2009)]. J. Chem. Phys. 2009, 131, 229901. (39) Scribano, Y.; Benoit, D. M. Iterative active-space selection for vibrational configuration interaction calculations using a reducedcoupling VSCF basis. Chem. Phys. Lett. 2008, 458, 384−387. (40) Carbonniere, P.; Pouchan, C. Modelization of vibrational spectra beyond the harmonic approximation from an iterative variation-perturbation scheme: the four conformers of the glycolaldehyde. Theor. Chem. Acc. 2012, 131, 1183. (41) Petrenko, T.; Rauhut, G. A new efficient method for the calculation of interior eigenpairs and its application to vibrational structure problems. J. Chem. Phys. 2017, 146, 124101. (42) Hepp, M.; Herman, M. Vibration-rotation bands in ethane. Mol. Phys. 2000, 98, 57−61.

4196

DOI: 10.1021/acs.jctc.9b00381 J. Chem. Theory Comput. 2019, 15, 4187−4196