Tensor Decomposition and Vibrational Coupled Cluster Theory

Article pubs.acs.org/JPCA

Tensor Decomposition and Vibrational Coupled Cluster Theory Ian H. Godtliebsen, Bo Thomsen, and Ove Christiansen* Department of Chemistry, Aarhus University, Aarhus, Denmark S Supporting Information *

ABSTRACT: The use of tensor decomposition in the calculation of anharmonic vibrational wave functions is discussed. The correlation amplitudes of vibrational coupled cluster (VCC) and vibrational configuration interaction (VCI) theories are considered as tensors and decomposed. A pilot code is implemented allowing a numerical study of the performance of the canonical decomposition/parallel factors (CP) for three and higher mode couplings in computations on water, formaldehyde, and 1,2,5-thiadiazole. The results show that there is a significant perspective in applying tensor decomposition in the context of anharmonic vibrational wave functions, with the CP tensor decomposition providing compression of data and a computational convenient representation. The calculations also illustrate how the multiplicative separability of the VCC ansatz with respect to noninteracting degrees of freedom goes well together with a tensor decomposition approach. Tensor decomposition opens for adjusting the computational effort spent on a particular mode-coupling according to the significance of that particular coupling, which is guaranteed to decrease to zero in the case of VCC in the limit of noninteracting subsystems.

1. INTRODUCTION A first-principles description of chemistry necessitates a simultaneous description of both electrons and nuclei using quantum mechanics. Following the spirit of the Born− Oppenheimer approximation we are left with separate but coupled quantum mechanical equations for the electrons and the nuclei. The solution of the electronic problem for a given nuclear structure defines an electronic energy. Considering the electronic energy as a function of the relative nuclear coordinates, we arrive at the concept of a potential energy surface (PES) providing the potential for the motion of the nuclei. Often a rather platonic picture of the nuclear motion is employed in which only the equilibrium molecular structure as defined from the minimum of the PES is considered.1 In the second most widespread view on nuclei in molecules, the nuclei perform harmonic motion in normal coordinates defined from a second order Taylor expansion of the PES. The normalcoordinate harmonic treatment has for many purposes been an excellent model for interpretation of basic features of molecular spectra. However, the PES is in general very complicated and couples all vibrational degrees of freedom, and a low-order Taylor expansion will have only limited scope of validity. To go beyond the above platonic views on molecules, however, one needs to explicitly solve the Schrödinger equation for nuclear motion, which is an immense task even in approximate calculations. The accurate and efficient description of PESs is alone a very difficult task and has been a very active field of research for years with good progress; see refs 2−5 and references therein. Under the assumption that one can solve the electronic problem and represent the PES, we here choose to focus on solving the vibrational Schrödinger equation. The vibrational degrees of freedom are denoted as modes and describe the © XXXX American Chemical Society

relative positions of the nuclei. The vibrational wave function is a many-mode function. In the need to go beyond the harmonic approximation, one of the simplest and most popular approaches is the method of vibrational self-consistent field (VSCF).6,7 In VSCF, each mode is optimized based on the average field of the other modes. To proceed beyond VSCF and this mean-field picture, one needs to account for the correlation between modes. The VSCF method provides the reference state for such more advanced methods including vibrational configuration interaction (VCI)3,8−12 and vibrational coupled cluster (VCC).13 VCI can be described as based on a linear expansion of the wave function in a space consisting of a reference state and states generated from excitations out of the reference. VCC is the exponential parametrization in the same space, which can be rigorously formulated in a manymode second quantization formulation.14 The free parameters of the models can be denoted as correlation amplitudes and relate more or less directly to the weights of other configurations than the VSCF reference state. Numerical studies have shown that VCC provides higher accuracy at the same level of excitation than VCI,15 and even though the equations of VCC are significantly more involved than those of VCI, due to VCC being nonlinear and nonvariational, recent work has shown how VCC can be implemented with the same computational scaling as VCI.16 Both VCI and VCC were implemented with a well-defined polynomial scaling in the number of modes, M, for a given VCI or VCC calculation Special Issue: Joel M. Bowman Festschrift Received: January 31, 2013 Revised: April 25, 2013

A

dx.doi.org/10.1021/jp401153q | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

