Variational Calculation of Vibrational Linear and Nonlinear Optical Properties

A variational approach for reliably calculating vibrational linear and nonlinear optical properties of molecules with large electrical and/or mechanical anharmonicity is introduced. This approach utilizes a self-consistent solution of the vibrational Schrödinger equation for the complete field-dependent potential-energy surface and, then, adds higher-level vibrational correlation corrections as desired. An initial application is made to static properties for three molecules of widely varying anharmonicity using the lowest-level vibrational correlation treatment ͑i.e., vibrational Møller–Plesset perturbation theory͒. Our results indicate when the conventional Bishop– Kirtman perturbation method can be expected to break down and when high-level vibrational correlation methods are likely to be required. Future improvements and extensions are discussed.


I. INTRODUCTION
The importance of the theoretical prediction of nonlinear optical ͑NLO͒ properties can be attributed to potential design of new materials with applications in communications, medicine, optical computers, and holography. 1 At the microscopic level the low-order NLO properties are governed by the first and second hyperpolarizabilities ͑␤ and ␥͒.6][17] They write the total vibrational ͑hyper͒polarizability as the sum of two terms -one is the correction to the electronic property due to zero-point vibrational averaging ͑P zpva ͒ and the other is the remaining pure vibrational term ͑P ͒.In zeroth order the BK perturbation treatment utilizes the harmonic approximation for the field-free potentialenergy surface ͑PES͒ and assumes that the electrical properties depend linearly on the normal modes.Thus, it is not entirely unexpected that, in the case of highly anharmonic systems, the double perturbation series in mechanical and electrical anharmonicity might be poorly behaved.Calculations show this to be so in the case of HF ͑Refs.18 and 19͒ and H 2 O ͑Refs.18͒ dimers.Even for some ordinary -conjugated NLO molecules, such as 1-formyl-6-hydroxyhexa-1,3,5-triene, the perturbation series for the ab initio static vibrational hyperpolarizabilities is ͑at least͒ initially divergent. 13,20,21ssuming that all vibrational frequencies are small compared to optical frequencies, i.e., within the infinite ͑optical͒ frequency approximation, there exists an alternative to the perturbation approach for calculating vibrational contributions [22][23][24][25] to the most important NLO processes.In this alternative finite field ͑FF͒ procedure, the total vibrational hyperpolarizability is expressed as the sum of a nuclear relaxation ͑P nr ͒ and a curvature ͑P c ͒ term.P nr is due to modification of the pure electronic response to a static electric field caused by field-induced relaxation of the equilibrium geometry.P c is the sum of the zero-point vibrational averaging correction P zpva plus a remainder P c−zpva , determined by the field-dependent change in the zero-point contribution that arises from nuclear relaxation.The sum of P nr and P c−zpva gives the BK P . 23Whereas P nr consists of loworder anharmonicity contributions, and P c−zpva contains all the higher-order anharmonicity terms.In this paper our primary focus will be on the latter.
One can readily express P zpva as a perturbation series in electrical and mechanical anharmonicity.Because of computational difficulty in evaluating the required anharmonicity parameters, previous FF calculations have been limited to the first term in that series.The time is now ripe, however, for extending the FF procedure to include all the terms by using modern variational treatments of the vibrational Schrödinger equation which parallel those of electronic structure theory.The work presented here represents the first attempt along that line.
In the particular case of diatomic molecules the total vibrational hyperpolarizability can be obtained either by the seminumerical method of Ingamells et al. 26 based on Numerov-Cooley integration, 27,28 or by combining the finite field 29,30 procedure with the standard numerical Numerov-Cooley method.This approach may be extended to polyatomic molecules using the discrete variable representation 31 ͑DVR͒ method, but it is limited to tetra-or penta-atomics.On the other hand, the vibrational self-consistent field ͑VSCF͒ [32][33][34][35] treatment can be applied to much larger molecules.Although the VSCF treatment includes mode-mode anharmonic coupling only in the mean-field approximation, this approximation can be subsequently improved using variation or perturbation methods.
In this initial investigation we will study static vibrational hyperpolarizabilities obtained from ⌬E zpva by the FF procedure.Our treatment will be based on VSCF, with and without a subsequent correction, suitably modified to include the effect of a static electric field.In the VSCF procedure each mode vibrates in the average potential generated by all other modes.By definition there is no correlation between modes.The simplest method for introducing such correlation is second-order Møller-Plesset perturbation theory ͑VMP2͒, which is the approach that will be followed here.7][38][39] However, in order to avoid confusion with vibrational coupled-cluster methods, 40 we will use the nomenclature VMPn, 41 where n is the order of perturbation theory.5][46] Although such approaches are, in principle, more accurate than VMP2, they have a greater computational cost and will be saved for later investigation.
Our presentation is organized as follows.Section II summarizes the theory and computational considerations.It is followed in Sec.III by a presentation of the VSCF and VMP2 results obtained for H 2 O, HOOH, and HSSH using numerical ab initio PES's.These three molecules cover a wide range of anharmonicity effects from weak to strong.In order to provide a reference point for the analysis of our results full vibrational CI ͑FVCI͒ calculations are also reported.However, in our FVCI calculations the ab initio PES is truncated at the quartic terms in a normal coordinate expansion ͑and, in some cases, further modified͒.For the truncated PES we also make a comparison with the BK perturbation theory.Finally, our conclusions and some future plans are given in Sec.IV.

