Calculation of static zero-point vibrational averaging corrections and other vibrational curvature contributions to polarizabilities and hyperpolarizabilities using field-induced coordinates

The zero-point static vibrational averaging (ZPVA) correction to the (hyper)polarizability is written in first-order as the sum of two contributions, one involving electric field derivatives, and the other a sum over normal coordinate derivatives of the zero point energy. It is shown that the sum over 3N-6 normal modes can be replaced by a single term using field-induced coordinates (FICs). A computational strategy that takes advantage of this simplification is presented and applied to a typical push-pull polyene NH2-(CH=CH)3-NO2. From the dependence of the first-order ZPVA on the field-dependent equilibrium geometry we also obtain other loworder static and dynamic vibrational curvature contributions to the (hyper)polarizabilities. The entire set of electronic and vibrational terms is partitioned into two different sequences, each of which exhibits rapid initial convergence for NH2(CH=CH)3NO2 at the Hartree-Fock/6-31G+p level. Including electron correlation at the MP2 level and the frequency-dependence of the ZPVA correction is discussed.


I. Introduction
It is now widely recognized [1,2] that vibrational contributions to polarizabilities and hyperpolarizabilities can be of major importance for many molecules, especially those of practical interest as nonlinear optical chromophores. In some cases the vibrational terms are much larger than the corresponding electronic properties. Ab initio calculations of the former have been carried out at different levels of approximation. The simplest level is the double harmonic approximation where the vibrational potential is assumed to be harmonic and the dependence of the electrical properties on nuclear displacements is taken to be linear.
Anharmonic effects are completely accounted for in the general perturbation treatment of Bishop and Kirtman (BK) [3][4][5]. A particular set of low-order perturbation terms constitutes what is known as the nuclear relaxation approximation [6,7]. Here nuclear relaxation (NR) refers to the fact that part of the vibrational effect is associated with the shift in equilibrium nuclear positions caused by an electric field. The NR approximation was originally applied to static properties using a finite field (FF) method [7] and, then, extended to dynamic properties in the infinite optical frequency (ω→∞) limit by Bishop, Hasan and Kirtman (BHK) [8][9][10]. Taking the limit ω→∞ is equivalent to the assumption that (ωa/ω) 2 is neglible compared to unity (ωa is a fundamental vibrational frequency); an assumption that has been found to give reasonably accurate results in the cases examined thus far [11,12]. Finally, a generalization of the NR approximation to include properties calculated at an arbitrary optical frequency, which is based on a combination of BK perturbation theory and the FF procedure, has now been presented [13].
It turns out that the NR approximation includes the leading vibrational term of each square bracket type that occurs in the BK perturbation treatment. From the results found by applying the BHK method [10,14,15] to prototypes for molecules of practical interest, mixed conclusions about the role of anharmonicity in nuclear relaxation have been reached. On the one hand, anharmonicity does not seem to be important for planar π-conjugated oligomers while, on the other hand, it is quite significant for polysilanes [14] and even more so for push-pull π-conjugated molecules [15].
The zero-point vibrational averaging (ZPVA) correction is not included in the NR approximation because it has a different origin. In the BK perturbation treatment this correction is, in fact, considered as an adjunct to the equilibrium electronic property. Nonetheless, this term gives rise to first-(and higher-) order perturbation contributions due to mechanical and electrical anharmonicity. These are the only first-order terms that are not contained in the NR approximation. Although the first-order ZPVA contribution to β [16][17][18] and γ [19][20][21][22] has been investigated for a few small polyatomic molecules, the importance of this contribution for medium-sized and large molecules is unknown owing to extensive computational requirements for its determination [19]. In this paper we present a simplification that makes ab initio Hartree-Fock and correlated calculations of the ZPVA correction much more feasible for systems larger than those previously treated.
Our new procedure for evaluating the static ZPVA correction includes the utilization of field-induced coordinates (FICs), which were introduced recently [10] in implementing BHK finite field NR calculations. In general terms these coordinates are specified by the shift in equilibrium nuclear configuration due to imposing a static electric field on the system. Subsequently, we defined a set of analytical FICs [15] and showed that a small number of them suffice to reproduce exactly the results of a complete 3N-6 normal coordinate treatment of the static and ω→∞ NR effect. Since the required number of FICs does not grow with the size of the system this leads to computational simplification.
We demonstrate here that a situation similar to the above exists for the static first-order ZPVA correction. This correction is the sum of two terms, one due to mechanical anharmonicity and the other to electrical anharmonicity. In Sec. II both are written in terms of the zero-point vibrational energy. The electrical anharmonicity term involves electric field derivatives of the force constants, while the mechanical anharmonicity term utilizes normal coordinate derivatives.
Instead of requiring all 3N-6 normal coordinate derivatives, the derivative with respect to a single FIC turns out to be sufficient for any component of the polarizability or hyperpolarizability. Again, there results a reduction in the computations involved, particularly for large molecules.
The reference geometry for a ZPVA calculation is, normally, taken to be the field-free equilibrium nuclear configuration. However, there has been a renewed interest [23][24][25] in selecting alternative reference points. Although this does not give rise to any computational savings, tests on diatomics [24,25] show a gain in accuracy -albeit quite small -can be achieved by shifting to an optimized effective geometry. In the treatment presented here we retain the conventional choice.
After discussing computational considerations in Sec. III, we present a case study in Sec.
IV of the static longitudinal polarizability and hyperpolarizabilities of the typical donor/acceptor polyene NH2-(CH=CH)3-NO2. This molecule was chosen not only because it illustrates the capabilities of the method but also because the first-order anharmonicity contribution to the NR hyperpolarizabilities is known [15] to be large. Although both Hartree-Fock and MP2 calculations of the ZPVA corrections are feasible with our procedure only Hartree-Fock results are reported for reasons discussed in Sec.III. These results were obtained using several basis sets from 6-31G to 6-31G+pd. Our values provide a first indication of the significance of the ZPVA correction for α, β, and γ in π-conjugated push-pull systems and of the relative importance of mechanical versus electrical anharmonicity.
If one can obtain the first-order ZPVA correction, then FF methods [26] can be employed to determine perturbation terms (called C-ZPVA for reasons that will be evident later) beyond those included in the NR approximation. These C-ZPVA terms contain the next (non-vanishing) higher-order vibrational contribution of each (square bracket) type that occurs in the BK perturbation treatment. Thus, they provide information on the convergence of the perturbation series. Using the FF method of Kirtman, Luis and Bishop (KLB) [26] for the first time, the C-ZPVA terms derived from the first-order ZPVA correction are calculated in Sec. IV. The KLB method produces not only static values but also values for various nonlinear optical (NLO) processes in the ω→∞ approximation. Although the static quantities have been obtained previously for a couple of small polyatomic molecules [19,27], NLO calculations have been reported only for diatomics [28].
In the last section (Sec. V) we present our conclusions and discuss possible future extensions. One of these is to account for the frequency-dependence of the ZPVA correction. In that regard two different scaling approaches are suggested. We also speculate on further computational simplifications of the procedure given here.

