Extended Coupled Cluster Approach for Molecular Properties : Study of H 2 O and HF Complexes

In this paper, we study stationary variant of extended coupled-cluster response approach for properties. This has been studied at the singles and doubles approximation using cubic-truncated functional. This approximation has been studied earlier around equilibrium for small molecules. In this paper, efficacy of this approximation has been shown using perturbative arguments. Further we have calculated dipole moments and polarizabilities of weakly interacting dimers of HF, H2O and H2O HF complex. Results of HF and H2O monomers have been presented at the same level for comparison. The results have been compared with experimental results, wherever available and other theoretical results.


Introduction
Coupled-cluster (CC) method [1,2] is known to introduce dynamical correlation in an efficient manner through infinite partial summation of important terms of the many-body perturbation theory.This method has been accepted as the method of choice for the calculation of energies and energy derivatives [3][4][5][6].Multi-reference variant of the method has been developed to describe the difference energies accurately [7][8][9].Non-variational form of the coupled-cluster method [1] is the most traditional form in which the size-extensivity of energies and energy derivatives is transparent.Although the non-variational theory does not have a (2n+1)-rule inherent in it, the (2n+1) rule for the cluster amplitudes can be incorporated by a Z-vector technique [10].This has made the non-variational coupled-cluster method viable for computation of the energy derivatives.Using only one extra set of perturbation-independent variables, solved through a linear equation, the energy derivatives with respect to different modes of perturbation can be obtained [11].This Z-vector technique can also be derived within the context of constrained variational method, as was done by Koch et al [12].Parallel to these developments, Pal and co-workers have also pursued a fully stationary approach [13,14] based on the expectation-value (XCC) [15] and extended coupled-cluster (ECC) functionals [16,17], truncated to a fixed power in amplitudes.Being fully stationary, the method enjoys the advantage of a (2n+1)-rule.We used a truncation scheme of the functionals, based on the fixed power in the amplitudes.We have shown that the double-linked ECC functional is more appropriate compared to the XCC functional and produces fully size-extensive derivatives, in particular, the non-linear properties.The symmetric XCC functional involving a conjugate left vector is shown to be linked after denominator factorization, but this linked functional is an infinite series in cluster amplitudes [18].A fixed power truncation of this functional of this in a finite cluster approximation leads to unlinked terms in the stationary equations and thus the properties or energy derivatives are not extensive.Within the XCC framework, a truncation based on the perturbative order can derive connectivity of properties, as was shown by .On the other hand, the ECC functional, derived by Arponen and co-workers [16,17], uses a bi-orthogonal left vector and is shown by Arponen to have a special double-linking structure.This double linking structure of the bi-orthogonal or ECC functional not only ensures a terminating series in cluster amplitudes, but also provides fully connected stationary equations, even when truncated to a fixed power in amplitudes.The ECC functional contains an additional set of left vectors to be determined and thus it has twice the number of amplitudes to be solved in a coupled manner.For the derivative calculation, it is thus more expensive than the normal coupled-cluster method with the Z-vector solved separately.The differences between standard CCSD and ECCSD energy derivatives have been analyzed recently in the review by Piecuch and Bartlett [19].The linear response of the ECC functional and the stationary equations resulting from the variation of the response functional have been derived by us.Since the response functional is also double-linked, the response equations are connected and thus higher order energy derivatives are fully extensive.The above features helped us establish ECC method as a desirable way to perform variational coupledcluster response for energy derivatives.A de-coupled approximation of ECC method has been recently formulated by Vaval [20] simplifying the problem.However, the de-coupling leads to loss in accuracy.Using a cubic truncation scheme of the coupled form of the ECC functional in singles and doubles approximation (ECCSD) and without taking into account orbital relaxation, we studied molecular properties of various small to medium-sized molecules around the equilibrium.Around the equilibrium, the relaxation of orbitals is not so important for molecular properties.In the absence of the relaxed orbitals, inclusion of singles is important in the linear response study.The effect of relaxation is more pronounced away from equilibrium as has been shown by us [21].Van Voorhis et al [22] have recently studied this functional only for energy calculation using some additional quartic terms in a doubles only approximation.
The objective of the present paper is twofold.We first show generally through perturbative arguments the efficacy of the method for equilibrium or near-equilibrium property study.We will show that the terms missing in the cubic -truncated ECCSD are of minor importance for dipole moment around the equilibrium.We discuss the minimum perturbation order in which higher than cubic terms contribute.Similar analysis is done for polarizability calculation to show the viability of cubic-ECCSD method.We then show that this cubic-ECCSD can describe dipole moments and polarizabilities of molecules possessing high dipole moments.We have chosen weakly bound dimers of HF and H 2 O to study this second objective.These systems have drawn attention following Crawford's suggestions [23] that in gas phase any polar ionic molecule with dipole moment greater than 2.5 D is capable of binding an extra electron in the dipole field.Similar phenomenon is observed for polar nonionic molecules, although it has been observed that to have similar binding (or positive electron affinity), the nonionic molecule must have higher dipole moment.Water dimer [24][25][26] and Hydrogen Fluoride dimer [27][28][29] have been studied quite extensively using different level of theories.Theoretical calculation of the spectra of some of the systems has also been subject of intense investigation in recent years [30][31].There exist an infinite number of bound states for the dimers of HF and H 2 O.The importance of these anionic states has been noted in the context of electron scattering studies [32].Detailed studies correlating vertical electron affinity (EA) of such systems with dipole moments have been carried out in different basis [33].It is well known in this context that while EA is very sensitive to basis set, dipole moment is relatively insensitive.Hence, the dipole moment is an easier probe to observe the dipole bound anions.While the vertical EA obtained by Koopmans' approximation is an underestimation of the exact vertical EA, there is no fixed trend of the dipole moment at the SCF and correlated levels.Hence, the calculation of dipole moment must be done with a highly correlated method.These systems are ideal benchmark systems for studying cubic -ECCSD method.Higher order molecular properties are also important since they indicate the possibility of induced dipole moment in the presence of perturbation field.Thus, these can be used as signatures of enhanced stability of dipole-bound anions in the presence of external field or other molecules [34].Knowledge of higher order properties [35] helps in designing environments to enhance dipole moments and this can, in turn, be used in designing molecular anions with higher stability.
The paper is organized as follows.In Section II, we elaborate the method used for the study.In Section III we present an analysis of the accuracy of the cubic-ECCSD method by discussing the minimum perturbation order at which the correction terms contribute.Section IV presents the results and a discussion on them.