II. THEORY AND METHODOLOGICAL CONSIDERATIONS
This paper is concerned with pure vibrational effects.Therefore, rotation and rotation-vibration couplings are ignored.In addition, we neglect the mass-dependent terms in the effective vibrational potential of the full Watson Hamiltionian 47,48 as is commonly done.Under these approxi-mations, the vibrational Schrödinger equation for nuclear motion of a molecule in the ground vibrational state can be expressed as where M is the number of vibrational modes ͑i.e., 3N-5 or 3N-6͒, Q i is a normal coordinate, and V is isotopically invariant.In general, anharmonic terms in V͑Q 1 , Q 2 , ... ,Q M ͒ will couple the normal modes and, thus, the exact 0 ͑Q 1 , Q 2 , ... ,Q M ͒ cannot be written as a simple product of single-mode vibrational wave functions.However, in the harmonic approximation we may write for the ground state where i 0 ͑Q i ͒ is a single-mode harmonic-oscillator function, and the corresponding vibrational energy, denoted by E har zpva , is just the sum of individual harmonic-oscillator contributions.
It will be of interest to compare the variational results obtained in this paper with the BK perturbation theory. 1 For that purpose we assume, as in previous work, 49,50 that V͑Q 1 , Q 2 , ... ,Q M ͒ can be expanded as a Taylor series in the normal coordinates and define the cubic and quartic terms to be first-order ͓H ˆ͑1͒ ͔ and second-order ͓H ˆ͑2͒ ͔ perturbations, respectively.In that event, the first-order correction to the ground-state vibrational energy, which is given by ͗ 0 ͑0͒ ͉H ˆ͑1͒ ͉ 0 ͑0͒ ͘, vanishes whereas the second-order correction is given by

͑3͒
The vibrational energy corrected through second order by means of Eq. ͑3͒, with a zeroth-order harmonic approximation, is hereafter denoted as E PT2 zpva .In comparison the exact result obtained from a full space ͑see later for definition͒ FVCI calculation, for a truncated Taylor-series expansion of V͑Q 1 , Q 2 , ... ,Q M ͒ limited to no higher than quartic terms, is written as E FVCI zpva .For the PT2 calculations we use the complete set of harmonic excited states that contributes to the second term of Eq. ͑3͒ ͑i.e., harmonic excited states ␣ with nonzero ͗ 0 ͑0 ͉H ˆ͑1͒ ͉ ␣ ͑0͒ ͒͘.This set is constructed from all 1-, 2-, and 3-mode excitations where the sum of single-mode quantum numbers is either 1 or 3 ͑the single-mode quantum number is 0 for the ground state; 1 for the first excited state; etc.͒.
In the FVCI calculations all the excitations where the sum of single-mode quantum numbers is less than 9 have been included.This restriction on the basis set was dictated by practical considerations.By following the calculations as the maximum value for the sum of single-mode quantum numbers is increased we found that, when the ab initio PES is truncated at quartic terms, the energy diverges for the more anharmonic systems.In order to circumvent this problem the truncated potential was modified as described in Sec.III B.
3][34][35] In the VSCF methodology the vibrational wave function has the same form as in Eq. ͑2͒, but the i 0 ͑Q i ͒ are effective single-mode anharmonic vibrational wave functions obtained by solving the single-mode vibrational Hartree-type equation, where the effective potential is