such an implementation. It may also inspire similar work in the electronic structure context, as decomposition of correlation amplitudes seems also promising in the electronic structure context. This type of analysis, however, greatly benefits from vibrational structure theory per construction giving a segmentation of the full set of correlation amplitudes into much smaller segments referring to specific indistinguishable modes. Thus, with M being the number of degrees of freedom, the full set of three mode couplings consists of (M(M − 1)(M − 2))/6 smaller sets of couplings, while for electrons in canonical formulations there is only one huge set of threeelectron excitations. In section 2, we outline the basic theory of VCI, VCC, and CP tensor decompositions, as well as present some initial considerations on the advantages of tensor decomposition with respect to compression of data and faster contraction schemes. Section 3 reports results for formaldehyde and 1,2,5-thiadiazole. Also, calculations on separated water molecules are presented to showcase the interplay of tensor decomposition and size extensivity. Finally, we conclude in section 4 with a summary and outlook.

defined in terms of the level of mode-coupling included in the wave function and the Hamiltonian. This should be compared with the exponential scaling of exact theory and straightforward approaches. Even with polynomial-scaling VCI and VCC implementations, the complexity still grows fast with the number of atoms. There is thus still a significant need to be cost-effective to gain increased absolute speed in the computations and perhaps even obtain further reductions in computational scaling. One of the main reasons for the high cost of VCI and VCC calculations is the vast number of correlation amplitudes needed to describe the higher order mode couplings of large molecules. To make VCC more computationally attractive it is thus imperative that we find ways to (i) reduce the amount of data required for describing higher-order mode−mode correlations and (ii) process them more efficiently. In recent years, the decomposition of tensors has become a subject of great interest in many different fields of science; see, e.g., ref 17 and references therein. This is due not only to the compressive power of tensor decompositions, minimizing the amount of data, but also to the inherent way tensor decompositions decouple data sets making for easier analysis. A number of different tensor decompositions exist, including but not limited to canonical decomposition/parallel factors (CP),18,19 Tucker decomposition,20 and the recent Tensor Trains (TT).21 Both CP- and Tucker-decomposition closely resembles the singular value decomposition (SVD) of matrices and can in that regard be viewed as higher-order singular value decompositions (HOSVD). Tensor decomposition techniques are now also being explored in various forms for the quantum treatment of electronic structure.22−24 In the electronic context, the decomposed quantities so far seem to be primarily the twoelectron integrals, an issue that is not so closely related to our context. However, with some additional steps, such decompositions also carry over to correlation amplitudes.22,24 We are not yet aware of studies applying tensor decompositions in the context of anharmonic wave functions or other types of molecular quantum dynamics. In this work, we want to investigate whether it can be fruitful to consider higher order correlation amplitudes as multidimensional tensors and use tensor decomposition to identify the most important contributions. It is interesting to see whether tensor decomposition can offer the flexibility to adjust the treatment of mode−mode correlation according to the needs for a particular coupling. For example, one could speculate that distant localized modes are less correlated than close-lying modes, and we shall also discuss the limit of noninteracting subsystems. We will study these aspects by applying CPdecomposition to correlation amplitudes obtained in real-life VCI and VCC calculations and investigate the accuracy obtained for different numerical options. In many wave function calculations, one includes either in a brute force manner a vast number of correlation amplitudes or use of physical insight to select the most important. Tensor decomposition has the potential of allowing calculations to formally include everything as in the brute force manner, but in the practical calculation, the numerical tensor procedures automatically can adjust the effort according to physical importance. Developing an implementation that fully benefits from tensor decomposition in all stages of the computations is a very significant challenge. It is the hope that the present study can provide motivation and important background information for

2. THEORY 2.1. Vibrational Coupled Cluster. Consider a system of M distinguishable degrees of freedom denoted modes and indexed by m. In vibrational coupled cluster theory (VCC), the wave function for such a system is given by an exponential ansatz |VCC⟩ = exp(T )|VSCF⟩

(1)

using a VSCF state as reference. Here, T is an excitation operator also known as the cluster operator. The VSCF state is defined by having one particular one-mode function (modal) occupied for each mode, and this set of occupied modals is optimized during the VSCF procedure. The remaining modals for each mode can be denoted virtual modals and indexed by a and b, while we reserve the index i for the occupied modals. The cluster operator T is given as a sum of excitation operators, describing each combination of one-, two-, three-, and higher-mode simultaneous excitations, etc.

Tm

∑

T=

(2)

m ∈ MCR[T ]

The mode combination range (MCR) of T specifies the set of mode combinations (MCs) included in the excitation space, where an MC is simply a set of modes, for example, m = {m0m1m2} for excitations in modes m0, m1, and m2. For a given n-mode excitation, the excitation operator can be written as a sum of concrete modal excitations with corresponding amplitudes. Tm =

