A simple state-average procedure determining optimal coordinates for anharmonic vibrational calculations

Chemical Physics Letters 610–611 (2014) 288–297

Contents lists available at ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

A simple state-average procedure determining optimal coordinates for anharmonic vibrational calculations Bo Thomsen a , Kiyoshi Yagi b , Ove Christiansen a,∗ a b

Department of Chemistry, University of Aarhus, Langelandsgade 140, DK-8000 Aarhus C, Denmark Theoretical Molecular Science Laboratory and iTHES, RIKEN, Hirosawa 2-1, Saitama 351-0198, Japan

a r t i c l e

i n f o

Article history: Received 28 April 2014 In final form 16 July 2014 Available online 25 July 2014

a b s t r a c t A simple methodology for calculating state average energies in the context of vibrational self-consistent field (VSCF) is suggested. The suggested state average energy is employed in the optimization of coordinates for anharmonic vibrational wave function calculations where an orthogonal matrix of transformation between normal and optimized coordinates is variationally optimized. The convergence to the exact limit for approximate vibrational configuration interaction and vibrational coupled cluster wave functions is studied, comparing the performance of state average optimized coordinates, ground state optimized coordinates and standard normal coordinates. Exploratory calculations are presented for water, formaldehyde, the water dimer and trimer, and ethylene. © 2014 Published by Elsevier B.V.

1. Introduction The choice of coordinates is an important factor in the description of the internal dynamics of molecules. For accurate computation of anharmonic vibrational frequencies the standard choices includes rectilinear normal coordinates and curvilinear internal coordinates. Although well tested and successful for many purposes these two general choices are not entirely perfect in all circumstances. The normal coordinates are calculated based on a second order approximation to the potential, and as such they become less relevant when the amplitude of the motion is large. The couplings between some vibrational motions will in accurate anharmonic calculations be found to be larger in the rectilinear normal coordinates than in curvilinear internal coordinates. On the other hand, the use of internal coordinates is not straightforward requiring generally more effort for selecting the set of coordinates for each particular system under study as well as for determining the rather complex kinetic energy operator. Early works [1–5] have suggested to instead define the coordinates as being optimal for the vibrational energy calculated with the real anharmonic potential. Although this direction did not receive much immediate follow up work, there has recently been a revival of interest in this general idea. One of the authors [6] has implemented and

∗ Corresponding author. E-mail addresses: [email protected] (B. Thomsen), [email protected] (K. Yagi), [email protected] (O. Christiansen). http://dx.doi.org/10.1016/j.cplett.2014.07.043 0009-2614/© 2014 Published by Elsevier B.V.

performed efficient numerical optimization of a set of coordinates using the vibrational self consistent field (VSCF) method for providing the ground state energy in an anharmonic potential and determining the optimized coordinates as an orthogonal transformation of the normal coordinates of the system. We denote this approach optimized coordinate VSCF, oc-VSCF. There are significant advantages of the use of these rectilinear optimized coordinates. In particular, it has been illustrated that the convergence of approximate anharmonic wave function to the exact limit is improved using the optimized coordinates relative to normal coordinates. The ground state optimized coordinates have been used in subsequent calculation of many states using different methods, including vibrational configuration interaction (VCI) [6], vibrational coupled cluster (VCC) [7] and vibrational quasi-degenerate perturbation theory (VQDPT) [8], with very encouraging results. In the mentioned recent studies the coordinates were made optimal for the VSCF ground state energy only. In the majority of applications of anharmonic vibrational theory, not only one vibrational state is of interest, but rather several states and transitions between them are in focus. It is an obvious question, whether optimization of the coordinates for other states than the ground state would give substantially different results. The second order Taylor expansion of the potential employed to determine normal coordinates, gives a realistic description for infinitesimal distortions from the equilibrium structure. Optimizing the coordinates relative to vibrational motion in the vibrational ground state includes the effects of a significant part of the potential energy surface. However, excited vibrational wave functions clearly

B. Thomsen et al. / Chemical Physics Letters 610–611 (2014) 288–297

involves larger distortions, and the optimal set of coordinates will depend on the state under study. If all states are independently optimized the whole procedure becomes quite elaborate. Furthermore it raises some theoretical concerns. E.g. the ground state and the excited state will have different sets of coordinates complicating calculation of vibrational transition probabilities. In this work we follow a chain of arguments leading to a very simple method for optimizing coordinates including more states based on a stateaveraging procedure. In this method, a set of virtual states of a standard VSCF calculation is used to define a state-average energy obtainable with a computational cost that is essentially the same as an ordinary VSCF calculation. Employing this method in coordinate optimization we obtain an approach denoted optimized coordinate state-average virtual VSCF (oc-sa-vVSCF) that is applicable in practice to any system in which oc-VSCF is applicable. The performance of the new set of coordinates is explored in calculations on water, formaldehyde, the water dimer and the water trimer, as well as ethylene. The letter is organized as follows. In the next section we describe the theory and the chain of arguments leading to our new method. In Section 3 we describe the computational details, before we in Section 4 present the results. Section 5 presents a summary. 2. Theory

exact results, and to obtain wave functions that in more concise manners represent the important physical effects. In the coming section we shall describe briefly the VSCF ground state energy as choice of E, before describing generalizations aiming at including more than one state in the optimization. It is important that the Hamiltonian is in a form that is invariant under the transformation. Obviously the exact Hamiltonian is invariant, but approximate potentials are not necessarily invariant. One example that is invariant is a Taylor expanded force field such as a quartic force field, motivating this choice in the particular examples. An example that is not invariant is a many-mode expansion of the potential truncated to an order less than total number of modes. 2.2. Vibrational Self-consistent field and beyond A cornerstone in this work is the use of VSCF. The VSCF energy ˆ i Ei = i |H|