͑5͒
Since the evaluation of V ¯i 0 ͑Q i ͒ depends on the single-mode wave functions, the VSCF equation ͓Eq.͑4͔͒ must be solved iteratively.The ground-state VSCF energy is, then, given as where the second term on the right-hand side is needed so that the coupling terms are counted just once.The M-dimensional V͑Q 1 , Q 2 , ... ,Q M ͒ can be expressed as a sum of 1-, 2-,…, M-mode coupling terms 36 While the computational cost of evaluating the integrals in Eq. ͑5͒ increases rapidly with the number of modes that are coupled, the corresponding contribution to the energy is expected to decrease.As a result most VSCF calculations take into account only 1-, 2-, and 3-mode coupling terms, [37][38][39] although there are now programs available that will treat coupling up to four modes for small molecules. 42In this paper we consider only the effect of 1-, 2-, and 3-mode coupling terms on the curvature ͑hyper͒polarizabilities.
7][38][39][40][41][42][43] In the VMP treatment the perturbation operator H ͑1͒ is given by The VMP1 energy, which is the sum of the zeroth-and firstorder terms, is just the VSCF energy.Thus, the first correction to the VSCF energy is given by the second-order of the perturbation theory and the total zero-point energy through second order will be denoted as E VMP2 zpva .For testing purposes, in addition to calculations with the "exact" numerical PES ͑obtained at a given level of ab initio electronic structure theory͒, we will also present VSCF and VMP2 calculations using the same truncated analytical PES with modified quartic terms as employed in the FVCI calculations described above.
As discussed in the Introduction we wish to obtain the static P c ͑P = , ␣, ␤, and ␥͒ by evaluating the change in the field-dependent E zpva associated with the shift from field-free to field-dependent equilibrium geometry. 25,51Using the BK square bracket nomenclature the contribution of the zeropoint term to P c ͑P c = P zpva + P c−zpva ͒ may be written as a perturbation series, 17 For convenience, we have omitted the superscript zpva on the right-hand side of Eq. ͑9͒.The roman superscript indicates the total order in anharmonicity ͑electrical ϩmechanical͒ and all even order corrections are zero.On the other hand, the remaining P = P nr + P c−zpva contributions to the static ͑hyper͒polarizabilities are given by 17 The sum of the lowest-order square bracket terms of each type in Eqs.͑10͒-͑12͒ constitutes the nuclear relaxation ͑NR͒ contribution.Thus, using the static ␤ as an example, the NR contribution is and ␤ c−zpva consists of all the other terms in Eq. ͑11͒.Calculations of P c−zpva have been done previously using fieldinduced coordinates [52][53][54] and the FF method of Kirtman, Luis, and Bishop 25 but only ͓P zpva ͔ I or E har zpva was utilized because of computational difficulties.In that event, the FF calculation of P c gives just the lowest-order nonvanishing anharmonic contributions of each square bracket type, e.g., Note the definition of ␤ ͑c−zpva͒͑I͒ ; ␣ ͑c−zpva͒͑I͒ and ␥ ͑c−zpva͒͑I͒ are defined analogously.On the other hand, when the fielddependent E PT2 zpva is employed in the FF method, the nexthighest-order nonvanishing perturbation contributions ͓i.e., ͓P zpva ͔ III + P In principle, the complete P c to all orders can be obtained using the FVCI method.However, the calculations realized in this paper are not exact because a truncated ͑and modified͒ PES was employed for H 2 O ͑HOOH and HSSH͒.Nevertheless, our FVCI treatment includes a large part of the vibrational correlation not present in the VMP2 calculations for the complete PES and, in comparison with PT2, it includes a portion of all higher-order terms.Thus, the VFCI results provide a valuable indication regarding the accuracy of the PT2 and VMP2 treatments as we will see.
In order to evaluate the convergence of the BK perturbation treatment it has been suggested that two different sequences should be examined separately: 13,20 ͑A͒ p e , ͓p zpva ͔ I , ͓p zpva ͔ III ,..., ͑B͒ P nr , P ͑c−zpva͒͑I͒ , P ͑c−zpva͒͑III͒ , ... .When the ͓P zpva ͔ I / P e and P ͑c−zpva͒͑I͒ / P nr ratios are smaller than unity the series are initially convergent.If the ratios are each much less than unity, then this suggests that P e + ͓P zpva ͔ I + P nr + P ͑c−zpva͒͑I͒ may be a good approximation to the complete property value.A very important advantage of the variation ͑variation perturbation in the case of VMP2͒ methods is that an accurate numerical PES can be used which is not approximated by a truncated power series in the normal coordinates.The calculated P c has no order-by-order correspondence with the BK square bracket formulas and the question of convergence with respect to anharmonicity does not arise ͑although, of course, in the MP perturbation treatment the convergence of the MPn series is an issue͒.Thus, it is expected that VSCF and VMP2, combined with a numerical PES, will give better results than the BK perturbation theory when the anharmonicity is large.
Our calculations were limited to the component parallel to the dipole moment of the electrical properties, which was determined by numerical differentiation using field values of ±0.0001, ±0.0002, ±0.0004, ±0.0008, ±0.0016, ±0.0032, ±0.0064, ±0.0128, ±0.0256, and ±0.0512 a.u.The smallest magnitude field that produced a stable derivative was selected using a Romberg method triangle. 55Field-dependent VSCF and VMP2 zero-point energies were computed using a modified version of the GAMESS ab initio quantum chemistry package. 56The field-dependent harmonic, PT2, and FVCI zero-point energies were evaluated with our own code.Cubic and quartic vibrational force constants were computed by our own code using the GAUSSIAN98 suite of programs 57 to obtain the quadratic force constants at the equilibrium and vibrationally displaced geometries.For the HF and MP2 calculations we utilized the basis sets developed by Sadlej. 580][61][62][63] The field-dependent geometry optimizations were done by our own program which rigorously enforces the Eckart conditions. 64here is an uncertainty in the calculated P c that arises because of the numerical differentiations involved.Usually the magnitude of this uncertainty will increase with the order of differentiation and is, therefore, largest for the hyperpolarizabilities.As an example consider the evaluation of ␥ c,FVCI .This requires a numerical fourth derivative ͑with respect to the field͒ of the change in E FVCI zpva due to nuclear relaxation.
The computation of E FVCI zpva , in turn, is based on the quartic PES, which means that numerical second derivatives of the analytical Hessian ͑with respect to nuclear displacements͒ are required.Thus, the overall calculation is affected by numerical errors in sixth-order derivatives.Despite considerable effort to achieve high accuracy in these derivatives the precision of the values obtained in the Romberg procedure is not always as high as we would like.Appropriate error bounds are included in the tables and future plans to alleviate this situation are discussed in the concluding section.
In order to make the calculation of the field-dependent E VSCF zpva and E VMP2 zpva we have linked the VSCF code implemented in GAMESS to the GAMESS subroutines that compute the field-dependent electronic energy.A second key modification of the GAMESS VSCF/VMP2 algorithm was made as described in the following.Sometimes the coupling terms that appear in the SCF effective potential lead to negative vibrational eigenvalues.When this happens the terms are reduced by a scaling factor ͑f Ͻ 1.0͒ in the VSCF GAMESS code to ensure that the eigenvalues are positive.In carrying out the subsequent VMP2 calculation this scaling is maintained, which means that the total V͑Q 1 , Q 2 , ... ,Q M ͒ is not the same as that obtained from the electronic energy calculations.As a result the derivatives with respect to the electric field turn out to be poorly determined.For that reason we correct the perturbation potential in the subsequent VMP2 treatment so that the total result corresponds to the exact PES, rather than the scaled PES.In order to obtain well-determined derivatives we also need to require that f be the same regardless of the field.To do that we choose the smallest f generated by considering all the fields.
In the VSCF GAMESS algorithm V͑Q 1 , Q 2 , ... ,Q M ͒ is expressed in terms of 1-, 2-, and 3-mode couplings in grid form. 36All anharmonic contributions in the pure ͑as opposed to VSCF͒ single-mode potential are included.In this work we have used a 16-point grid, but we have checked that the changes in the results are negligible when a 32-point grid is utilized.The one-dimensional VSCF equation given by Eq. ͑4͒ is solved in the GAMESS program using the collocation method of Yang and Peet. 65