∑ t μm τμm m

m

μm

(3)

Here, μm denotes a specific excitation in the MC m. For example, for the three mode excitation m = {m0m1m2}, μm is a compound index defining the three virtual modals excited to, a0, a1, and a2. The time-independent Schrödinger equation for the VCC wave function is given by H |VCC⟩ = E VCC|VCC⟩ B

(4)



Article

the amplitudes is as a vector. The complete set of amplitudes in vectorized form t can for a VCC[3] wave function be written as shown in Figure 1, where we have defined tm as the vector holding all amplitudes corresponding to excitations in the mode combination m.

Transforming this with exp(−T) and projecting onto the reference state, we get the VCC energy equation E VCC = ⟨VSCF|exp( −T )H exp(T )|VSCF⟩ = ⟨VSCF|H |VCC⟩

(5)

Similarly projecting the transformed equation against the manifold of excitations included, we get the VCC amplitude equations m em μm = ⟨μ |exp( − T )H exp(T )|VSCF⟩

= ⟨μm |exp( −T )H |VCC⟩ = 0

(6)

The VCC wave function is found by solving these nonlinear equations to determine the VCC amplitudes such that the VCC error vector em μm is zero within a threshold. The definition of the MCR[T] used in eq 2 allows for the use of very general excitation spaces. However, a simple and systematic way of including excitations is to include all excitations involving up to n or less modes simultaneously. This produces a hierarchy of models VCC[n], where n is the highest number of modes included in an excitation. We will in particular focus on VCC[3], but the concepts and implementation are general. 2.2. VCI. The ansatz for the vibrational configuration interaction (VCI) wave function, is a linear combination of a reference state, usually a VSCF state, and excitations out of the reference |VCI⟩ = c0|VSCF⟩ + C|VSCF⟩

Figure 1. Vector packing of VCC[3] amplitudes using Nn virtual modals for mode mn.

Looking at the three-mode amplitude vectors tm0m1m2, it is quickly realized that these are vectorizations of three index quantities. On that notion, they are easily repacked into thirdorder tensors with each vibrational mode in the MCs corresponding to a mode of the tensors. repacking

t m0m1m2 ⎯⎯⎯⎯⎯⎯⎯⎯→ ; m0m1m2

For illustrative purposes, a tensor repacking of this kind is shown in Figure 2, where the virtual modals are numbered 1, 2, 3, etc., reserving modals 0 to be the occupied modal in the VSCF reference. 2.4. Canonical Decomposition/Parallel Factors Tensor Decomposition. The canonical decomposition/parallel factors, abbreviated as CP-decomposition, is one among many possible types of tensor decompositions.17 The CP-decomposition is the factorization of a tensor into a weighted sum of vector outer products. Let ? ∈ I1× I2 × ... × IN be a general Nth order tensor. The CP decomposition of ? is then

(7)

where the operator C is defined as

∑

C=

∑ c μm τμm m

m

m ∈ MCR[C ] μm

(8)

cm μm}