II. Theory and computational considerations
The ZPVA correction to the electronic property P is just the difference between the average value of the property in the ground vibrational state and the value at the reference geometry, which is taken here to be the field-free equilibrium structure R0 e e zpva P P P − = (1) As noted in the Introduction the perturbation treatment of P zpva yields two first-order terms that may be written as [26,29]: In Eqs.
Again all quantities in Eq.(5) are evaluated at R0 or, equivalently, Q = 0. Similarly, one can write Eq. (4) in terms of electric field derivatives of E zp : where n = 1 for the dipole moment, µ ; n = 2 for the linear polarizability, α ; n = 3 for the first hyperpolarizability, β ; and n = 4 for the second hyperpolarizability, γ . For convenience, the designation of the components of the electric field, F, has been suppressed. Eqs. (5)-(6) are more compact than Eqs. (3)-(4), in either case, however, one must compute the complete Hessian in order to evaluate the term [P] 0,1 or [P] 1,0 . That is because We define the FIC associated with the property P by the relation - If P e is the dipole moment then, in the notation of our previous paper [15], where the subscript har denotes the harmonic component of the FIC. Similarly, for the polarizability χ P =χ 2, har ab , whereas for the first and second hyperpolarizability χ P =χ 3,har abc and χ P =χ 4,har abcd respectively. There is a different χ P associated with each independent tensor component (a total of 3 for the dipole moment, 6 for the linear polarizability, 10 for the first hyperpolarizability and 15 for the second hyperpolarizability). Following the same procedure as in Ref. [15], one can transform from normal coordinates to the set consisting of with ω χ P being the circular frequency obtained from the Hessian defined by χ P . Eq. (8) has an advantage over Eq. (5) in that, for each component of P e , only one anharmonicity parameter is required rather than the entire set of 3N-6 derivatives ∂E zp /∂Qa. correction differs from the NR contribution in that regard because the latter can be obtained by means of a geometry optimization [7,8] without explicit calculation of the vibrational force constants [10,14].
From the effect of a finite field on the static P zpva one can obtain the static and dynamic (in the ω→∞ approximation) C-ZPVA vibrational contribution P c-zpva . The procedure [26] is exactly analogous to that used in calculating static and dynamic (ω→∞) NR vibrational contributions from the pure electronic properties. In either case one determines the change in property values arising from the relaxation of the equilibrium nuclear geometry due to an applied field. On the other hand, it can also be shown [7,30] that the sum P c-zpva = P zpva + P c-zpva is associated with the shape of the field-dependent vibrational potential (evaluated at the fielddependent equilibrium geometry). That is why P c is usually referred to as the curvature term and why we have chosen the notation P c-zpva = P c -P zpva .
Let us define - as the difference in the ZPVA correction caused by changing from R0 to the field-dependent equilbrium geometry RF. In the latter case the property is determined in the presence of the field.
For each property this difference can be expanded as a power series in the field [26]: ... 6 2 ...
The expansion parameters in Eqs.  (25) in the square bracket perturbation theory notation of BK (the Roman numeral superscript gives the total order in electrical and mechanical anharmonicity). In this notation the subscripts on the rhs are The terms on the rhs in Eqs. (20)- (25) are of exactly the same (square bracket) type as the corresponding NR contributions but they are each two orders of perturbation theory higher; for example, It is easy to show [4] that the terms of intermediate order vanish. Thus, from FF calculations based on Eqs. (11)-(13), we can test the convergence of the square bracket quantities in BK perturbation theory for the first time in a polyatomic molecule. Finally, it is important to note that the non-vanishing ZPVA correction terms also increase by two-orders of perturbation theory in each step -in other words, [P zpva ] II = [P zpva ] IV = ... = 0.