III. RESULTS AND DISCUSSION
The purpose of the calculations discussed below is to validate the variational approach for determining curvature contributions to nonlinear ͑and linear͒ optical properties.It should be noted that our results are the first ever obtained for contributions beyond the lowest ͑nonvanishing͒ order of anharmonicity in any polyatomic molecule.A second goal is to assess the accuracy of the VSCF and VMP2 methods, in particular.We begin by calibrating with the nearly harmonic molecule, H 2 O, and then move on to the more anharmonic cases of HOOH and HSSH.In this paper only static properties are considered and, for convenience, we restrict ourselves to the diagonal component of the ͑hyper͒polarizabilities along the molecular dipole axis.
Our general approach is as follows.Firstly, we check the initial convergence of perturbation series ͑A͒ and ͑B͒.The good behavior of both series is necessary in order to obtain reliable results for vibrational ͑hyper͒polarizabilities using the BK perturbation theory.Secondly, we test the new variational approaches as well as the BK perturbation theory by comparing with full vibrational CI, i.e., FVCI, results.In order to make the FVCI treatment feasible we work with an approximate PES, which is truncated at the quartic terms.However, that is only for testing purposes and it should be emphasized that our final VSCF and VMP2 calculations are not subject to this limitation.Indeed, a point of interest is the magnitude of the contributions that are not captured by the quartic approximation for the PES.