Method
The bi-orthogonal CC functional, also known as the ECC functional, was proposed by Arponen [16,17].Using a bi-orthogonal set of parameters for the left and right vectors, the expectation value of an operator in ECC formalism can be expressed as, where Ψ and Ψ ' are parameterized differently and bear a bi-orthogonal relation.Ω is a linear operator which includes hole particle de-excitation operators.It can be seen that the differentiation of Eq. ( 1) with respect to Ω operator leads to the non-variational CC equations of Cizek [1].However, Arponen and co-workers generalized this to use a full exponential parameterization for Ψ' and a double similarity transformation to arrive at what is known as the ECC functional.The functional after the transformation has the form where, L denotes that the T operator to the right of H is linked or connected to H, and DL indicates the double-linked structure ensuring that the left vertex Σ must either be connected to H or to two different T vertices.The double-linking form of the series ensures that it is a naturally terminating series unlike the expectation-value functional and is thus a more useful functional.The cluster amplitudes of ∑, T operators as well as those of the derivative ∑ (1) and T (1) operators are obtained by a stationary principle.The use of stationary principle has attendant (2n+1) rule simplifying the calculation of non-linear properties.The double-linked nature guarantees that the terms generated by the deletion of a vertex from the energy diagrams remain connected.Thus, the stationary equations resulting from the variation of the functional (2) with respect to the cluster amplitudes are themselves connected.These yield the amplitudes of the cluster operators ∑, T operators.The amplitudes are denoted as σ (0) and t (0) amplitudes.The following set of equations is solved to obtain the cluster amplitudes.
To obtain the derivative amplitudes, one uses an additional set of equations similar to above equations.The derivative energy functional E (1) is obtained from the graphs of the energy functional by replacing either the Hamiltonian or one of the cluster operators by the derivatives of the operators.One can easily see that the derivative energy functional also, by construction, has a double-linked structure.The expression for the first derivative of the energy can be written as, where, the right vertices i.e., T and its derivative T (1) , are explicitly connected to the Hamiltonian derivative or to the Hamiltonian (as the case may be) and similarly the left vertices i.e., ∑ or its derivative ∑ (1) , will either be connected to Hamiltonian derivative or to the Hamiltonian (as the case may be) or to two different T vertices.To obtain the derivative amplitudes, the derivative functional is made stationary with respect to the amplitudes of ∑ and T operators.For example, the first derivative amplitudes of ∑ and T operators are obtained by the solution of the following set of equations.
The double-linked structure of the E (1) leads to the connectivity of the terms of the equations (5a) and (5b) and thus the derivative cluster amplitudes are connected.The resulting higher order properties (up to first hyper-polarizability) are thus fully size-extensive.
In the present work, we have used singles and doubles approximation for the cluster amplitudes without any orbital relaxation.Though the functional is naturally terminating, the actual truncation of the functional is at the total of seventh power of the cluster amplitudes in doubles only approximation.
For practical purpose, we truncated the series of energy and energy derivatives up to cubic in cluster amplitudes for the energy.This truncation has been considered adequate for an accurate representation of non-linear properties in our earlier study.A perturbative analysis is generally a helpful guide to identify the importance of different terms.We will show in the next Section that around equilibrium this approximation is justified by making a perturbative analysis of the terms missing in the cubic-ECCSD approximation.