The set of amplitudes {c0, minimization of the VCI energy E VCI =

⟨VCI|H |VCI⟩ ⟨VCI|VCI⟩

are optimized through a

(9)

The minimization of eq 9 can also be reformulated into an eigenvalue equation

HC = CE

(11)

R

?≈

(10)

(2) (N ) ∑ λr a(1) r ◦a r ◦...◦a r r

(12)

where ◦ denotes outer product. This is illustrated in Figure 3 for a third order tensor. The rank of a tensor ? is defined as the smallest number of vector outer products that when summed exactly reproduces the tensor. The rank of ? is denoted R? and is per definition equal to the smallest number of terms that make eq 12 become an equality. The CP-decomposition with R = R? is also known as the rank decomposition of ? . The notion of tensor rank is comparable to that of matrix rank, but the properties are quite different. There is, for example, no direct algorithm to determine the rank of a tensor except for few special tensors, and the problem of finding tensor rank is in general an NP-hard one.25 In practice, finding the rank of a tensor is done by fitting higher and higher rank decompositions until the relative fit is below a given fit threshold. The fit of a decomposition is usually defined as the Frobenius norm analogue of the difference between the original tensor and the decomposed form. Moving on, we will now introduce the CP decomposition into our wave function ansätze

where C, the eigenvectors of the Hamiltonian, holds the optimal VCI states, and E is a diagonal matrix with the VCI energies. In the following, we will focus on the ground state obtained as the lowest root. To better compare with VCC we choose intermediate normalization such that c0 = 1. 2.3. Correlation Amplitudes As Tensors. A tensor is defined as an N-dimensional array, i.e., a structure of elements with N-indices. We define the number of dimensions as the order of a tensor, and thus, an N-dimensional tensor will have an order of N. Following this definition, a vector would also be a first order tensor, and a matrix a second order tensor. All tensors of order three or higher order are in general known as higher-order tensors but can also be referred to as just tensors. The dimensions of a tensor are also denoted modes. For a general overview of tensor ideas and nomenclature, we refer to ref 17. In the VCC (and VCI) equations discussed previously, we have not been explicit on how to index and store the m amplitudes tm μm (and cμm). A conceptual simple way to think of C



Article

Figure 2. Tensorial packing of 3 mode amplitudes for m = {m0m1m2} with number of virtual modals N0 = N1 = N2 = 4.

Figure 3. CP tensor decomposition of a third-order tensor where for simplicity the λr factors have been absorbed into the vectors.

amplitudes. This can first of all be done by checking the fits of each individual three-mode decomposition

and discuss how to investigate its performance and what the perspectives of tensor decompositions are. 2.5. Tensor Representation of the VCC/VCI Wave Function Amplitudes. With the three- and higher-mode amplitudes given in tensor form, we can apply a CPdecomposition CP decomposition

; m0m1m2 ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→ ; m0m1m2 ,CP

ε ; m0m1m2 ,CP = ||; m0m1m2 ,CP − ; m0m1m2||

(14)

or it can, e.g., be probed through the fit of the complete set of original and decomposed amplitudes ε tCP = ||tCP − t||

(13)

(15)

Here, tCP will also include the nondecomposed two-mode and one-mode amplitudes. We use throughout norms that by default are analogues of the Frobenius norm, i.e., the square root of the sum of all elements squared. One thing is the direct error in the amplitudes but we are also interested in the decomposition effect on other properties of the wave function. Defining TCP as the cluster operator using the decomposed amplitudes tCP, we can now define an approximate VCC wave function as

giving us a set of decomposed amplitudes. In this study, we will explore the extent to which such decompositions can be accurate enough, as well as investigate numerical structures in typical data. For this purpose, we solely investigate already converged VCC (or VCI) amplitudes. Thus, we first run a conventional VCC (or VCI) calculation, and thereafter, we make decompositions of this set of amplitudes. We have thus not fully integrated the tensor decomposition into our optimization and contraction algorithms (see also later discussions). Creating such an implementation will be a topic of future research once it has been demonstrated that it is an attractive option. To measure the quality of the decomposed amplitudes, we need to introduce some error measures. The most straightforward way to quantify this is to look directly at the error in the

|VCCCP⟩ = exp(T CP)|VSCF⟩

(16)

It is now possible to measure the effect of the decomposition on, e.g., the VCC energy. First, we define the energy of the decomposed state as CP E VCC = ⟨VSCF|H |VCCCP⟩

D

(17)



Article

imation, but it should be kept in mind that we do not in any case know the exact analytical form of potentials, etc. As mentioned earlier, the CP-decomposed form of the VCC amplitudes can be used to decrease the scaling of contractions between the Hamiltonian and the amplitudes. With the Hamiltonian in the SOP form given in eq 24, consider now the one mode contraction

in complete analogy with the definition of the VCC energy in eq 5. We can now introduce the absolute and relative error in the energy CP ε Eabs CP = |E VCC − E VCC|

(18)

ε ErelCP = ε Eabs E VCC CP

(19)

VCC

VCC

VCC

( ma00am1...1...ajm...ja...Km−K1 −1 =

Finally, we can also check how far the decomposed state is from being an optimized VCC state by checking the error reintroduced by the decomposition into the VCC error vector em μm. One can argue that if this equation is fulfilled to satisfactory accuracy by an amplitude vector of CP form, we still have a valid set of cluster amplitudes. We write the approximation to the error vector as e μmm,CP = ⟨μm |exp( −T CP)H |VCCCP⟩

(20)

ε cCP = ||cCP − c||

ar(,majj) =

0

1

K − 1)

(26)

∑ hambo j

t jmj

j j

(m )

br , bjj

bj

(27)

where the only elements of ) needed in the contraction are the elements corresponding to the specific mode mj that is contracted. The resulting tensor of the contraction would then in turn be given as

(22)

R