A. H 2 O
The results for H 2 O based on Hartree-Fock ͑HF͒ and MP2 electronic structure calculations are given in Table I for P e , ͓P zpva ͔ I , P nr , and P ͑c−zvpa͒͑I͒ .The importance of these vari-ous lower-order contributions to the ͑hyper͒polarizabilities and the effect of basis set and electron correlation upon them have already been analyzed in detail by several authors. 24,51,64,66Here we simply wish to note that, at both the HF and MP2 levels, the ͓P zpva ͔ I / P e and P ͑c−zvpa͒͑I͒ / P nr ratios are always less than 0.24 and 0.18, respectively.Hence, the perturbation sequences ͑A͒ and ͑B͒ are initially convergent, although ͓P zpva ͔ I and P ͑c−zvpa͒͑I͒ must be taken into account for quantitative accuracy.In Table II we present P c,VSCF and P c,VMP2 results calculated for the exact numerical ab initio PES.This table also contains reference P c,FVCI values obtained for the truncated PES as well as P c,har , P c,PT2 , P c,VSCF , and P c,VMP2 results for that same truncated potential.Here the superscript har, PT2, etc. refers to the method for determining the E zpva from which the electrical property was calculated.Note that P c,har = ͓P zpva ͔ I + P ͑c−zpva͒͑I͒ and P c,PT2 TABLE I. Values of P e , ͓P zpva ͔ I , P nr , and P ͑c−zvpa͒͑I͒ for static as well as for the diagonal component of ␣, ␤, and ␥ in the direction parallel to .The calculations were done using the POL basis set of Ref. 58.All quantities are in a.u.Error bounds were determined as the minimum difference between the final Romberg values for different fields.For water, one would expect P c,har to be a reasonable approximation to the total P c .The truncated PES may be used to assess that proposition.From Table II it can be seen that P c,har differs from P c,FVCI by less than 24% using the truncated HF PES.For the truncated MP2 PES the differences are generally much smaller except for the dipole moment which, itself, is quite small.Although our results show that P c,har gives a reasonable approximation the errors that remain are significant, at least for the truncated HF PES.
Given that the errors in the harmonic approximation for the vibrational curvature properties can be significant, it is of interest to examine P c,PT2 which contains the next-highestorder nonvanishing BK perturbation contributions.Using the PT2 method for the truncated HF PES leads to much reduced errors, the maximum being 9% for ␤.For the MP2 PES the error in ␤ actually increases, in contrast with the other properties, although it is still less than 9%.Overall, the P c,PT2 calculation accounts fairly well for anharmonicity.Even though the anharmonic contributions are substantial, the ratio ͑P c,PT2 -P c,har ͒ P c,har is always less than 0.20, which confirms that the BK perturbation expansion for curvature is initially well behaved.
Next we turn to the variational approach, again using the truncated PES to assess its accuracy.In all instances, then, the P c,VSCF values are improved over P c,har .In fact, the VSCF values are somewhat better overall than those obtained by the PT2 procedure, although the two are quite similar.A further improvement is found upon going from VSCF to VMP2.For H 2 O, VMP2 appears to be an excellent alternative to FVCI since the difference between them is negligible in all cases for the truncated PES.
We have also investigated the influence of terms in the PES that are not included in the quartic power-series expansion.The ability to take these terms into account is a major advantage of the VSCF and VMP2 treatments.A comparison of the 3-mode VMP2 results calculated using the truncated HF PES with those obtained from the exact numerical HF PES shows that the additional terms are roughly of the same importance as vibrational correlation ͑i.e., going from VSCF to VMP2͒.Finally, the effect of 3-mode coupling was considered.Table II shows that the difference between 2-mode and 3-mode coupling is less than 2% for exact PES.The same result is obtained if the truncated PES is used instead.Summing up for the water molecule, ͑i͒ P c,har values constitute a reasonable approximation, ͑ii͒ P c,PT2 and P c,VSCF contain an important part of the anharmonicity effect ͑errors reduced to a maximum of 9%͒, ͑iii͒ P c,VMP2 takes into account most of the remaining ͑after VSCF͒ anharmonicity due to mode-mode coupling ͑error smaller than 1%͒, ͑iv͒ the terms not included in the quartic PES are roughly as important as vibrational correlation, and ͑v͒ the effect of 3-mode coupling is small ͑less than 2%͒.

