THE JOURNAL OF CHEMICAL PHYSICS 142, 144115 (2015)

Automatic determination of important mode–mode correlations in many-mode vibrational wave functions Carolin Königa) and Ove Christiansenb) Department of Chemistry, Aarhus University, DK-8000 Aarhus C, Denmark

(Received 14 January 2015; accepted 13 March 2015; published online 14 April 2015) We introduce new automatic procedures for parameterizing vibrational coupled cluster (VCC) and vibrational configuration interaction wave functions. Importance measures for individual mode combinations in the wave function are derived based on upper bounds to Hamiltonian matrix elements and/or the size of perturbative corrections derived in the framework of VCC. With a threshold, this enables an automatic, system-adapted way of choosing which mode–mode correlations are explicitly parameterized in the many-mode wave function. The effect of different importance measures and thresholds is investigated for zero-point energies and infrared spectra for formaldehyde and furan. Furthermore, the direct link between important mode–mode correlations and coordinates is illustrated employing water clusters as examples: Using optimized coordinates, a larger number of mode combinations can be neglected in the correlated many-mode vibrational wave function than with normal coordinates for the same accuracy. Moreover, the fraction of important mode–mode correlations compared to the total number of correlations decreases with system size. This underlines the potential gain in efficiency when using optimized coordinates in combination with a flexible scheme for choosing the mode–mode correlations included in the parameterization of the correlated many-mode vibrational wave function. All in all, it is found that the introduced schemes for parameterizing correlated manymode vibrational wave functions lead to at least as systematic and accurate calculations as those using more standard and straightforward excitation level definitions. This new way of defining approximate calculations offers potential for future calculations on larger systems. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4916518]

I. INTRODUCTION

The interpretation of vibrational spectra of sizable (bio-)systems typically requires theoretical simulation.1,2 The current computational tools that are feasible for the calculation of vibrational spectra of sizable systems are often not accurate enough to reliably predict such spectra. High-accuracy methods based on a wave function description of the nuclear motion have traditionally only been realistically applicable to small to medium-size molecules. This is due to the steep computational scaling of such methods with system size. This steep scaling is obtained for both the (i) generation of the potential energy surface needed in order to construct the vibrational Hamiltonian and (ii) the correlated vibrational wave function itself. Our focus here is on the second problem. More precisely, we aim at more efficient and flexible formulations for constructing many-mode correlated vibrational wave functions. The basic vibrational self-consistent field (VSCF) method3–5 has been implemented in a form that allows large systems to be treated.6 For unambiguous interpretation of vibrational spectra, however, it will often be necessary to go beyond the VSCF mean field description. Using VSCF as a reference, we will formulate different estimates for the importance of different mode–mode correlations in the wave function calculation for a given Hamiltonian. Please note that a)Electronic mail: [email protected] b)Electronic mail: [email protected]

0021-9606/2015/142(14)/144115/19/$30.00

throughout this manuscript, we consider the potential energy surface to be given and do not prune the mode–mode couplings within the Hamilton operator. The estimates are then used to automatically set up parameterizations for correlated manymode vibrational wave functions of exponential vibrational coupled cluster (VCC)7,8 and vibrational configuration interaction (VCI)7,9–12-type. For truncated VCI and VCC, a loworder polynomial computational scaling with the number of degrees of freedom can be achieved. This has been illustrated by automatic derivation and analysis.13,14 Still, the polynomial scaling quickly leads to costly computations. For VCC, it has been shown that for large model systems that are actually decoupled, the computational cost can be reduced significantly by exploiting this decoupling in the parameterization.15 However, it has not been defined how to automatically exploit this in real interacting systems where there surely is a continuous range of interactions from weak to strong. The present study provides one answer to addressing this problem in the future. The correlated vibrational wave function can be expressed in terms of an expansion in Hartree products constructed as products of one-mode wave functions, the so-called modals. The latter are typically obtained in VSCF calculations. One of these modals per mode is occupied in the reference (typically VSCF) state for each of the distinguishable degrees of freedom. If we allow for every possible occupation of every mode, a large number of Hartree products can be generated, corresponding to an exponential increase in the available state space with system size. Expansion in all possible Hartree products