(m0m1...mK −1 =

CP

⟨VCI |H |VCI ⟩ ⟨VCICP|VCICP⟩

∑ λr b(rm )◦b(rm )◦...◦b(rm 0

1

◦a(rmj)◦b(rmj+1)◦...◦

j − 1)

r

b(rmK −1)

(23)

Absolute and relative errors in the VCI energies for the decomposed amplitudes can also be defined as in eqs 18 and 19. 2.6. Tensor Decompositions and Fast Contractions. In this subsection and the following, we will briefly describe some perspectives that make tensor decompositions interesting. First, we shall argue that the simplicity of the CP form can be very attractive from a computational speed perspective. An important step in recent research aiming at applying VCC methods to larger systems has been to develop computer automatized procedures for deriving and implementing the detailed VCC equations with well-defined polynomial computational scaling with respect to the number of modes, M.16 For example, VCC[3] with a three mode Hamiltonian operator has an M 4 computational scaling. In this type of VCC implementation, the numerical costly step is evaluation of many different contractions involving cluster amplitudes (and similar for VCI). Here, representing the tensors by their decomposed form could be very rewarding for computational efficiency. For efficiency, we choose to cast the Hamiltonian in a sum over product (SOP) form m

∑ λr b(rm )◦b(rm )◦...◦b(rm

Recasting ) to this form allows the one-mode contraction to be rewritten as a set of R one-mode vector contractions instead

and by constructing an approximate VCI wave function using the decomposed amplitudes, we can define a corresponding approximate energy according to the VCI energy expression in eq 9.

t

(25)

R

) m0m1...mK −1 =

In the case of a VCI wave function, the three- and highermode amplitudes can in complete analogy to the VCC case be packed in tensor form and decomposed. The error in the VCI amplitudes is given as

∑ ct ∏ hm,t

m m ...m ...m

) a00a1...1 bj...jaK −K1 −1

for 1 ≤ K ≤ M, which is equivalent to the forward contraction used in the fast VCC contraction engine presented in ref 16 . Suppose now we, instead of the full ) , had the CPdecomposed form, thus having it as a sum over R component rank-one tensors

(21)

μ

H=

j j

t jmj

r=1

,CP ε e mm,CP = ||e m − em μm|| μm

CP E VCI =

j

bj

and using this definition we can define

CP

∑ hambo

(28)

If the decomposition in eq 26 is exact, the resulting tensor in eq 28 is also exact, i.e., equal to a tensor calculated through the contraction of the full ) in eq 25. The full contraction in eq 25 would have a steep computational scaling of O(nK+1), whereas doing the onemode contraction with the decomposed form as in eq 27 only scales with a mere O(Rn2). Thus, tensor decomposition has the perspective of being very useful in VCC and VCI implementations. This perspective is a decisive motivation for the present work. However, here we solely investigate the numerical possibilities of a CP-representation of the correlation amplitudes of VCC and VCI. We have thus still to make an implementation that exploits the potential benefits, as such a task is very complicated since the VCC implementation in itself is very complex. 2.7. Tensor Decompositions and Compression of Data. In doing a decomposition of the wave function amplitudes, it is also paramount to ask a question regarding data size, i.e., how does the amount of elements before and after decomposition compare? We define N; m as the number of elements in an amplitude tensor ; m . It is trivially given by the product of the dimensions D ;i m of ; m

(24)

N

N; m =

as this allows for fast computer implementations of vibrational correlation theories.16 This Hamiltonian form is an approx-

∏ D ;i i=1

E

m

(29) dx.doi.org/10.1021/jp401153q | J. Phys. Chem. A XXXX, XXX, XXX−XXX


Article

For the CP-decomposed form ; m,CP, the number of elements needed to be stored is easily shown to be