B. HOOH
The electronic and vibrational contributions to the electrical properties of HOOH computed for the component parallel to the dipole moment ͑i.e., the C 2 axis͒, using both the HF and MP2 PES, are presented in Tables I, III, and IV.It turns out ͑cf.Table I͒ that the vibrational contribution to the hyperpolarizabilities is very important in this case.For both ␤ and ␥ the sequence ͑B͒ is, initially, either divergent or barely converging.This ill behavior of the BK perturbation treatment is also reflected in the ͑P c,PT2 − P c,har ͒ / P c,har ratios ͑cf.Table III͒ obtained using both truncated PES's.For example, in the case of ␥ the value of the latter ratio is 1.6 ͑1.3͒ for the HF ͑MP2͒ potential.
The FVCI energy computed using the PES truncated at quartic terms diverges when the vibrational basis set is systematically increased due to the fact that the truncated PES is  41 That is to say, a large displacement in certain directions causes the energy to go below the equilibrium geometry value.We found a fairly simple, though by no means unique, way around this difficulty for HOOH.The first step was to remove the 3-mode and 4-mode coupling terms from the truncated PES.Then, performing large displacements from the equilibrium geometry, it was determined that the major source of the problem was due to simultaneous motion along three pairs of normal coordinates.In each instance one member of the pair was the symmetric OH stretching mode and the other was either the symmetric or antisymmetric OH bend or the torsion.Hence, we increased the quartic force constants of the type ‫ץ‬ 4 V͑Q͒ / ‫ץ‬Q i 2 ‫ץ‬Q j 2 for these three pairs, which is guaranteed to increase the potential.At the same time the diagonal quartic force constants for all four modes were also increased.For simplicity the additive constant was taken to be the same in all cases.The value we used was 9.42ϫ 10 −3 ͑7.68ϫ 10 −3 ͒ a.u. at the HF ͑MP2͒ level.This is the smallest value that eliminates the divergence and, in that sense, the modification is minimal.
The values for the curvature contribution to the static electrical properties of HOOH obtained using the truncated HF and MP2 PES's, modified as described above, are given in Table IV.The same 2-mode coupling approximation was used as for the truncated ͑unmodified͒ MP2 PES in Table III.Thus, a comparison of results in the two tables provides a measure of the effect due to modifying the quartic force constants.From the VSCF calculations we see that the modification reduces all electrical properties by roughly 1 / 3.Although that is substantial we judge that the modified potential is still suitable for a qualitative assessment of the reliability of the PT2, VSCF, and VMP2 methods for HOOH.The VSCF values turn out to be in reasonable accord with VFCI in all cases except for the second hyperpolarizability determined for the HF PES, which is in error by 59%.On the other hand, there is no such agreement for the harmonic or the PT2 treatment which means that neither can be relied upon.Note that PT2 uniformly gives the wrong sign.From the difference between the VSCF and VMP2 results in Tables III and IV it is clear that vibrational correlation is important.For c and ␣ c the VMP2 method gives a substantial improvement over the VSCF value ͑cf.FVCI in Table IV͒; however, hyperpolarizabilities are generally not improved.We conclude that, in our final calculations with the exact ab initio PES reported in Table III, the VMP2 results for the effect of vibrational correlation should be regarded as meaningful for and ␣, though not ␤ or ␥.
At the VSCF and VMP2 levels the differences in Table III between the truncated PES and the exact PES are small, less than 13%, except for the VSCF second hyperpolarizability obtained from the MP2 PES.For HOOH, the influence of terms not included in the quartic power-series expansion of the potential is clearly less important than vibrational correlation.From Table III it is seen that the effect of 3-mode coupling for the exact HF PES is also small, as in H 2 O, except perhaps for the second hyperpolarizability.This exception may or may not be meaningful since it corresponds to the one property where there is a larger difference between the VSCF/VMP2 and FVCI values in Table IV.We believe that higher-order coupling terms are, for the most part, negligible.At this point we want to reemphasize the advantage of VSCF and VMP2 in that all anharmonicity is included except for coupling terms of high dimensionality which, in our case is M ജ 4.
For the HOOH molecule our major results are ͑i͒ P c,har values are very poor, ͑ii͒ P c,PT2 results are even worse since they give the wrong sign, ͑iii͒ P c,VSCF yields a big improvement and, with the exception of the second hyperpolarizability, accounts for anharmonicity fairly well, ͑iv͒ for c and ␣ c P c,VMP2 includes most of the vibrational correlation correction, but for ␤ c and ␥ c there is no improvement over the VSCF values, ͑v͒ the effect of terms not included in the quartic PES is less important than vibrational correlation, and ͑vi͒ 3-mode coupling is potentially significant only for ␥ c .TABLE IV.Curvature contributions to and to the diagonal component of ␣, ␤, and ␥ in the direction parallel to for HOOH.The zero-point vibrational energy E zpva is given as well.Calculations of the PES are at the HF/POL and MP2/POL levels and all quantities are in a.u.Error bounds were determined by the minimum difference between the final Romberg values for different fields.