Cubic-ECCSD Approximation
To study the validity of cubic-ECCSD approximate, We consider the important terms at the next higher degree in the polynomial E 0 .At the quartic approximation, the important terms in the functional arise from (Σ 2 (V T 1 3 ) L ) DL , (Σ 2 (V T 1 2 T 2 ) L ) DL and (Σ 1 (V T 1 3 ) L ) DL in singles and doubles approximation, where V is the two particle perturbation term of H.The constraint of double-linking helps in the reduction of terms.However, we show that these expensive quartic terms are at high perturbation order.Using the restricted Hartree-Fock determinant as the zeroth order wave function (which is a very good zeroth approximation around the equilibrium),one can derive the dominant perturbation order at which the amplitudes of Σ and T operators contribute from the equations determining these amplitudes.The leading terms in the equation of t 2 and σ 2 amplitudes are of the first order and those of the t 1 and σ 1 amplitudes are of the second order as in the standard single reference CC theory.Thus, one can see that the quartic terms in the ECCSD functional are at least seventh or higher order in perturbation.The corresponding terms in the equations of σ 2 / t 2 amplitudes will be only one order lower, as these are derived from deletion of one T 2 /Σ 2 vertex, which is of the first order in perturbation.Similarly, in the equations for one body amplitudes, the contributions of these terms are two orders lower.In case, where the RHF determinant is a good zeroth description, as is the case of equilibrium geometry, these high order terms will not have any significant contribution to the amplitude equations and consequently energy functional.However, away from equilibrium, where there is a high degree of quasi-degeneracy and the single reference approximation breaks down, the RHF determinant is not a good zeroth order approximation to the wave function and thus, the contributions of higher order terms may be significant.However, the cubic ECCSD is expected to work reasonably for geometries around equilibrium.
A similar analysis can be done to the linear response functional of the cubic ECCSD method and the response equations.We first observe that the leading term in eqs.(5a-b) for the one body t 1 (1) amplitudes comes from the integrals of dipole moment operator µ ( which is one body H (1) term of the E (1) functional in non-relaxed approximation).Thus, these amplitudes have zero perturbation order in V.The leading perturbative order contribution to the two body response amplitudes, on the other hand, comes from inhomogeneous term µt 2, obtained as a result of stationarity of the (σ 2 (µt 2 ) L ) DL term in E (1) .Thus, this is at the first order in perturbation V. We note that in our non-relaxed model, there is no term corresponding to the second derivative of H, i.e.H (2) .Quartic terms in ECCSD response functional E (1) arise from the response of the terms mentioned above for the quartic E (0) .This reveals that the most important (Σ 2 (V T 1 (1) T 1 T 2 ) L ) DL term is at least at the fifth order and our analysis shows that the stationarity of this term with respect to Σ 2 / T 2 vertex and T 1 vertex contributes at the fourth and third order respectively.Other quartic terms have contribute at even higher order.Thus, the quartic terms can play more important role in the response equations compared to the calculation of energy.However, if the RHF is a good zeroth order approximation, these are still at sufficiently high order to be of importance.Our results of dipole moment and polarizability of some systems with high dipole moments will demonstrate the confidence in this cubic-ECCSD model.