The above discussion is based on the assumption of a separable Hamiltonian. Note, that obtaining such a Hamiltonian in the vibrational context will be related to the choice of coordinates, meaning that the coordinates used to describe A and B must be local to A and B. For some types of applications like solids and specific forms of the PES, an analysis circumventing this assumption may be useful and relevant, for example, leading to the suggestion of modified VSCF methods for ensuring size-extensiveness in delocalized coordinates for extended systems.27 Our perspective here is that even with local coordinates the size-extensivity issue is nontrivial and important. We will explore this issue later in a comparison of VCI and VCC for noninteracting systems, where we will study the interplay of size-extensivity and tensor decomposition. 2.9. Implementation. The above ideas have been implemented in the Molecular Interactions, Dynamics and Simulations Chemistry Program Package (MidasCpp)28 that includes the necessary PES, integral, and VCI and VCC wave function functionalities. As should be clear from the above, the program flow follows standard VCI and VCC calculations, only after convergence new analysis functions are called for analyzing the correlation amplitudes, their decompositions, and the resulting energies and errors. To construct the CP-decompositions of the amplitude tensors, we use an alternating least-squares scheme (CPALS), employing singular-value decomposition to construct pseudo inverses. To find the approximate rank of an amplitude tensor we use the obvious idea of fitting multiple CPdecompositions with increasing ranks until one is good within a given fit threshold, as noted by, e.g., ref 17. Using this method to find the best lowest rank of a tensor, we denote FINDGOOD-CP. For more details on CP-ALS and FIND-GOODCP algorithms, see Supporting Information. To control the accuracy in absolute fit, the FIND-GOOD-CP procedure needs an input fit threshold (FT). The accuracy and rank of the tensor decomposition is studied as function of FT. If a very loose (large) FT is set, the tensor may ultimately be represented as a zero tensor, as the absolute norm error in doing so is less than FT. This essentially corresponds to screening away the excitations for this particular modecombination. However, if FT is chosen very small, many terms will be retained in the CP sum, eventually compromising efficiency and compression.

N

N; m,CP =

∑ D ;i

m

R; m + R; m (30)

i=1

For successful compression of the data after decomposition, the quantities in eqs 29 and 30 at least need to satisfy N; m,CP ≤ N; m , which is true if ⎢ ∏N D i m ; R ; m ≤ ⎢ N i = 1i ⎢⎣ ∑ D ; m + i=1

⎥ ⎥ 1 ⎥⎦

(31)

where ⌊⌋ denotes flooring to the nearest integer value. Equation 31 also includes the special case where the amount of data is exactly the same before and after decomposition, and we are in general content as long as the data size is not increased. It is obviously attractive if we can obtain a large compression of the data, as well as achieving the data in a form convenient for efficient calculations as described in the previous subsection. 2.8. Size-Extensivity and VCC and VCI Compared. An important issue for any approximate theory is its behavior with increasing number of degrees of freedom. We will in later sections present some numerical studies relating to this aspect, and we shall therefore make a few theoretical points. Consider a system consisting of two noninteracting fragments A and B. The Hamiltonian for such a system is a sum of the Hamiltonian of the two subsystems, HAB = HA + HB. The exact wave function for the compound system satisfies the Schrödinger equation, HAB|AB⟩ = EAB|AB⟩. If we consider the two fragments separately, the exact wave function for subsystem A satisfies HA|A⟩ = EA|A⟩, and similarly for B. In second quantization (SQ) formulation,14 we may write the fragment wave functions |A⟩ = WA|vac⟩ and similarly for B, with WA and WB as wave operators. The wave operators create the appropriate products of modes and make the appropriate linear combination. The exact wave function for the compound system is |AB⟩ = WAB|vac⟩ = WAWB|vac⟩, and the exact energy is EAB = EA + EB. Approximate wave functions can be classified according to whether they satisfy this criteria or not. The essence here is that the wave function is manifestly separable as a product when the Hamiltonian is additively separable. Approaches that have the correct scaling properties with respect to the number of systems are commonly denoted sizeextensive13,26,27 (following the nomenclature of ref 26). In the case of a separable Hamiltonian, VSCF is separable. In continuation, for systems consisting of the two noninteracting subsystems A and B, we will have the cluster operator TAB = TA + TB

3. RESULTS In this section, we present sample calculations for formaldehyde, 1,2,5-thiadiazole, and water. In the remainder of this work, we will simply use thiadiazole to denote 1,2,5-thiadiazole. 3.1. Computational Details. The potential energy surface calculations were all driven by MidasCpp, utilizing different electronic structure programs to do the explicit single-point calculations. The formaldehyde and water PESs are previously calculated surfaces (see ref 29) obtained as fourth-order Taylor expansions limited to up to three-mode couplings. The single points used for the expansion were calculated with DALTON30 using the CCSD(T)/aug-cc-pVTZ and CCSD/d-aug-cc-pVTZ for formaldehyde and water, respectively. For thiadiazole, the PES was calculated using the so-called multiresolution scheme,31 where the generation of single points was dictated by the adaptive density-guided approach (ADGA).32 In this approach, the VSCF vibrational wave

(32)