142, 144115-1

© 2015 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 15 Apr 2015 09:40:25

144115-2

C. König and O. Christiansen

in a given basis leads to the full vibrational configurational interaction (FVCI) ansatz. It is exact (in the given basis) but typically only applicable to very small systems. One, hence, needs to truncate the expansions for true applications. For this purpose, we structure the vibrational state space. First of all, we group all Hartree products in which the same modes are excited, i.e., not in their reference state. For instance, a configuration in which all modes except those labelled a, b, c are in reference state is assigned to the mode combination (a, b, c). Each mode combination m then gives rise to (N − 1)dim(m) Hartree products, where N is the number of modals per mode (which is for simplicity here assumed to be chosen equally for all modes) and dim(m) is the number of modes in m. Subsequently, we can write both, the VCI and the VCC ansatz, in terms of a mode combination range (MCR) defining the manifold of considered mode combinations in the wave function parameterization. A typical choice of MCR considers all mode combinations up to a certain level, i.e., number of simultaneously excited modes. Here, we aim at softening the borders between different levels and instead choose the mode combination range of a correlated many-mode vibrational wave function based on importance estimates. Clearly, the wave function parameters for a given mode combination are strongly related to the description of the physical coupling between the modes of the mode combination. What we are aiming at corresponds to steering the parameterization of the wave function to focus on the most important mode–mode correlations. This idea is very far from new and is inherent in many different approaches. For a similar task, the selection and screening of mode–mode couplings for the parameterization of potential energy surfaces, several schemes have previously been suggested.16–22 The importance of these mode–mode couplings may be estimated based on pair coupling estimates (e.g., root mean square deviations of the two-mode potentials),16,17 higher-order coupling estimates,18 local overlap criteria,21 or focus on certain “active” modes of interest.19 For the parameterization of correlated manymode vibrational wave functions, which is in the focus of the present work, selection and screening in one way or the other have also been widely used and are in the most general sense part of almost all anharmonic wave function calculations. Many authors23–31 have used different treatments of particular configurations to speed up perturbational and variational calculations. The importance of the configurations is in most cases assessed by perturbational measures for the importance of individual configurations.23,24,26,27,29 The alternative block-diagonalization scheme by Rauhut has been argued to be reliable also in case of (almost-)degenerate states.25 Similar schemes also found application for screening of configurations in multi-configurational self-consistent field theory.32 Other methods, such as those of Benoit,33 perform selections on the level of mode–mode couplings. In this work, the Hamiltonian is considered given, and we calculate measures to automatically decide on the wave function parameterization using a threshold noting that in comparison with the above references: (i) We define the wave function parameterization consistently with respect to full mode combinations, rather than single configurations, and integrate this parameterization in VCI and VCC wave functions.

J. Chem. Phys. 142, 144115 (2015)