C. HSSH
The diagonal component of the electronic and vibrational ͑hyper͒polarizabilities of HSSH in the direction of the dipole moment ͑i.e., the C 2 axis͒ is presented in Tables I, V, and VI.The results for the HF PES in Table I show that in the case of , ␣, and ␥ the electronic contribution is dominant, but for ␤ the vibrational contribution is essential to reproduce the correct sign and order of magnitude.Nevertheless, the ͑A͒ and ͑B͒ sequences are initially convergent for all properties.For the truncated PES the ͑P c,PT2 − P c,har ͒ / P c,har ratios are all less than 0.21 except for ␤, in which case the ratio is 0.79.It should be noted that the ␤ c and ␥ c values for HSSH are three orders of magnitude smaller than those of HOOH.These relatively small values are often accompanied by large uncertainties as seen in some instances in Tables V and VI.
The FVCI energy calculations have the same divergence problems as the preceding system.In order to obtain convergence we follow exactly the same procedure as before.The problematic pair of normal modes and the modified force constants are exactly analogous to those of HOOH.For HSSH the additive constant is 2.05ϫ 10 −3 a.u.
From Table VI we see that HSSH is similar to H 2 O in that the VMP2 method gives accurate results versus the reference FVCI treatment for the truncated and modified HF PES ͑with 2-mode coupling͒ and appears to be consistently better than VSCF.A comparison of VMP2 and VSCF values in Table VI indicates that vibrational correlation is important for ␤ c and ␥ c .The relatively large probable errors in our values make this point somewhat uncertain, but it is reinforced by the data in Table V.Unlike H 2 O and HOOH we find that, in the case of HSSH, PT2 yields distinctly better results than VSCF for ␤ c and ␥ c .This may be seen from the truncated PES calculations by comparing the two methods with FVCI in Table VI and with VMP2 in Table V.The P c,PT2 values, in turn, show a modest improvement over P c,har for and ␣ but, given our precision, we cannot say whether the same is true for ␤ and ␥.Apparently anharmonicity does not play an important role in determining or ␣.On the other hand, the fact that P c,har results for ␤ and ␥ are in fairly close agreement with those obtained by VMP2 for the exact PES ͑cf.Table V͒ is due to an accidental cancellation of anharmonicity contributions.This conclusion follows from the substantially larger differences found between the harmonic results and the VSCF values, which include intramodal anharmonicity.
Finally, the data in Table V establish that terms in the PES which are not included in the quartic power-series expansion have a very small effect on c and ␣ c ͑less than 2%͒ but, in general, these terms have an important effect on the curvature hyperpolarizabilities.This is consistent with the analysis given in the previous paragraph.We also note that the differences between the 2-mode and 3-mode property values are negligible for c and ␣ c ͑less than 1%͒ and not   450± 30ϫ 10 0 320± 80ϫ 10 0 470± 70ϫ 10 0 520± 50ϫ 10 0 a This value is not converged but from the Romberg triangle we know it is between 100ϫ 10 −2 and zero.

204108-8
Torrent-Sucarrat, Luis, and Kirtman very important for ␤ c and ␥ c .This result suggests that the 2-mode results obtained for the exact PES provide a satisfactory treatment of the mode-mode coupling.Summing up for HSSH, ͑i͒ as in the case of H 2 O the VMP2 method gives accurate results, ͑ii͒ for ␤ and ␥ vibrational correlation is important and P c,PT2 is better than P c,VSCF , ͑iii͒ vibrational anharmonicity is insignificant for and ␣, ͑iv͒ the effect of terms not included in the quartic PES is small for and ␣ but generally substantial for ␤ and ␥, and ͑v͒ the 3-mode contribution is either negligible or of minor importance.