as the solution of the VCC equations, and therefore, exp(T) = exp(TA)exp(TB). Thus, the exponential ansatz in VCC gives rise to a wave function that is manifestly separable, in accord with exact theory. The VCC theory outlined above is, like its electronic counterpart, not a variational theory. This may in some cases introduce difficulties. However, while VCI is variational, the linear expansion is not manifestly separable for truncated wave function expansions. VCI at a given excitation level will thus lose accuracy with increasing number of systems, or alternatively, if the same accuracy is sought, the excitation level must be increased as the system increases. Therefore, sizeextensivity is crucial for larger systems. F



Article

else is indicated, all calculations in this section are done at the VCC[3] level with this distribution of amplitudes. Figure 5 contains data for different rank-R decompositions for all three-mode tensors in formaldehyde. Turning our attention to the individual three-mode fits, we see an overall improvement of the fit for all decomposed tensors going from a rank-1 approximation up to rank-20. After this the quality of some fits start to fluctuate, but this behavior is anticipated as we are down at CP-ALS convergence level for most of the amplitude tensors. Figure 5 also clearly shows the close relationship between the error in the combined decomposed VCC amplitudes and the individual three-mode decomposition fits. The size of the errors are comparable; so, for example, at a rank-20 decomposition where the largest error in the individual three-mode fits is in the order of 10−7, the amplitude error is also in the order of 10−7. The amplitude error is a bit larger than the individual errors as it accumulates the errors from all the fits. Moving on to the energy errors, we see that these are several orders of magnitude smaller than the error in the amplitudes, as the third-order VCC amplitudes only have influence on the energy through the correlation energy stemming from three-mode couplings of the PES. Altogether, we see that we can obtain good accuracy by representing the amplitude tensors with CP-decompositions, but the poor convergence properties and numerical issues of the CP-ALS algorithm can be a hindrance when pushing for very high accuracy in the fits. To overcome these issues other algorithms might be needed. In any case, one should note that, e.g., the energy errors found for reasonable large ranks (say 20) are far lower than many other sources of error. So far we have not taken into account that the three-mode amplitude tensors may not be tensors of the same rank and therefore should be fitted to different rank approximations for the combined decomposition of the VCC amplitudes to be optimal. In Table 1 is shown a summary of FIND-GOOD-CP calculations employing different FT. The amplitude error corresponds nicely to the FT in all calculations as it should, noting that FT is given for the individual three-mode couplings, and it is the full norm that is given in the table. Note that very tight FT compression of the data set is sadly not achieved, as seen for the result using an FT of 10−8. Still, we have obtained a valid and reasonably accurate CP representation of the same data, which can be convenient in some respects. The results in Table 1 also illustrate that even for relatively loose FT, such as 10−4, quite acceptable errors in the energies are found. To give some further insight into the origin of the mean rank, we give in Table 2 the individual ranks for the 20 three modecouplings, together with characteristic data for each MC. There is a clear correlation between the MC statistics and the rank obtained, where MCs with large norms in general lead to decompositions with larger ranks. However, it is also notable that in some cases MCs with comparable statistics may lead to different rank approximations. This is most clearly seen by comparing MCs 16 and 18. Overall, the varied ranks obtained shows that it is attractive to have a flexible procedure that can adjust the rank according to need. In the current setup, it is simple to achieve with the present FIND-GOOD-CP algorithm since we already calculated the data we want to compress. Obviously, this issue becomes more complicated if the tensor decomposition is to be built into the whole optimization procedure. In the above error measures, we only investigated the ground state energy. It is also relevant to check that other properties

function density is used as an importance weight factor, in an adaptive construction of the PES. All single points were supplied by MOLPRO,33 the one-mode part at CCSD(T)-F12/ VTZ-F12 with frozen core and the two-mode part at CCSD(T)-F12/VDZ-F12 with frozen core. In addition, the PES includes three-mode couplings, extrapolated using MOLPRO gradients at the MP2/VTZ frozen core level, provided at the points in the two-mode surface. For more information on the extrapolation procedure, see ref 34. All vibrational correlation calculations were done in MidasCpp using an optimized ground state VSCF state as reference and using for each mode the energetically lowestlying VSCF one-mode eigenfunctions (modals) as basis functions for the subsequent VCC and VCI calculations. The calculations on formaldehyde and thiadiazole all use 8 such VSCF modals per mode except for the data in Table 3, and the water calculations were carried out with 11 modals per mode. The VSCF states in the thiadiazole calculations were found using a B-spline basis.35 For formaldehyde and water, a basis consisting of the 11 first harmonic oscillator functions obtained from a harmonic approximation of the one-mode potentials was used to optimize the VSCF states. The convergence threshold for the CP-ALS is in all calculations set to 10−10 change in fit, and the relative cutoff threshold for constructing the pseudoinverses is set to 10−12. See Supporting Information for further discussion on these thresholds. 3.2. Formaldehyde. Formaldehyde has 6 one-mode combinations, 15 two-mode combinations, and 20 distinct three-mode combinations. Using 8 VSCF modals to describe each mode means that we have 7 virtual modals per mode as excitation space as one modal is used as the occupied modal. Thus, the three-mode amplitude tensors are third order tensors with 7 entries in each mode. This makes for a total of 343 elements in each tensor. In Figure 4 is shown the distribution of correlation amplitudes stemming from a VCC[3] calculation. We see