(ii) The importance measures applied are different and require only one-mode integrals over the VSCF wave function and the chosen Hamiltonian for their evaluation. We recall that the parameterization of the applied potential energy surface is not affected by the approaches proposed in this work. Our importance estimates for the mode–mode correlations in many-mode vibrational wave functions will be derived from upper bounds to the Hamiltonian matrix elements and a perturbational analysis of the VCC wave function, respectively. The perturbational analysis of the VCC equations has previously been employed to derive new VCC methods (VCC[2pt3] in Ref. 34 and VCC[3pt4] in Ref. 35) that provide intermediate steps in terms of accuracy and computational cost. These approximations were, however, not tailored to the individual mode combinations. Another crucial factor for an efficient ansatz for the vibrational wave function is the choice of vibrational coordinates.36–38 The most common choice of coordinates are normal modes. Normal modes, however, are often very delocalized. For larger systems, normal modes cannot generally be assumed to separate into non-interacting local subsets. A large number of other vibrational coordinates have been suggested in different contexts. Among these are internal coordinates, localized modes,25,38–40 optimized coordinates,41–44 and local modes45–53/local monomer modes.54–57 All of the above choices have in common that they tend to be somewhat more local than normal modes. This locality might either be achieved by construction,45–57 by optimization with respect to a localization criterion,25,39 or by optimization of the VSCF energy in an anharmonic force field leading to a partial localization of the optimized coordinates.41–44,58 The locality of the latter is, however, not guaranteed by construction, but has been observed so far. Localized coordinates have first been suggested for analysis of size extensivity breaking in VCI calculations25 or harmonic calculations.39 Recently, they have also found application in VCI38 and VSCF40 calculations. Here, we will focus on the comparison between normal modes and optimized coordinates. Both are rectilinear, and hence enable easy formulation of the kinetic energy operator, in contrast to the curvilinear internal coordinates. They can in all aspects of our calculation be treated similarly. The partial localization found in optimized coordinates so far is automatically obtained when applying an anharmonic potential in the optimization process. The locality of optimized coordinates can therefore be considered to be a “physically reasonable” locality. Since VSCF is exact for uncoupled systems, we expect the variational optimization of the coordinates to lead to minimization of the correlation between the vibrational modes in an overall average sense. The main drawback of optimized coordinates currently is that, to be rigorous, the optimization procedure requires potentials that are invariant upon rotation.43 This holds true for any Taylor expansion, but unfortunately not for the widely used n-mode representation5 of the potential energy surface [cf. Eq. (4)]. Moreover, the determination of optimized coordinates requires the availability of an invariant potential energy surface prior to the optimization of the coordinates. In the present scheme, this precludes computational savings in the generation of the potential energy surface in optimized coordinates based on the local character of these coordinates.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 15 Apr 2015 09:40:25

144115-3

C. König and O. Christiansen

J. Chem. Phys. 142, 144115 (2015)

This article is organized as follows. Section II starts with a very brief introduction to the methods of vibrational structure theory applied in this work. These will then be used to derive Hamiltonian-based as well as perturbative importance measures for mode–mode correlations in Secs. II A and II B. After stating the computational methods in Sec. III, we validate the importance measures for the examples of formaldehyde and furan in Sec. IV. Besides evaluations of the measures itself (Sec. IV A), we will report zero-point energies (Sec. IV B) as well as infrared (IR) spectra (Sec. IV C) obtained with both VCI and VCC parameterizations. Subsequently, the performance of the proposed screening of the mode combinations for correlated many-mode vibrational wave functions in normal modes is compared to that in optimized coordinates for the water dimer, trimer and hexamer in Sec. V. Finally, we conclude on our results in Sec. VI. II. THEORY

m

that minimizes the expectation value ⟨Φi |H | Φi⟩ in a given basis of one-mode functions. Here, i is a composite index containing the levels i m for each mode in the reference state. The VSCF procedure yields, besides optimized occupied modals, a number of unoccupied, virtual modals, which can subsequently be used to expand correlated vibrational wave functions. In analogy to the product form of the wave function, we can choose to express the individual terms (Htm H ) of the Hamiltonian as a product of one-mode terms (h t m ) for mode m. This leads to a sum-over-product form of the Hamiltonian,61,62   ht m . (4) cm H , t H mH = {t |m t =m H }

Vibrational-structure theory aims at solving the vibrational Schrödinger equation within the Born–Oppenheimer approximation. It requires the knowledge of the potential energy surface of the system. The generation of the potential energy surface is itself a major challenge to current computational methods but not in the focus of the present work. Often, the potential energy surfaces are chosen to be in the so-called nmode representation,5 where we shall use the formulation of Refs. 59 and 60. It relies on a hierarchical expansion of the potential energy surface in terms of sub-potential energy surfaces of reduced dimensionality where the contributions already covered by lower-dimensional terms have been subtracted,  V= V¯ m. (1) m∈MCR[V ]