IV. CONCLUSIONS
The particular vibrational contribution to dipole moments and ͑hyper͒polarizabilities known as curvature arises from electrical and mechanical anharmonicity.Although this contribution can be quite significant its evaluation from field-dependent zero-point vibrationally averaged properties has, thus far, been limited to a first estimate because of computational difficulties.We have presented a new variational approach that circumvents these difficulties.Our approach utilizes the self-consistent solution of the vibrational Schrödinger equation for the field-dependent PES.Just as in the analogous electronic structure treatment this initial VSCF model can be improved by adding higher-level VMP2, VCI, etc. corrections.Here, in calculations for the exact PES, we limit ourselves to VSCF and VMP2.The fact that the exact PES can be taken into account is a noteworthy advance.In this initial study we consider only static dipole moments and ͑hyper͒polarizabilities, which are determined from the energy property by differentiation with respect to the field.In order to obtain satisfactory electrical field derivatives it is necessary to modify the methodology of the GAMESS quantum chemistry program employed for the VSCF and VMP2 calculations.
Our new approach has been applied to three molecules: H 2 O, HOOH, and HSSH.As expected beforehand anharmonic effects are small for H 2 O and larger for the other two molecules.This provides an opportunity to compare P c,VSCF and P c,VMP2 with each other, and with the P c,har and P c,PT2 values of the BK perturbation theory, for varying degrees of anharmonicity.To our knowledge the calculations reported here represent the first time that vibrational contributions beyond P c,har ͓i.e., beyond ͓P zpva ͔ I and P ͑c−zvpa͒͑I͒ ͔ have ever been obtained for other than a diatomic molecule.For benchmark purposes a full vibrational CI ͑FVCI͒ treatment was also carried out.However, to make such calculations feasible the PES's had to be truncated after the quartic terms in a normal coordinate power-series expansion.Additional modifications were required to obtain convergence with increasing size of the basis set in the case of HOOH and HSSH.Nonetheless by combining all of the information we were able to assess the performance of different levels of theory.We emphasize again that the final VSCF and VMP2 reported here are based on the exact numerical ab initio PES.
For the molecules we have studied the initial convergence of sequences ͑A͒ and ͑B͒ give a reliable indication of the importance of anharmonicity.Thus, the magnitudes of the ratios ͓P zpva ͔ I / P e and P ͑c−zvpa͒͑I͒ / P nr are the determining factors.Accordingly, H 2 O is nearly harmonic for all properties; HOOH is quite anharmonic, particularly for ͑hyper͒polarizabilities; and HSSH lies in between.Consistent with this classification we find that P c,har , which corresponds to the first term in sequences ͑A͒ and ͑B͒, represents a good approximation for all properties of H 2 O as well as for and ␣ of HSSH.In these same instances the PT2 method, which corresponds to the second term in sequences ͑A͒ and ͑B͒, leads to significantly improved results.For these cases the BK perturbation method is a reliable procedure for treating the effect of the anharmonicity.However, the remaining cases involve large anharmonicity and the BK perturbation method breaks down.In highly anharmonic situations the new variational approach is required to ensure reliable values.P c,PT2 values may even have the wrong sign ͑although fortuitous cancellation of large anharmonicity effects can sometimes occur͒.P c,VSCF , which is roughly comparable in accuracy to P c,PT2 when the anharmonicity is weak, turns out to be a reasonable starting point even for large anharmonicity.Nonetheless, it is insufficient in the latter case to obtain quantitative results, for which vibrational correlation is needed.VMP2 is the lowest-level vibrational correlation treatment and it works well as long as the anharmonic coupling terms are not too large as they are for ␤ c and ␥ c of HOOH.In such highly anharmonic cases one must chose a higher-level method such as VCC or VCI.One of our future goals is to develop and implement the necessary methodology to compute vibrational contributions to NLO properties using these higher-level methods.
Generally speaking, the effect of 3-mode coupling turns out to be small, though nonnegligible when anharmonicity is very important ͑e.g., ␥ of HOOH͒.Of course, the molecules considered in this paper are small and that picture may change when larger species having more vibrational degrees of freedom are examined.The effect of terms not included in the quartic PES is roughly comparable to, or smaller than, the effect of vibrational correlation.These terms are automatically taken into account by using a numerical PES as we have done here.However, the efficient extension of current numerical methodology to higher-level vibrational treatments and/or larger systems, with or without truncation schemes, remains an open question to be examined in further work.
One problem that has already been encountered is the imprecision due to numerical noise.In the future we plan to implement the method of Kirtman, Luis, and Bishop ͑KLB͒, 25 for dynamic as well as static properties.Following KLB the zero-point vibrational energy will be replaced with the zero-point vibrational average of , ␣, and ␤.As a result the order of numerical differentiation will be reduced by one, two, and three, respectively.͑M.T.͒ thanks the Generalitat of Catalunya for financial help through CIRIT Project No. FI/01-00699.Another author ͑J.M.L.͒ thanks the Generalitat de Catalunya for financial help through the Gaspar de Portolà program.We are grateful to R. Benny Gerber and Galina M. Chaban for illuminating discussions regarding the VSCF method.

TABLE II .
Curvature contributions to and to the diagonal component of ␣, ␤, and ␥ in the direction parallel to for H 2 O.The zero-point vibrational energy E zpva is given as well.Calculations of the PES are at the HF/POL and MP2/POL levels and all quantities are in a.u.Error bounds were determined as the minimum difference between the final Romberg values for different fields.

TABLE III .
Curvature contributions to and to the diagonal component of ␣, ␤, and ␥ in the direction parallel to for HOOH.The zero-point vibrational energy E zpva is given as well.Calculations of the PES are at the HF/POL and MP2/POL levels and all quantities are in a.u.Error bounds were determined as the minimum difference between the final Romberg values for different fields.

TABLE V .
Curvature contributions to and to the diagonal component of ␣, ␤, and ␥ in the direction parallel to for HSSH.The zero-point vibrational energy E zpva is given as well.Calculations of the PES are at the HF/POL level and all quantities are in a.u.Error bounds were determined as the minimum difference between the final Romberg values for two different fields.

TABLE VI .
Curvature contributions to and to the diagonal component of ␣, ␤, and ␥ in the direction parallel to for HSSH.The zero-point vibrational energy E zpva is given as well.Calculations of the PES are at the HF/POL level and all quantities are in a.u.Error bounds were determined by the minimum difference between the two final Romberg values for two different fields.