(3)

for a given set of coordinates and a given vibrational Hamiltonian, ˆ is optimized under variation of the orthonormal one-mode funcH, tions imm (Qm ) entering into the VSCF direct product wave function ansatz i (Q) =

2.1. Coordinate optimization

289

M 

smm (Qm ) = i11 (Q1 )i22 (Q2 ). . .iMM (QM ).

(4)

m=1

In coordinate optimization of orthogonal coordinates we relate two sets of orthogonal coordinates to each other by an orthogonal transformation of the coordinates Q˜ = U T Q

(1)

As input we use a set of M = 3Natom − 6 (M = 3Natom − 5 for a linear molecule) normal coordinates as Q = (Q1 , Q2 , . . ., QM )T , which in turn is given in terms of the orthogonal transformation of mass-weighted cartesian displacement coordinates corresponding to diagonalizing the mass-scaled molecular Hessian. We then seek to identify the rotation matrix U such that a given success quantity denoted E is minimized (or maximized). The concrete procedure used in previous works has been an iterative Jacobi-sweep algorithm. Here, inside each iteration pairs of coordinates are rotated in a series of Givens rotations. Thus, the total rotation matrix becomes

 n=1



M(M−1)/2

M(M−1)/2

U = U(1) U(2) · · ·U(niter ) =

(1)

Un

n=1



M(M−1)/2 (1)

Un · · ·

n=1

iter ) U(n n

The i vector is an M-dimensional index vector defining which onemode function for each mode is ‘occupied’ in the VSCF state. In the process of going from having the potential available at a restricted set of grid points, to a representation in a vibrational structure program, it is computationally convenient to work with a Hamiltonian in the sum over products (SOP) form, H=

Nt   t=1

(5)

m∈mt

Here Nt is the number of terms in the operator, mt is a mode combination(MC) denoting the set of active modes for the term t. The active modes for a term are the set of modes with one-mode opert ators hm,o different from the unit operator in that particular term. The superscripts m, ot defines which particular one-mode operator is present for a given mode m and term t. An example of a Hamiltonian operator in the form of Eq. (5) is 1 ∂ + V Taylor 2 ∂Qm2 2

H=−

(6)

m

(2)

where the products over n are over all M(M − 1)/2 pair rotations. Further details of the optimization procedure will not be described here. We refer to Ref. [6] for general aspects and the first implementation, and to Ref. [7] for the embarrassingly parallel implementation used here. The essential feature is that fast evaluation of a success quantity E is needed for efficiently performing the niter M(M − 1)/2 one-dimensional optimizations of the rotation angles in the Jacobi-sweep procedure. That is, the rotation angles are chosen to minimize E. The outcome of the procedure obviously depends on the quantity E to be optimized. So far the VSCF ground state energy has been used, making the procedure comply to the quantum mechanical variational principle. It should be noted, that the exact energy is invariant under the choice of coordinates. However, the quality of approximate energies and wave functions will depend on the choice of coordinates. Thus, the optimization is a way to improve approximate treatments by providing better convergence to the

t

hm,o

ct

where VTaylor is a Taylor expanded potential, for example a quartic force field. This is the concrete Hamiltonian used in the calculations presented in this work. The variational procedure leads to the non-linear VSCF equations F m,i imm (Qm ) = m m (Qm ) im im

(7) imm (Qm )

from here on denoted modals. for the one-mode functions The effective mean-field operator for mode m, Fˆ m,i , is essentially obtained by averaging the Hamiltonian over all other modes. In the formulation of Ref. [9] the mean field operator can be written as F m,i =



t∈{tact }

(ct



hm,t )hm,t , im im

(8)

m ∈mt \m

here only the terms that are active for the particular mode m is included in the mean field operator, thus {tact |t : mt ∩ {m} = {m}}.

290

B. Thomsen et al. / Chemical Physics Letters 610–611 (2014) 288–297

The total energy is obtained from Eq. (3) using the one-mode functions that satisfy the VSCF eigenvalue equations and becomes with the Hamiltonian in Eq. (5) E=

  ct

t

hm,t . im im

(9)

m∈mt

There are some variation of the formulation of the mean-field operator in the literature relative to exclusion or full or partial inclusion of inactive terms. This does not affect the total wave function and the total energy, and only the m energies by a constant. Accordim ingly it has no consequence for the final predictions of the theory. With exclusion of inactive terms and the mean-field operator of Eq. (8) an Mn computational scaling is achieved per VSCF ground state energy calculation for an M-mode system with a given Hamiltonian including up to n-mode couplings, see Ref.[9]. Excluding inactive terms and employing screening methods, the computational scaling is reduced towards linear for large systems. In this section we have so far described variational modal optimization of the VSCF energy. In the oc-VSCF we in addition optimize the coordinates Q˜ for a particular state. One should here note that the limitations of VSCF in this case is turned into a benefit. If an exact wave function is used, all sets of coordinates would be equally good. Optimization of the coordinates for the VSCF energy has the attractive feature, of optimizing the coordinate such that the modal-optimized direct product wave function gives minimum energy. E.g. if the system is described by a Hamiltonian that is separable into one-dimensional uncoupled anharmonic oscillators, VSCF is exact. If such a system is represented in a set of coordinates that does not provide this additive separability in the Hamiltonian, VSCF is not exact. However, with variational coordinate optimization, oc-VSCF is able to rotate the coordinate to obtain the exact minimum energy and wave-function. Thus, the fact that VSCF is not highly accurate is acceptable from the point of view of the coordinates optimization. Variational optimization of the coordinates for the VSCF ground state will lead to a lower, and thus better VSCF energy for the variational VSCF wave function. For the given Hamiltonian, the coordinate optimization thus reduces the correlation problem for the ground state. The oc-VSCF approach has been successfully used to optimize coordinates for the ground state only, and these coordinates have subsequently been used in advanced vibrational structure calculations such as VCI, VQDPT and VCC with success. Of course, there may be cases where the VSCF is simply qualitatively unreasonable already for the ground state and in these cases e.g. a vibrational Multiconfigurational self-consistent field (VMCSCF) approach[10] would be more appropriate. Nevertheless, for very many molecules VSCF is qualitatively fully correct and optimizing the ground state energy using VSCF fully satisfactory. There is, however, another valid concern. Most often we are interested in not only the ground state, but rather in several or many states and the energy differences and transition properties between them. Correspondingly, it is a relevant question whether extending the oc-VSCF idea to simultaneously treating more states in the optimization procedure would give substantially different results. 2.3. Multi-state coordinate optimizations 2.3.1. General aspects on State average energies The most obvious extension to the many-state case is to employ the procedures above to optimization of the excited states on a state by state fashion. There are, however, a few fundamental problems in that. First, the different states would be described by different coordinates. Second, we have not in any of our optimization steps ensured that if applied to an excited state, the resulting state is orthogonal to the ground state. Third, the excited state may with

much higher probability than the ground state, be qualitatively poorly described by a single configuration as employed in VSCF. Fourth, if we single out a state, we have to define ways to keep the focus on this state in the optimization. For example, if we want to optimize the coordinates for the state (1,0,0) we need to ensure that the coordinate optimization does not simply switch the coordinates so as to optimize another state, such as (0,1,0). Although perhaps a technical problem, it illustrates that it is tedious to keep track of many vibrational states in such optimizations. Indeed, the explosion in the number of available states with increasing size of the molecule as well as with increasing energy, is itself a discouragement of such state by state procedures. These things together make it very complicated to both obtain and apply the resulting coordinates using a state by state approach. We therefore proceed to consider some kind of average of the energy (or other measures) to make a compromise, such that calculations on several states are possible, with one set of coordinates. Thus, we introduce a state-average energy Eave =



wI EI

(10)

I

where EI is the energy of state I and wI is the weight in the averaging. We further would like to enforce that the success quantity, here the state average energy, should be invariant under interchange of two coordinates. The energy of each state could be obtained in separate VSCF calculations. However, the different VSCF states still have different modals and may not be orthogonal to the ground state, which may be both a practical and theoretical concern. The most logical and balanced way seems to be to achieve the result through a stateaverage VMCSCF approach, sa-VMCSCF. Thus, each state would be described by a wave function I =



CrI r

(11)

r

and energy EI = I |H|I 

(12)

Here r denotes direct product configuration similar to the previous VSCF i state, only now with several different occupations allowed in the summation in Eq. (11). Using the sa-VMCSCF energy in the coordinate optimization we obtain what we would denote the oc-sa-VMCSCF results, as a fully consistent variational treatment of both modal optimization and coordinate optimization for several states. As such this is a theoretically very attractive generalization of the oc-VSCF procedure. Multi-configurational methods in the vibrational context has received increasing attention [10–14]. However, we are not aware of any state-average VMCSCF implementation. Furthermore, the computational cost will increase as the system size increase due to active space in VMCSCF (see Refs.[10,13]) for discussions), but also because the number of states will increase. E.g. we may define a modal space and include all possible configurations from this space in the wave function of Eq. (11) and all resulting states in the average energy of Eq. (10) (an ultimate type of complete-active-space (CAS) strategy). While this is very attractive, it would be a rather costly procedure. If we wish to include all states that can be constructed from having each mode in one of its lowest Nm modals, the number of states scales exponentially with the size of the system Ntot,states =