The corrected sub-potential energy surfaces corresponding to a certain mode combination m are also denoted “bar-potentials” (V¯ m). Here, the term “mode-combination” denotes a set of modes contained in the corresponding sub-potential. The set of all mode combinations considered, i.e., the MCR, defines the level of approximation. In case that all possible mode combinations are included, the n-mode representation is exact. Using this formulation of the potential energy surface, the vibrational Hamiltonian can similarly be expanded in terms of its MCR as    H= H mH = cm H , t Htm H , m H ∈MCR[H ]

VSCF method, aiming at the Hartree product, i.e., the product of one-mode functions ({φimm }, modals) for all modes (qm ),  Φi = φimm (qm ), (3)

m H ∈MCR[H ] {t |m t =m H }

(2) where m H denotes a mode combination contained in the representation of the Hamilton operator. The Hamiltonian for the mode combination m H can further be expanded as a sum over all operator terms (t) being dependent on those modes present in m H . In other words, in order to contribute to H m H , the mode combination mt of the operator term t needs to be equal to m H . This restriction is abbreviated by {t|mt = m H }. cm H , t represents the expansion coefficient for the term t with the corresponding potential function Htm H being a function of all modes in m H . The goal is to find the best approximation to the eigenvalue problem for this Hamiltonian. The most basic approach is the

m ∈m H

For instance, using a n-mode representation of the potential energy surface with a polynomial expansion up to third order in the overall exponent, the operator term H m H , for m H = (a, b), employs the one-mode operators a, a2, b, and b2. These are then combined to the product terms a·b, a2·b, and a·b2. The weighted sum of these terms then forms H m H . This particular form is advantageous for the calculation of multi-dimensional integrals needed in correlated vibrational structure calculations, as they can be evaluated as a product of one-dimensional integrals. The sum-over-product form of the potential energy surface may be obtained as a numerical representation of an existing analytical function, see, e.g., Refs. 61 and 62. We prefer to generate it directly from electronic structure single-point calculations using Taylor expansions,59 or fitting using static or adaptiveiterative grid methods,60 or combinations thereof. In second quantization (SQ), the one-mode terms in Eq. (4) read63  tm m = (5) hSQ hrt mm, s m arm† m as m, r m, s m

where



hrt mm, s m =

  tm m t m m s . dqm (φrmm (qm ))∗ hFQ φ s m (qm ) = r m hFQ (6)

m m Here, φrmm (qm ) and φ m s m (qm ) are modals for level r and s , tm respectively, of mode m, and hFQ is the contribution of mode m to the operator term t in first quantization (FQ). The creation m and annihilation operators, arm† m and a s m , respectively, occupy m or deoccupy a certain modal (r and s m , respectively) for mode m. In the single Hartree-product ansatz of VSCF, each mode is exposed to the mean field of the other modes. The corresponding one-mode mean-field operator may be written in second-quantization as6,34,35  m† m m† m m ⟨Φi| aim† Fm,i = m ar m HFQ a s m a i m |Φi⟩ ar m a s m r m, s m

=

 r m, s m

⟨Φi|



  m† m iact. armm , HSQ , a sm† m |Φi⟩ ar m a s m + Fm,i , (7)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 15 Apr 2015 09:40:25

144115-4

C. König and O. Christiansen

J. Chem. Phys. 142, 144115 (2015)

where Φi refers to the reference state, and    m′ iact. hit m cm H , t Fm,i = ′ m ′, ,i m H |m

Automatic determination of important mode-mode correlations in many-mode vibrational wave functions.

We introduce new automatic procedures for parameterizing vibrational coupled cluster (VCC) and vibrational configuration interaction wave functions. I...
2MB Sizes 0 Downloads 6 Views