III. Computational considerations
For our initial study we chose the typical push-pull polyene NH2-(CH=CH)3-NO2. Only the dominant longitudinal component of the polarizability and hyperpolarizabilities was considered. A recent investigation [15] of the NR contribution to the longitudinal β and γ found that first-order anharmonic terms are quite important in this molecule. Thus, it is natural to question whether the first-order ZPVA correction might be considerable as well.
In order to determine the FICs from Eq. (7)  for an MP2 calculation except that only ∂µ/∂Qa is available analytically.) An alternative way to obtain the higher-(and lower-) order property derivatives is to start with the analytical force and differentiate with respect to the field as many time as necessary [32]. For ∂γ/∂Qa, however, this method requires fourth derivatives whereas in the above procedure no higher than third derivatives are necessary.
For the evaluation of [P] 1,0 by means of Eq. (6) the key quantity is the field-dependent E zp or, equivalently, the field-dependent Hessian. The latter is determined analytically at the Hartree-Fock level and, in principle, also at the MP2 level in GAUSSIAN98 [31]. In practice, after MP2 calculations were completed we found that the analytical results did not reproduce literature values [18] for H2O. Although, the numerical MP2 treatment was successful in this regard it is much more time-consuming and, in addition, numerical fifth derivatives are required for γ zpva . This is the reason why only Hartree-Fock values are reported for the ZPVA correction. Even with the analytical Hartree-Fock Hessian, numerical fourth derivatives are necessary for γ. An alternative is to differentiate the analytical ∂α/∂Qa but, then, one of the three differentiations would involve the 3N-6 normal coordinates and this is not desirable. One must be careful in evaluating ∂ 4 Faa/∂F 4 to stay within the window of field values where the fourth derivative is stable. In this study it was found that fields of ±0.0050 a.u. and ±0.0100 a.u. were satisfactory.
We checked the stability by repeating our calculations with F = ±0.0040, ±0.0080 a.u.
The advantage of using FICs is that [P] 0,1 is obtained by evaluating the single term in Eq.
(8) rather than a sum over 3N-6 normal coordinates. This term contains ∂E zp /∂χP, which is calculated from Eq. (9) by the method of finite nuclear displacements. The magnitude of the displacement used was 0.04 a.u. and the stability of the derivative was checked by repeating the calculation with the magnitude of the displacement doubled.
For the Hartree-Fock C-ZPVA contribution we need to determine the ZPVA correction at the field-free, and at several field-dependent, equilibrium geometries. The geometry optimizations were carried out for the fields 2 k x F0, with k=0-4 and F0 = 8 x 10 -4 a.u., following the procedure of Ref. [10] where the field-free Eckart conditions are rigorously enforced. Then the coefficients of the power series expansions in Eqs. (14)- (20) were obtained using the Romberg technique [33].
In order to compare P c-zpva with P nr we also calculated the NR vibrational polarizability and hyperpolarizabilities at both the Hartree-Fock and MP2 levels. The same method was employed as in our recent paper on analytical FICs [15]. For π-conjugated molecules, previous correlated treatments of vibrational (hyper)polarizabilities are almost non-existent [34,35] even in the double harmonic approximation.
Four different basis sets were employed, the smallest being 6-31G [36]. The latter has been widely used in studies of the electronic and vibrational hyperpolarizabilities of π-conjugated molecules. We, then, added diffuse p or sp functions on the heavy atoms (N,C,O) with exponents that were optimized by maximizing α. This optimization technique has been used previously, for example, in studies by Kirtman and Hasan [37], Spackman [38], and Werner and Meyer [39]. In practice the exponents were varied manually in steps of 0.01. The p exponents so were chosen to maximize α ; the p exponents were not reoptimized. The 6-31G +pd basis was first introduced by Hurst, et al. [40] to calculate polarizabilities and hyperpolarizabilities of polyacetylene oligomers and has been extensively employed since then. A 6-31G+sp basis was found to be most efficient for polysilane oligomers in Ref. [37]. for the longitudinal component of P e , P zpva , P nr and P c-zpva . Although there is a significant variation with basis set, for the three larger bases the variation is modest. The 6-31G basis seems satisfactory for a qualitative analysis -there are a couple of large differences ( > 50%) from 6-31G+p but these occur only when the contribution of that term to the total property value is quite small. We judge the 6-31G+p basis to be appropriate for a semiquantitative treatment. The maximum difference from 6-31G+sp or 6-31G+pd is about 30% (again, only for small contributions) and, in most instances, it is much less. For this reason we did not undertake finite field C-ZPVA calculations for the +sp or +pd basis.
In examining the convergence of the BK perturbation treatment [3][4][5] we suggest that it is best to look separately at two different sequences which combined give the total property: Here P c-zpva (I) , for example, is used to indicate the fact that this C-ZPVA term is derived from [P zpva ] I (although it only contains contributions higher than first-order). From Table I it is evident that the RHF value of the first-order static ZPVA correction is always much smaller than its static electronic counterpart as expected. The relevant ratios are less than .02 for α, .08 for β , and .07 for γ. Thus, sequence (A) initially converges quite rapidly for NH2-(CH=CH)3-NO2. The same is true of sequence (B) where, in every case, the ratio P c-zpva (I)/P nr is less than 0.12. Finally, we note that the C-ZPVA contribution to each dynamic process is smaller in magnitude than the corresponding static value.
Our MP2 results for the electronic and nuclear relaxation contributions (in a.u.) are presented in Table II. Although correlation has a large effect, the qualitative conclusions we have drawn from the RHF calculations remain unaltered (where appropriate). The difference between the 6-31G and 6-31G+p basis sets is somewhat larger for MP2 -typically, of the order of 50%but that is still sufficient for qualitative purposes. Assuming that the 6-31G+p basis will give semiquantitative accuracy, as it does for RHF, the MP2 calculations using the +sp and +pd bases were omitted for the sake of computational convenience. From the 6-31G+p results we see that the ratio of P nr to the corresponding static P e consistently decreases when correlation is added at the MP2 level. Nonetheless, this ratio remains greater than, or comparable to, unity for all the properties mentioned in connection with the RHF treatment except β nr (-ω; ω, 0)ω→∞.