M 

Nm = N M

(13)

m=1

where the last equality applies if all modes have the same number of modals, Nm = N. This will become a problem with increasing M.

B. Thomsen et al. / Chemical Physics Letters 610–611 (2014) 288–297

In fact, this scaling is so severe it will be a significant problem even if fast state-specific VSCF calculations were used for the energies in Eq. (10). We therefore choose presently to look for other ways ahead, ultimately arriving at an approximate approach including all these states implicitly, but without actually suffering the cost of treating them individually.

291

2.3.3. State average energies using virtual VSCF states and factorizable weights We choose the weights to be factorizable as wr =

M 

wrmm

(16)

m=1

and appropriately normalized to sum to one, meaning

 2.3.2. State average energies using virtual VSCF states To make things simple from the outset we require we only want to deal with one set of coordinates and one set of modals, common for all states. We choose to describe also the excited states by the modals optimized in the ground state VSCF calculations. These modals are the VSCF eigenfunctions obtained as the solutions of Eq. (7) differing from the ground state. We denote these states as the virtual VSCF states, as they are obtained by exciting individuals modes from the ground state modal to a virtual modal. Let each state be described by a vector of modal indices per state, r = r = (r1 , · · ·, rm2 , · · ·rM ). The energy of each state is the simple expectation value Er = r |H|r 

(14)

Assigning weights to individual states defined by the index vector r we have a general virtual VSCF state average approach. For the reasons mentioned above, it becomes prohibitive to explicitly sum over the many available states for larger systems, and we have to specialize to either a limited set of states or special types of weights (next section). In this section we shall be explicit on a simple virtual VSCF state average approach including an average of the ground state and all fundamental energies with equal weights. In this case we have