Results and Discussion
In this Section, we present the results of properties up to the polarizability of water dimer, Hydrogen Fluoride dimer and water -HF complex.All these systems have been well studied at equilibrium ground state structure.There are some attempts to study the dipole moments of these systems, due to the dipole-bound nature of the systems.However, most of the studies are at the SCF or MP2 level.We have made a detailed study of these complexes in two different basis sets.One of the basis sets consists of double-zeta valence (9s5p/4s primitive Gaussian type orbitals contracted to 4s2p/2s) augmented by a set of uncontracted p polarization functions with exponent 0.70 on hydrogen and d polarization functions with exponent of 1.211 and 1.6 on oxygen and fluorine respectively.In addition to this, a diffuse p function was used on oxygen and fluorine with exponents 0.0845 and 0.09.Thus, the final basis consists of 4s3p1d.This is denoted as basis set A. Another basis set, denoted as basis B, is the Sadlej's optimized basis [36] set for property evaluation.This basis consists of contracted [5s3p2d/3s2p] functions.We report the property values of their respective monomers in the same basis sets for comparison, which highlights the extent of enhancement of properties in dimers.The energy values of each system have also been reported, which highlights the quality of the basis.For comparison, experimental values, wherever available, have been quoted.However, these are vibrationally corrected ones.Only for HF monomer, good estimate for the extrapolation of the experimental dipole moment to the non-vibrating limit is available and is reported here [37].For comparison, we have also reported the MP2 values of properties in both the basis sets.Only for HF, where detailed studies have been made, we have reported some CC results available in Sadlej basis.Some general conclusions are in order.We observe that the dipole moment values are reduced consistently by inclusion of electron correlation.Thus the dipole moment at MP2 level is lower than the SCF dipole moment and the dipole moment at ECCSD level is lower than that at MP2 level.In addition, in the more extensive basis B, dipole moment is reduced at each level (i.e SCF, MP2 and ECCSD).However, for polarizability use of more extensive basis set enhances the values at each level.It is also observed that, at MP2 level polarizability values are enhanced compared to SCF values.Further inclusion of electron correlation by ECCSD reduces the polarizability.

Monomer
The water molecule has been chosen to be in the xy plane with y axis as the molecular axis.The geometry used in this paper is the equilibrium geometry of water i.e.R(OH)=0.95732 0A and (∠HOH=104.5 0).Table 1 reports the property values of water molecule.Comparing the dipole moments in two different basis sets we observe that the basis A overestimates the dipole moment (2.102 D) compared to the experimental value (1.855D) 38 .However, in Sadlej's basis (1.860 D) it is in good agreement with that obtained by experiment.For the polarizability values reported in this paper, in the basis set A, the values are underestimated where as in the more extensive basis B, the values are in good agreement with the experimental values.

Dimer
A number of theoretical studies have been reported for the water dimer using ab initio calculations.
For the dimer, several conformations such as linear, bi-furcated and cyclic shapes are local minima.
However, very recently more elaborate calculations by Marsden et al [39] showed that the bifurcated structure is a transition state.The most stable structure of the water dimer is a linear structure.We report the calculations carried out at this structure of the water dimer having a trans linear configuration (C s ) with a linear hydrogen bond as described in ref 40.The R O-O is 2.98 0 A, with proton donor water monomer in XY plane while the proton acceptor water monomer in YZ plane.Table 2 reports property values in both the basis sets for water dimer.The RHF wave function for water dimer in C s symmetry has configuration 1a' 2 2a' 2 3a' 2 4a' 2 1a'' 2 5a' 2 6a' 2 7a' 2 8a' 2 2a'' 2 .We compare our results of the dipole moment with the experimental values [26].First, we note the significant increase in the values of dipole moment and polarizabilities of the dimer in relation to the monomer in a given basis.Comparing dipole moment results in basis A, we observe that the value at ECCSD level (2.803D) is 0.157 D higher than the experimental dipole moment (2.64 D) of the dimer obtained by Dyke et al 26 .At the MP2 level in the same basis the dipole moment (2.830 D) is 0.19D higher compared to the experimental value.However, the more extensive basis of Sadlej (basis B) produces dipole moment of 2.775 D at the SCF level itself.With the inclusion of electron correlation the ECCSD calculation reduces the dipole moment to 2.658D, which is just 0.018D higher compared to the experimental value.Thus, in basis B the dipole moment obtained by ECCSD method is in close agreement with the experimental value.For this system, we have also reported two components of polarizability as indicative values.Although there are no experimental results, we can infer from the experience of the results of monomer that the values of polarizabilities in the basis B are closer to what would be the correct ones.