Figure 4. Distribution of VCC[3] amplitudes of formaldehyde.

that the three-mode amplitudes are many, and a lot of them are relatively small, which adds to our hopes of representing and compressing them with tensor decomposition while maintaining sufficient accuracy for the final wave function. Except when G



Article

Figure 5. R-rank decompositions of VCC[3] 3-mode amplitudes of formaldehyde. Upper figure shows error in amplitude vector in full line, error in the VCC error-vector in dotted line, error in absolute energy (au) in thick dashed line, and error in relative energy in thin dashed line. Middle figure displays absolute fits of each individual decomposed 3-mode tensor. Lower figure displays whether or not convergence was reached in the CP-ALS procedure after 2000 iterations. Both upper and middle plot uses a logarithmic scale for the y-axis.

Table 1. FIND-GOOD-CP Data for VCC[3] Calculations in Formaldehyde amp. err.a

FT −3

10 10−4 10−5 10−6 10−7 10−8

1.83 2.51 2.97 3.20 3.53 3.51

× × × × × ×

evec err.b

−3

10 10−4 10−5 10−6 10−7 10−8

9.58 1.64 2.69 3.00 3.67 4.07

× × × × × ×

E abs. err. (au)c

−5

10 10−5 10−6 10−7 10−8 10−9

5.49 2.10 1.53 3.10 5.20 7.71

× × × × × ×

E rel. err.d

−8

10 10−10 10−12 10−14 10−16 10−15

2.08 7.95 5.79 1.18 1.97 2.93

× × × × × ×

mine

meanf

maxg

comp.h

0 0 1 3 5 9

0.40 1.55 3.55 7.30 12.65 20.80

1 4 7 16 25 37

YES YES YES YES YES NO

−6

10 10−9 10−11 10−12 10−14 10−13

a

Correlation amplitude error. bError-vector error. cAbsolute energy error. dRelative energy error. eMinimum rank of the set. fMean rank of the set. Maximum rank of the set. hCompression according to eq 31, assuming mean rank as the rank for all decompositions.

g

Table 2. Data for Decomposition of VCC[3] Amplitudes in Each Three-Mode MC of Formaldehyde amplitude statistics MC 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 a

norm 5.40 1.40 2.80 1.71 4.17 9.89 8.16 1.06 9.28 1.32 4.80 1.13 1.25 2.89 9.51 1.04 2.89 2.31 8.72 4.26

× × × × × × × × × × × × × × × × × × × ×

max elem. −5

10 10−5 10−4 10−4 10−4 10−4 10−5 10−3 10−5 10−3 10−4 10−3 10−2 10−4 10−3 10−3 10−3 10−3 10−4 10−4

4.40 9.78 1.75 1.64 4.07 9.18 4.30 1.01 7.48 9.22 4.64 1.07 1.25 2.81 9.50 7.50 2.21 2.28 6.13 3.39

× × × × × × × × × × × × × × × × × × × ×

−5

10 10−6 10−4 10−4 10−4 10−4 10−5 10−3 10−5 10−4 10−4 10−3 10−2 10−4 10−3 10−4 10−3 10−3 10−4 10−4

min elem. 1.00 8.67 8.55 7.36 2.06 8.41 5.66 3.79 4.38 3.25 7.17 2.05 1.35 1.22 7.66 1.93 6.82 5.22 4.59 4.09

× × × × × × × × × × × × × × × × × × × ×

−15

10 10−18 10−14 10−14 10−15 10−12 10−13 10−13 10−14 10−12 10−17 10−13 10−11 10−13 10−12 10−13 10−12 10−12 10−12 10−14

FT

Tensor Decomposition and Vibrational Coupled Cluster Theory

Recommend Documents