Eave =



This is crucial for the simplicity and practicality of the approach presented here. With this approach we can obtain a very simple way of calculating an average energy for many states. For example, if we for a mode m include the Nm lowest modals with equal weights we have a factor wrmm =

(18)

wr =

M

1

m=1

Nm

Eave =



wr r |H|r 

r



r |H|r 

(15)

=

ct

···

r1

r2

   ct

(

m∈mt

  ct

t

rM

ct

m∈mt

hrm rm wrmm )

  (

m/ ∈mt

hrm rm wrmm =

m∈mt rm



t

hm,o m |r 

M 

wrmm

m=1

wrmm

m/ ∈mt

wrmm )

rm

  ct

t



m∈mt

r t

m,ot

m,ot

r |

m [hm,o rm rm wrm ])

(

rm

  t

Here m are the one-mode modal energies obtained by solvrm ing the VSCF effective eigenvalue equations in Eq. (7). The ground state has an index occupation vectors with only zeroes, i = 0. The fundamental state energies are E1m for excitation to a state described by an index vector with one for mode m and zero otherwise. In the above equations we have used the following rewrite of the fundamental state energy that is possible for one-mode excited states as the fundamentals, Em = 1m |H|1m  = E0 + (1m |H|1m  − 0 |H|0 ) = − m ). The average E0 + (1m |F m,0 |1m  − 0m |F m,0 |0m ) = E0 + (m 1 0 energy in Eq. (15) is simple to calculate since after the usual VSCF energy calculation a very fast summation of the M one-mode zeroth-order excitation (m − m ) is performed. We denote the set 1 0 of coordinates obtained from this approach as oc[FU] to be contrasted with the coordinates in the following section. Overtones can be included in a similar fashion, while the inclusion of general combinations is more cumbersome. We therefore turn to a different approach.

wrmm =

   

t

=

M 

m=1

t

m

m

(19)

which is obviously just one divided by the total number of states that can be constructed based on the set of modals allowed to be occupied. The weighting can be more advanced, but in the absence of better guidelines we believe the simple strategy suggested by Eqs. (16)–(18) is sufficient. With the direct product form of the weights we can rewrite the state average energy as

=

 1  m 1 (E0 + (E0 + m − m )) = E0 + (1 − m 1 0 0) M+1 M+1

1 Ntot,states

=

r

 1 E1m ) = (E0 + M+1

m

1 Nm

for that mode. Using a similar pattern for all modes we have in this particular case

r

=

(17)

rm

=

wr r |H|r 

wrmm = 1

t h˜ m,o

(20)

m∈mt

By pre-evaluating for each mode and each one-mode operator the numbers t h˜ m,o =



t

m hm,o rm rm wrm

(21)

rm

the evaluation of the state average energy can be efficiently evaluated by summing over all terms in the Hamiltonian. This evaluation is as simple as the evaluation of the VSCF energy for a given set of modals. Thus, the extra cost is the evaluation of the MNop transformed integrals assuming Nop operators per mode, and the summation for the total energy scaling as Mn for a Hamiltonian with up to n-mode couplings. In conclusion, the evaluation of the state-averaged VSCF energy is thus obtainable at very little extra cost after a standard VSCF calculation. We also note in passing that Eave is, as is the case for the VSCF energy, a size extensive quantity. The inspiration for the simple state-average VSCF ansatz comes from the v-VSCF approach to temperature averaging [15] and we may denote it state average virtual-VSCF, sa-vVSCF. In the case of

292

B. Thomsen et al. / Chemical Physics Letters 610–611 (2014) 288–297

temperature averaging the weights are chosen in accord with Boltzman populations for individual states, but factorizable as in Eq. (16) following a few assumptions and approximations. Such temperature dependent weights could also easily have been used here with similar computational costs. However, even up to room temperature and beyond, highly excited states are essentially not thermally populated, so this is not an attractive manner for optimizing the coordinates for higher-lying modes (e.g. X H stretch fundamentals). In particular, X H frequencies was found to be a case, where the optimized coordinates were particularly interesting in comparison to standard normal coordinates. We therefore choose to proceed with the average weights above, that in the simplest case can be defined by one common integer, N, defining the number of levels taken into account in the averaging per mode. If the Nm for different modes are not identical, care is required in avoiding ambiguous results arising from rotation of coordinates in a similar manner as commented before for the state-specific VSCF calculation. We therefore generally use a common integer for all states, such that the state average energy is guaranteed to fulfill the criteria of being invariant under interchange of two coordinates. The optimized coordinates can in this case be denoted oc[N]-saVSCF. The case oc[1] correspond to only the ground state energy, and is thus identical to the previous ground state only oc-VSCF coordinates. The others corresponding to state averaging including higher and higher excited levels. E.g. oc[2] includes all fundamentals and all combinations with excitation to the lowest excited modal per mode, while oc[3] in addition includes overtones, and all combinations excited at most to the second excited state per mode, etc.

3. Computational details For all molecules we employ quartic force-fields. Quartic force fields are known to be not fully satisfactory for highly accurate calculations of vibrational transition energies. However, they are attractive from the point of view of coordinate optimization due to their invariance under coordinate transformation, and are thus well-suited for this study. We will restrict ourselves to consider only the performance of different optimized coordinates measured in terms of the convergence of vibrational wave function methods towards converged reference values. It is not the purpose of this paper to produce absolute high accuracy for comparison with experiments. Instead we will explore the performance of the methodology to several different systems including semi-rigid molecules as formaldehyde and ethylene, as well as water and the water dimer and water trimer. The quartic force field has been calculated at the MP2 level with the cc-pVTZ basis set for water and with the cc-pVDZ basis sets for others (i.e., formaldehyde, ethyelene, water dimer, and water trimer) using the SINDO program [16,17]. The electronic structure calculation has been carried out using Gaussian 09 [18]. We note that the water dimer and water trimer are considered in the model [19–21] where the modes corresponding to intermolecular degrees of freedom are frozen, and the vibrational problem is solved for the remaining coupled intra-molecular degrees of freedom. Thus, the water dimer is effectively a six-dimensional problem and the water trimer a nine-dimensional problem. We stress that there are real and significant couplings between the individual molecules of the dimer and trimer in this model. The state-average vVSCF calculation and oc[N] coordinate optimization has been implemented in MidasCpp [22]. The VSCF state average calculations increases the cost of a VSCF calculation by a neglible fraction maintaining Mn order computational scaling [9], and the optimization algorithm follows our previously described embarrassingly parallel implementation [7]. MidasCpp was also

used for the subsequent VCI and VCC calculations. All vibrational structure calculations used a basis set consisting of the lowest eight harmonic oscillator eigenfunctions for each mode in the molecule. The exponents for the harmonic oscillator functions were generated from the one-dimensional quadratic terms of quartic forcefield. For all systems we determine oc[N] coordinates where N =1,2,3,4 is used in the calculation of the state average energy. Following the determination of the coordinates and the potential in the optimized set of coordinates, we perform a series of vibrational wave function calculations, at the level of VCI[3], VCC[2], VCC[3], and VCC[3pt4]. A modal basis set of eight modals per mode were used for water, formaldehyde and water dimer, while six modals per mode were used for water trimer and four for ethylene. All excitation energies are calculated using vibrational response theory. The results are compared to suitable references calculated in this work. For water, formaldehyde and the water dimer we use oc-FVCI, while for water trimer and ethylene we use oc-VCC[5].

4. Results 4.1. Coordinates We shall first discuss qualitatively the changes in coordinates from normal coordinates to the optimized coordinates using either the VSCF ground state energy, or a state-average VSCF energy. For this purpose we shall consider the transformation matrices between normal coordinates and the various optimized coordinates. These transformation matrices are given in Fig. 1, where a diagonal matrix indicates that no rotations were possible. In general we observe more or less block-diagonal transformation matrices. Thus we observe substantial changes in the coordinates for some, but not all, degrees of freedom already when moving from normal coordinates to the oc[1] level. This mirrors the changes seen in previous works [6–8] namely that the X H stretches are more localized than their symmetric normal coordinate counterparts. Furthermore the bends of the water trimer are transformed from the coordinates delocalized over all three molecules to those localized on individual molecules when the rotations of the oc[1] model are applied. The oc[FU], oc[2] and the higher oc[N] models change the coordinates further, as seen in Fig. 1. In particular, the intermediateenergy vibrations start mixing to a larger extent as we move up the hierarchy. The notable exception is the water monomer which remains unchanged as we add the fundamentals or increase N, suggesting that the localization of the O H stretches is optimal for the state average energies as well. Formaldehyde undergoes some changes in oc[FU] or oc[N] with increasing N, displaying a stronger mixing of the C O stretching and H C H bending motions, though this is barely noticeable in a visual examination of the coordinates. In the case of ethylene, the C H stretching modes are localized around a single bond using the oc[1] and subsequent methods, as seen for the O H stretches of water and water clusters (see below). To provide a more quantitative measure of the changes we turn to Fig. 2, which gives the norms of the coordinate vector elements for each atom. It is clearly seen that the C H stretching vectors (modes 9–12) of all optimized coordinates have a dominant component from a single atom (= hydrogen). Furthermore, adding more states in the state average measure is found to blend the out of plane vibrations, as seen in Fig. 1. The inclusion of fundamentals in oc[FU] localizes the motion of mode two and three to three atoms, and adding the combinations in oc[2] localizes the motions of mode two, three and four. A visual representation of the vibrational modes is depicted in Fig. 3 for modes two and nine. The rest of the coordinates are given in the supplementary information in Figs. S1 and S2.

B. Thomsen et al. / Chemical Physics Letters 610–611 (2014) 288–297

293

Fig. 1. Graphical illustration of the transformation matrices between normal and optimized coordinates for (a) water, (b) formaldehyde, (c) ethylene, (d) water dimer and (e) water trimer in the oc[1](leftmost), oc[FU](left middle), oc[2](right middle) and oc[3](rightmost) models. The matrices are ordered by the harmonic frequencies of the vibrational modes increasing from left to right and top to bottom. The absolute value of each element is represented by a color going from white(0) to black(1).

We have so far seen that the inclusion of the higher energy states leads to vibrational modes with more local nature, which may be what is generally expected. Nonetheless, our oc[N] models may well yield non-local coordinates, since the localization itself is not the measure of optimization in our method. In fact, we do see an opposite effect in the case of water clusters, where the bending modes change from being localized to each water in the oc[1] and oc[FU] models to being delocalized in the oc[2] and higher models. The coordinates are illustrated in Fig. 4 for two selected modes of water dimer. The figure shows that the bends become delocalized in the oc[2], though the O H stretches are kept localized to each O H bonds. The reason for the delocalization of the water cluster bending motions in oc[N] may be traced to the combination states included in the state-average energy. In the oc[2] model, for example, the singly excited combinations are included in the average in addition to the ground state and the fundamentals. Here, note that the combination state has an opposite character to that of the underlying coordinates, i.e., the combination in terms of local bending coordinates corresponds to a simultaneous excitation of the two waters,

and vice versa. Therefore, introducing the combinations into the averaging adds energy contributions that, since they originate from states with delocalized motions, may increase delocalization of the coordinates in the variational optimization, and this is indeed the effect found by comparing oc[FU] and oc[2] in the water dimer. The potential delocalization of coordinates caused by the combinations in oc[N] may or may not benefit the subsequent correlation calculations, as we examine in the next section.

4.2. Energies In Table 1 we have collected the average deviations of the calculated vibrational excitation energies relative to the reference FVCI results. As for the wave function model, the errors decrease in the sequence of VCI[3] > VCC[2] > VCC[3] > VCC[3pt4] as expected. Repeating the principal findings of previous work for formaldehyde and the water clusters [7], the results based on oc[1] are found with smaller errors relative to those based on normal coordinates. The same trend is very pronounced for ethylene. Thus, it is confirmed

(c) The distribution of motion in optimized coordinates for the fundamentals only model.

Mode 12

Mode 11

Mode 10

Mode 9

Mode 8

Mode 7

Mode 6

Mode 5

Mode 3

Mode 4

Mode 12

Mode 11

Mode 10

Mode 9

Mode 8

Mode 7

Mode 6

Mode 5

Mode 4

Mode 12

0

Mode 11

0

Mode 10

0.2

Mode 9

0.2

Mode 8

0.4

Mode 7

0.4

Mode 6

0.6

Mode 5

0.6

Mode 4

0.8

Mode 3

0.8

Mode 2

1

Mode 1

1

Mode 3

(b) The distribution of motion in optimized coordinates for N = 1.

Mode 2

(a) The distribution of motions of the atoms in normal coordinates.

Mode 2

Mode 12

0 Mode 11

0 Mode 10

0.2

Mode 9

0.2

Mode 8

0.4

Mode 7

0.4

Mode 6

0.6

Mode 5

0.6

Mode 4

0.8

Mode 3

0.8

Mode 2

1

Mode 1

1

Mode 1

B. Thomsen et al. / Chemical Physics Letters 610–611 (2014) 288–297

Mode 1

294

(d) The distribution of motion in optimized coordinates for N = 2.

Fig. 2. The coordinate localization on single atoms in ethylene for normal coordinates(a), oc[1](b), oc[FU](c) and oc[2](d). Each color corresponds to a different atom. A weight of one would indicate that the motion is completely localized on the given atom. The modes are ordered energetically from left to right. (The reader is referred to the web version of this article for a full color version of the figure.)

 i ) and oc[2](Q i ) coordinates for the second and ninth vibration of ethylene. Fig. 3. Normal(Qi ), oc[1](Q˜ i ), oc[FU](Q

B. Thomsen et al. / Chemical Physics Letters 610–611 (2014) 288–297

295

 i ) and oc[2](Q i ) coordinates for the second and fifth vibration of water dimer. Fig. 4. Normal(Qi ), oc[1](Q˜ i ), oc[FU](Q

that there is an overall faster convergence toward FVCI for wave functions in terms of oc[1] compared to normal coordinates. A further reduction of error is observed for ethylene using the coordinates derived from the state-average energy (except for VCC[3PT4]; see below). Taking VCC[2] as an example, we see that the errors are reduced by 14% and 20% in oc[FU] and oc[2], respectively, relative to oc[1]. The improvement is found to be remarkable especially for the out-of-plane modes with the reduction of error being as large as 28% and 44%. It is also notable that the C H stretching modes are improved as well, even though the coordinates remain the same as oc[1], i.e., local C H stretches. This result clearly exemplifies the advantage of optimized coordinates based on the state-average energy, here seen clearly for the localized out-of-plane modes. Proceeding to oc[3] and oc[4] exhibits little improvement or sometimes an increase in the error, for example, for the C H stretches. This may be due to a mismatch between the

states included in the average and the ones we are calculating (i.e., fundamentals), a match of the states in the FVCI and in oc[2], and so on. Nevertheless, they still remain more accurate than oc[1]. While oc[N] (N ≥ 2) is an improvement for ethylene, the opposite is found to be true for the water clusters. For both clusters holds that the average errors for the fundamentals increase substantially by an order of magnitude proceeding from oc[1] to oc[2] in VCC[2] and VCC[3], though the absolute value itself is smaller than that for formaldehyde and ethylene. We have also found a similar trend for the overtones and combinations of water dimer, which indicates that the mismatch between the states to be calculated and those included in the average is not the only cause. Rather, the error is interpreted to be related to that these coordinates give larger mode coupling terms. Indeed, inspecting the average force constants of QFF, given in Table 2, supports this picture. Recall that oc[N] (N ≥ 2) incorporates combinations of the bending modes in

Table 1 Average deviations (in cm−1 ) relative to converged reference (see Section 3) values for vibrational excitation energies, of the fundamentals unless otherwise noted, for different set of coordinates. Normal coordinates are denoted nc in the table, and the optimized coordinates are label by the value of N, going from 1 up to 4. VCI[3] nc a

Water Formaldehyde Ethylene Ethylene(C Hc ) Ethylene(OoPd ) Water-dimer Water-dimer(overtonese ) Water-dimer(combf ) Water-trimer

VCC[2] oc[1]

b

0.00 1.69 31.15 27.03 31.46 0.20 0.67 11.36 3.52

oc[FU] b

0.00 0.91 22.93 24.91 19.84 0.12 0.46 1.46 0.25

b

0.00 0.90 18.59 20.67 16.77 0.14 0.46 1.45 0.26

oc[2]

oc[3] b

0.00 0.91 16.69 19.19 16.23 0.58 2.19 2.61 1.24

oc[4] b

0.00 0.89 15.92 18.07 15.68 0.60 2.27 2.73 1.26

b

0.00 0.90 15.99 18.34 15.68 0.59 2.29 2.78 1.26

VCC[3]

Formaldehyde Ethylene Ethylene(C Hc ) Ethylene(OoPd ) Water-dimer Water-dimer(overtonese ) Water-dimer(combf ) Water-trimer a b c d e f

nc

oc[1]

oc[FU]

oc[2]

oc[3]

oc[4]

18.53 2.17 13.78 15.86 15.75 0.51 4.16 72.00 2.51

3.30 1.50 10.27 8.42 13.81 0.21 4.67 7.80 0.21

3.29 1.47 8.86 7.27 10.08 0.21 4.70 7.78 0.21

3.35 1.42 8.22 7.09 7.74 1.03 9.98 7.65 1.19

3.46 1.40 8.13 7.28 7.62 1.04 10.31 7.84 1.21

3.52 1.40 8.22 7.63 7.59 1.04 10.41 7.96 1.20

VCC[3PT4]

nc

oc[1]

oc[FU]

oc[2]

oc[3]

oc[4]

nc

oc[1]

oc[FU]

oc[2]

oc[3]

oc[4]

0.40 3.48 3.77 3.65 0.01 0.07 11.40 0.41

0.08 2.35 3.27 2.05 0.00 0.01 0.40 0.00

0.08 2.08 3.02 1.45 0.00 0.01 0.40 0.00

0.08 2.13 3.00 1.85 0.03 1.11 1.51 0.05

0.09 2.01 2.90 1.83 0.03 1.14 1.53 0.05

0.09 2.02 2.91 1.84 0.03 1.13 1.53 0.05

0.17 1.50 2.78 0.88 0.00 0.01 1.10 0.14

0.03 0.75 0.99 0.74 0.00 0.00 0.07 0.00

0.03 0.77 0.87 0.75 0.00 0.00 0.07 0.00

0.03 0.84 0.88 1.06 0.01 0.38 0.36 0.02

0.03 0.85 0.85 1.07 0.01 0.39 0.36 0.02

0.03 0.85 0.83 1.08 0.01 0.39 0.36 0.02

The fundamentals, first overtones, and combination tones. VCI[3] corresponds to FVCI for water, i.e. the excitation energies are independent of the choice of coordinates. C H stretching modes. Out-of-plane modes. The first overtones (double excitation in one mode). The first combination tones (single excitions in two modes).

296

B. Thomsen et al. / Chemical Physics Letters 610–611 (2014) 288–297

Table 2 The mean absolute size of the force constants of the quartic force field for the water dimer in normal coordinates (nc), oc[1], oc[FU] and oc[2] given in Hartree. Mode coupling

Order

1 1 1 2 2 2 3 3 4

2 3 4 2 3 4 3 4 4

Mode coupling

Order

1 1 1 2 2 2 3 3 4

2 3 4 2 3 4 3 4 4

a b

All modesa nc

oc[1]

oc[FU]

oc[2]

1.14 × 10−4 1.69 × 10−6 3.58 × 10−8 4.58 × 10−11 1.07 × 10−6 2.27 × 10−8 8.38 × 10−8 8.90 × 10−9 5.06 × 10−8

1.14 × 10−4 3.11 × 10−6 6.06 × 10−8 1.45 × 10−6 1.62 × 10−7 7.91 × 10−9 7.55 × 10−8 9.93 × 10−10 9.35 × 10−8

1.14 × 10−4 3.11 × 10−6 6.06 × 10−8 1.49 × 10−6 1.54 × 10−7 7.77 × 10−9 7.54 × 10−8 9.98 × 10−10 9.32 × 10−8

1.14 × 10−4 3.09 × 10−6 6.06 × 10−8 2.21 × 10−6 1.67 × 10−7 7.97 × 10−9 1.93 × 10−7 5.54 × 10−9 2.41 × 10−7

Bending modesb nc

oc[1]

oc[FU]

oc[2]

2.97 × 10−5 1.28 × 10−7 6.17 × 10−10 1.84 × 10−11 2.26 × 10−7 1.02 × 10−8 5.08 × 10−9 6.94 × 10−10 7.43 × 10−10

2.97 × 10−5 1.27 × 10−7 6.08 × 10−10 9.95 × 10−8 2.33 × 10−7 1.22 × 10−8 1.21 × 10−9 6.46 × 10−11 6.29 × 10−11

2.97 × 10−5 1.27 × 10−7 6.01 × 10−10 1.51 × 10−7 2.32 × 10−7 1.21 × 10−8 1.21 × 10−9 6.52 × 10−11 6.34 × 10−11

2.98 × 10−5 4.99 × 10−8 2.69 × 10−10 1.36 × 10−6 2.61 × 10−7 1.27 × 10−8 6.88 × 10−9 1.22 × 10−9 1.30 × 10−9

The average over all combinations of modes. The average over the combinations of modes that involve one or more bending modes (Q1 and Q2 ) in the indices.

the state average energy to be optimized, and as a result there is an increased mixing of bending modes of separate water molecules. Inspection of the force constants in the different coordinates shows that the one- and two-mode terms are similar or slightly increasing in size from normal coordinates (nc) to oc[1] and oc[FU], while the three- and four-mode terms are decreasing in magnitude. However, proceeding to oc[2] shows that the three- and four-mode terms increase in magnitude again to be more similar to the level obtained in normal coordinates. This originates in particular from the coupling terms involving the bending coordinates, in which the strength is 1–2 orders of magnitude larger. This means that the coupling patterns are clearly more complicated in the oc[2] coordinates than in oc[1] and oc[FU]. While optimizing an average energy may optimize the coordinates in some sense, the delocalization of the bending modes in oc[2] spreads out the couplings such as to require longer wave function expansions for both VCI and VCC. The present result shows that VCC favors locality coded in the coordinates, which should be further linked to the size-extensivity of the VCC wave function; see Ref. [7] for more details. We emphasize that in contrast to oc[2], oc[FU] is found to be hardly deteriorated, which is reasonable since the locality of the coordinates is in line with that of optimizing the fundamental excitations. Thus, the nature of the states included in the state average optimization is important. Finally, let us turn to the VCC[3PT4] results also presented in Table 1. It is seen that oc[1] achieves sub 1 cm−1 accuracy in this model, oc[FU] keeps the same quality, but oc[N] (N ≥ 2) show an increase in the error. As the level of theory is closer to FVCI and FVCI is independent of the coordinates, the numerical effects are smaller compared to VCC[2] and VCC[3], but the trends are still in line with the above discussion. As we include more states in the average, the performance of VCC[3pt4] does not behave uniformly better, but it is encouraging to see that oc[FU] gives overall a good performance comparable to oc[1].

how the method can be implemented with only very modest increase in computational cost, compared to standard VSCF calculations. In exploratory calculations for water, formaldehyde, the water dimer, the water trimer, and ethylene, it was confirmed that optimized coordinates based on variational optimization of state-average VSCF energies may lead to faster convergence of the vibrational wave function to the exact limit compared to corresponding normal coordinate expansions. The coordinates obtained from optimization of the state-average calculations including many other states than the vibrational ground state in some cases differ significantly from the coordinates optimized for the ground state only. In some cases (ethylene out of plane vibrations) these differences lead to improved convergence of the vibrational excitation energies, while in other cases (water-dimer and water-trimer) the opposite is found for some choices of defining state average energies. The choice of states included clearly affects the performance of approximate wave function calculation. The oc[FU] approach provided reasonable good results for all the fundamentals of the systems considered, in the sense of supporting higher accuracy VCC calculations, and maintaining locality. Acknowledgements K.Y. acknowleges financial support from Incentive Research Projects, RIKEN. O.C. acknowledges support from the Lundbeck Foundation, the Danish National Research Foundation, the Danish Council for Independent Research and the Danish e-infrastructure cooperation (DeiC). Appendix A. Supplementary Data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett. 2014.07.043.

5. Summary A simple methodology for calculating state average energies in the context of VSCF has been suggested and applied in the context of determining optimized coordinates. It has been shown

References [1] T.C. Thompson, D.G. Truhlar, J. Chem. Phys. 77 (1982) 3031. [2] R. Lefebvre, Int. J. Quant. Chem. 23 (1983) 543.

B. Thomsen et al. / Chemical Physics Letters 610–611 (2014) 288–297 [3] N. Moiseyev, Chem. Phys. Lett. 98 (1983) 233. ´ R.B. Gerber, M.A. Ratner, J. Phys. Chem. 90 (1986) 3606. [4] Z. Bacic, [5] A. Hidalgo, J. Zúniga, J.M. Francés, A. Bastida, A. Requena, Int. J. Quant. Chem. 40 (1991) 685. [6] K. Yagi, M. Kec¸eli, S. Hirata, J. Chem. Phys. 137 (2012) 204118. [7] B. Thomsen, K. Yagi, O. Christiansen, J. Chem. Phys. 140 (2014) 154102. [8] K. Yagi, H. Otaki, J. Chem. Phys. 140 (2014) 084113. [9] M.B. Hansen, M. Sparta, P. Seidler, O. Christiansen, D. Toffoli, J. Chem. Theor. Comp. 6 (2010) 235. [10] S. Heislbetz, G. Rauhut, J. Chem. Phys. 132 (2010) 124102. [11] F. Culot, J. Livin, Theor. Chim. Acta 89 (1994) 227. [12] H. Meyer, F. Le Quere, C. Leonard, F. Gatti, Chem. Phys. 329 (2006) 179. [13] S. Heislbetz, F. Pfeiffer, G. Rauhut, J. Chem. Phys. 134 (2011) 204108.

297

[14] W. Mizukami, D.P. Tew, J. Chem. Phys. 139 (2013) 194108. [15] M.B. Hansen, O. Christiansen, D. Toffoli, J. Chem. Phys. 128 (2008) 174106. [16] K. Yagi, K. Hirao, T. Taketsugu, M.W. Schmidt, M.S. Gordon, J. Chem. Phys. 121 (2004) 1383. [17] K. Yagi, SINDO 3.4, RIKEN, 2013. [18] M.J. Frisch, et al., Gaussian 09 Revision D.01, Gaussian Inc., Wallingford, CT, 2009. [19] G.R. Low, H.G. Kjaergaard, J. Chem. Phys. 110 (1999) 9104. [20] D.P. Schofield, H.G. Kjaergaard, Phys. Chem. Chem. Phys. 5 (2003) 3100. [21] T. Salmi, V. Hänninen, A.L. Garden, H.G. Kjaergaard, J. Tennyson, L. Halonen, J. Phys. Chem. A 112 (2008) 6305. [22] O. Christiansen, et al., MidasCpp (Molecular Interactions, Dynamics and Simulation Chemistry Program Package in C++), University of Aarhus, 2014.