V. Conclusions and future work
The first-order ZPVA correction to the (hyper)polarizability is the sum of two terms, one due to mechanical anharmonicity, i.e.
[P] 1,0 . We have shown that [P] 1,0 may be written as an nth-order derivative of E zp with respect to an electric field, whereas [P] 0,1 depends upon first derivatives of E zp with respect to vibrational coordinates. By using FICs rather than normal coordinates the sum over normal modes in the expression for [P] 0,1 can be reduced to a single term thereby reducing the computational effort so that calculations can be done on large molecules. Given a method of determining the first-order ZPVA correction one can, then, use the finite field approach of KLB to obtain a first estimate of higher-order vibrational contributions, P c-zpva , not included in the nuclear relaxation term, P nr .
Contributions to both static and dynamic (in the ω→∞ approximation) NLO processes are obtained in this manner.
Hartree-Fock first-order ZPVA calculations have been carried out on the model push-pull molecule NH2-(CH=CH)3-NO2 using GAUSSIAN98 to analytically determine the fielddependent Hessian. An MP2 treatment is feasible as well but awaits successful implementation of the analogous analytical procedure. Starting with the Hessian all subsequent differentiations were done numerically. For comparison purposes the NR vibrational contribution was also computed at both the Hartree-Fock and MP2 levels by means of the FIC procedure that we recently described [15]. In addition to 6-31G several basis sets with optimized polarization functions were also employed.
For NH2-(CH=CH)3-NO2 it was found as expected that P nr is often greater than, or roughly the same as, the corresponding static electronic property. Likewise, P c-zpva is often comparable to, or exceeds, P zpva . On the other hand, one can logically partition the entire set of electronic and vibrational terms that contribute to the (hyper)polarizability into two different sequences each of which shows rapid initial convergence. Finally, it turns out that electrical anharmonicity is much more important than mechanical anharmonicity for P zpva of the molecule studied here.
The treatment we have presented pertains to the static ZPVA correction. One simple method for estimating the frequency-dependence is to assume that the correction scales (multiplicatively) with frequency in exactly the same manner as the (field-free) equilibrium value of the electronic property. The limited data available for small molecules [16,19,20,41,42] yields mixed results with dispersion errors at optical frequencies typically, but not always, less than 10%. Further tests of this approximation should be carried out for larger molecules.
Although Eq. (6) cannot be employed to obtain the frequency-dependent [P] 1,0 , Eq. (8) can be utilized for the frequency-dependent [P] 0,1 . If the latter term makes the dominant contribution to the ZPVA correction this immediately gives a good estimate for the frequency dispersion. Although the opposite is true for the molecule studied here the dispersion of [P] 0,1 might mimic the frequency-dependence of P zpva more accurately than the electronic property does so. This possibility remains to be investigated.
The only correlation method discussed in this paper is MP2. We avoided the more efficient DFT approach because common exchange-correlation functionals do not correctly describe [43] the electric field polarization of quasilinear π-conjugated systems like NH2-(CH=CH)3-NO2. On the other hand, DFT does appear to give an accurate Hessian for such molecules. This suggests that a simplification might be achieved by using the DFT E zp in Eq. (8) for [P] 0,1 . 2.74x10 6 3.85x10 6 3.74x10 6 3. χ . The high accuracy of this procedure was verified for the 6-31G+p basis by carrying out a finite field BHK calculation. χ . The high accuracy of this procedure was verified for the 6-31G basis by carrying out a finite field BHK calculation.