Monomer
We report the properties of Hydrogen fluoride at its equilibrium geometry (R H…F =0.91696 0 A).Hydrogen fluoride has been well studied for molecular properties in different basis as well as at different geometry (equilibrium as well as at stretched geometry) [5,41].Kondo et al [41] have studied the energy and property surface for HF using orthogonally spin adapted CC response method.We compare our results with the experimental values [37] as well as the finite field results using orthogonally spin adapted CC linear response results [41] in basis B. For HF monomer, however, good estimate of vibrational correction to experimental dipole moment is known.This is more amenable to comparison with our calculated values.We have reported this extrapolated experimental dipole moment.Table 3 reports the properties using basis A and B. In basis A, dipole moment at the SCF level 2.083D is reduced to 1.975D by inclusion of electron correlation.In the basis B, which is an optimized basis for properties, having two polarized functions for fluorine atom, i.e. 2d functions, the SCF value is reduced to 1.923 D. The inclusion of electron correlation reduces it further to 1.798D.Thus, dipole moment in basis B is in close agreement with the experimental value (1.8265 D).It is underestimated by 0.0285D only.Compared to the orthogonally spin adapted linear response results of Kondo et al [41], our value is 0.016D higher towards the experimental value.Observation of the parallel (α zz ) and perpendicular (α xx ) components of polarizability reveals that as in the case of water, both components of polarizability are enhanced by the electron correlation compared to the SCF value.The value of the polarizability is increased, as more polarized functions are added, making it in better agreement with the experimental value.However, we observe that a larger basis is necessary to get reliable values for the perpendicular component of polarizability.In basis A, for example, the parallel polarizability value (4.45 au) is roughly 2 au lower compared to the experimental (6.40 au) value.On the other hand, the perpendicular polarizability (1.82 au) is about 3.26 au less than the experimental value.The error of perpendicular polarizability is more pronounced in terms of percentage error.On the other hand, the perpendicular component of polarizability (5.12 au) improves significantly in basis B showing considerable basis set effect.Thus, the perpendicular polarizability component is more sensitive to basis set compared to the parallel polarizability component.Hydrogen fluoride dimer has been studied in the same two basis sets as its monomer.The geometry used is a linear C s structure described in the reference 42.The distance between the two fluorine atoms as R F..F = 2.7625 0 A, with zero bond angle between them.H-F bond length of the proton donor HF monomer is 0.9205 0 A, while that of the proton acceptor HF monomer is 0.91918 0 A. Thus, the proton donor monomer is bit more stretched (0.0037 0 A) compared to the HF monomer at its equilibrium geometry.The RHF wave function for HF dimer in linear C s symmetry has configuration 1a' 2 2a' 2 3a' 2 4a' 2 5a' 2 6a' 2 1a'' 2 7a' 2 8a' 2 2a'' 2 .In Table 4, we report dipole moment and polarizability in basis A and B at the ECCSD level.In basis A Dipole moment at the ECCSD level is 3.778 D. In the most extensive basis B, the dipole moment value calculated at the ECCSD (3.462 D) level is 0.462 D higher compared to the experimental value.However, in both the basis sets dipole moment value is overestimated compared to the experimental value [28] (3.0 D).Since the geometry used by Dyke et al [28] for the experimental determination of dipole moment of HF dimer is quite different than the one used in this paper, there is discrepancy between the theoretical and experimental dipole moment value.However, the basis set effect is more for the polarizability component perpendicular to the plane of the molecule (α zz ) than for the in-plane components.The effect of correlation increases the in-plane as well as perpendicular components of polarizability.From the monomer results, we assume that the polarizability results obtained by basis F will be closer to the experimental values.

Hydrogen fluoride -Water complex
It is interesting to study the HF-H 2 O complex along with their respective dimers in the same basis sets.The geometry of the complex is the planer equilibrium structure with C 2V symmetry [43] with O-F bond distance to be 2.71 0 A assuming the monomer structures to be unchanged.It has been observed that the complex has larger dipole moment compared to their respective dimers.There is no unique ground state geometry for this complex too.It is well-known that the dipole moment of the HF-H 2 O depends on the angle between Oxygen and Fluorine.We have assumed this angle to be zero in our study.The RHF wave function for HF-water complex in 1 A 1 ground state has configuration 1a 1 2 2a 1 2 3a 1 2 4a 1 2 5b 2 2 6a 1 2 7a 1 2 1b 1 2 2b 2 2 2b 1 2 .Table 5 reports the dipole moment and polarizability values in basis A and B. We compare the dipole moment with the experimental value available [43].We see that the dipole moment at the ECCSD level in basis A is 4.655 D, which is quite high compared to the experimental value (3.82 D).In basis B dipole moment is reduced to 4.395D, which is 0.0575D higher compared to the experimental value.For polarizability, we do not have any experimental values available for comparison.The observations are in similar lines with the previous observations.

Table 1 . Property values for H 2 O in different basis sets a .
b: See Ref. 38.

Table 2 .
Dipole moment and polarizability values for water dimer.

Table 3 .
Properties of Hydrogen Fluoride in different basis.

Table 4 .
Properties of HF dimer in different basis.

Table 5 .
Properties of HF-Water complex.