Simple Finite Field Method for Calculation of Static and Dynamic Vibrational Hyperpolarizabilities: Curvature Contributions

In the static field limit, the vibrational hyperpolarizability consists of two contributions due to: ͑1͒ the shift in the equilibrium geometry ͑known as nuclear relaxation͒, and ͑2͒ the change in the shape of the potential energy surface ͑known as curvature͒. Simple finite field methods have previously been developed for evaluating these static field contributions and also for determining the effect of nuclear relaxation on dynamic vibrational hyperpolarizabilities in the infinite frequency approximation. In this paper the finite field approach is extended to include, within the infinite frequency approximation, the effect of curvature on the major dynamic nonlinear optical processes.


I. INTRODUCTION
The design of materials with large nonlinear optical ͑NLO͒ properties is currently of great interest, 1-3 mainly because of potential utilization in a variety of optical and electro-optical devices.At the molecular level these NLO properties are determined by the dynamic hyperpolarizabilities.Although one might, at first, think that dynamic hyperpolarizabilities are primarily electronic in origin, there is a growing body of evidence that, for materials with large NLO properties, the vibrations play an important role.Indeed, a number of cases exist [4][5][6][7] where the vibrational hyperpolarizability far exceeds the electronic hyperpolarizability.
A perturbation treatment of dynamic vibrational hyperpolarizabilities has been given by Bishop and Kirtman 8,9 ͑BK͒.This treatment is based on the general sum-over-states formulas 10 for the total hyperpolarizability given in terms of vibronic energies and dipole moment matrix elements.The vibrational and electronic contributions are, then, separated by applying a canonical, or clamped nucleus, approximation 11 wherein the electronic and vibrational motions are considered sequentially rather than simultaneously.BK express the resulting vibrational terms using a double perturbation expansion in orders (n,m) of electrical and mechanical anharmonicity, respectively.In low order, i.e., (n,m) ϭ(0,0), ͑1,0͒, ͑0,1͒, ͑1,1͒, and ͑2,0͒, this leads to a set of compact expressions 9 for the dynamic vibrational hyperpolarizabilities.
There have been two impediments to including anharmonicities in dynamic vibrational hyperpolarizability calculations.First, the BK compact formulas were, until now, 30 complete only through first order ͑i.e., nϩmϭ1͒.Probably of more importance, particularly for NLO materials which usually involve large molecules or polymers, the anharmonicity constants that are required are often computationally burdensome to evaluate.
It is possible to circumvent these difficulties in two special cases through the use of finite field ͑FF͒ methods whereby various molecular properties are determined as a function of one or more static applied electric fields.One special case is the vibrational hyperpolarizability in the static limit, 22,24 which satisfies the relation vibrational hyperpolarizabilityϩZPVA ϭnuclear relaxation contribution ϩcurvature contribution.

͑1͒
Here the vibrational hyperpolarizability is the quantity discussed in BK, and ZPVA refers to the zero-point vibrational averaging correction for the electronic hyperpolarizability ͑which depends upon the nuclear coordinates͒.On the righthand side ͑rhs͒ of Eq. ͑1͒, the nuclear relaxation contribution 24 arises from the change in the electronic energy, or dipole moment, due to the field-induced relaxation of the molecular geometry.The origin of the curvature contribution 24 is the change in zero-point vibrational energy caused directly by the field and indirectly by the geometry relaxation.For a diatomic molecule it has already been demonstrated 31 that Eq. ͑1͒ is valid.The ZPVA term, in particular, is part of the curvature contribution.An extension of the proof for diatomics to an arbitrary polyatomic a͒ Permanent address: Institute of Computational Chemistry and Department of Chemistry, University of Girona, 17017 Girona, Catalonia, Spain.
molecule, 30 building on the treatment given recently by Luis et al., 16 is presented in the next paper.The relationship between the alternative approaches implied by the two sides of Eq. ͑1͒ has recently been reviewed. 32he other special case is the infinite frequency limit.In that limit the nuclear relaxation contribution is obtained by considering the effect of field-induced geometry relaxation on the linear polarizability ͑␣͒ and the first hyperpolarizability ͑␤͒, rather than on the electronic energy or dipole moment.As Bishop, Hasan, and Kirtman ͑BHK͒ have shown, 33 this yields the leading perturbation terms of each type 34 in the formulas for the most common NLO processes.In the few examples that have been examined, 35 it has been found that the vibrational NLO property at a typical optical frequency is not very different from the value at the infinite frequency limit.Thus, the latter constitutes a useful approximation.
The primary purpose of the current paper is to extend the FF/infinite frequency method to include the analog of the curvature contribution in Eq. ͑1͒.We find that instead of the electronic ␣ and ␤, which yield the nuclear relaxation contribution, one must now use the ZPVA correction to these properties.Our treatment then follows along exactly the same lines as BHK and gives entirely analogous results, as will be seen from the analysis carried out in the next section.
The leading terms in the ZPVA correction for ␣ and ␤ are first order in electrical or mechanical anharmonicity.However, the required anharmonicity constants have repeated indices and, therefore, can be determined with only slightly more computational effort than is necessary for the harmonic parameters.This is discussed in Sec.III, along with other approximate computational simplifications.In addition, the finite field approach has certain limitations from a theoretical and interpretive point of view.Methods to reduce these limitations are also considered in Sec.III and, finally, we close with a summary of our results.

II. ANALYSIS
We let the equilibrium geometry in an applied electric field, F, be denoted by R F , while P e (FЈ,R F ) is the value of the electronic property P e calculated at R F in the presence of a field FЈ.In the following, F and FЈ will always be the same, although this is not required.The field-dependent vibrationally averaged value of P e is given by

͑3͒
Here R is an arbitrary geometry and ͉0 F ͘ is the field- dependent ground-state vibrational wave function.The BHK FF/nuclear relaxation method is based on the first term on the rhs of Eq. ͑2͒.One defines the difference and expands this quantity as a Taylor series in F. For P e ϭ e , ␣ e , ␤ e , this leads to In Eqs.͑8͒-͑10͒ we have used the standard notation, e.g., ␥(Ϫ ; 1 , 2 , 3 ), to designate the frequencies of the oscillating electric fields ͑in the order F ␣ ,F ␤ ,F ␥ ¯͒ and, as usual, ϭ ͚ i i .The value obtained for each quantity is that at the field-free equilibrium geometry R 0 .Although all the calculations are done with static fields, Eqs.͑9͒ and ͑10͒ yield dynamic NLO properties in the nuclear relaxation infinite frequency (nr/→ϱ) approximation.Thus, ␤ ␣␤␥ nr (Ϫ;,0) →ϱ contributes to the Pockels effect, ␥ ␣␤␥␦ nr (Ϫ;,0,0) →ϱ to the Kerr effect, and ␥ ␣␤␥␦ nr (Ϫ2;,,0) →ϱ to dc-second harmonic generation ͑dc-SHG͒.Analytical expressions for the terms included in the nr/→ϱ approximation can be obtained 16 from a double expansion about (0,R 0 ) of the electronic energy, V(F,Q), in terms of the field vector Fϭ(F x ,F y ,F z ) and the normal coordinate displacements, Q.For a fixed F this double expansion yields R F and, then, subsequent variations of F x ,F y ,F z give the coefficients in Eqs.͑5͒-͑7͒.Exactly the same expressions can be derived from the BK perturbation treatment by taking the lowest order terms of each type that survive after letting become infinite.These terms are listed in Table I of BHK and the connections with the double expansion method are given in Table I of Ref. 16.The two approaches together yield a definitive interpretation of the quantities in Eqs.͑9͒ and ͑10͒.
Next we consider the second, or ZPVA, term in Eq. ͑2͒.In this case application of the finite field method as above yields the vibrational curvature/→ϱ contribution to the various NLO processes.Unlike nuclear relaxation, the curvature contribution is not limited to low orders of perturbation theory.It is convenient, however, to begin with the lowest order terms in Eq. ͑3͒ which are given by ⌬ P ZPVA ϭ͓ P͔ 0,1 ϩ͓ P͔ 1,0 , ͑11͒ with Now we take the difference where ⌬ P ZPVA (0,R 0 ) is defined by Eq. ͑3͒ with Fϭ0, and expand that quantity as a power series in the field͑s͒.This leads to a set of relations analogous to Eqs. ͑5͒-͑7͒ except that one must use ⌬ P ZPVA (F,R F )Ϫ⌬ P ZPVA (0,R 0 ) instead of (⌬ P e ) R F .Similarly, in Eqs.͑8͒-͑10͒ we substitute ⌬ P ZPVA for P e and simultaneously replace ''nr'' by ''curv.''In replacing ''nr'' with ''curv'' it should be understood that the resulting quantity is the difference between the total curvature contribution and that due to ⌬ P ZPVA .Of course, for ␣, ␤, etc., ⌬ P ZPVA vanishes at the infinite frequency limit.
Until recently, 30 the compact BK perturbation formulas of order ͑0,2͒ have not been available and formulas for total order nϩmу3 are still not known.One can, nonetheless, verify Eqs.͑15͒-͑17͒ with recourse to known formulas.If the static perturbation expression is available ͑as it always is͒, then the infinite frequency formula is readily obtained without knowing the general frequency-dependent result.As far as the diagonal tensor components ͑and the mean values͒ are concerned, the static and infinite frequency expressions are identical except for multiplicative constants, which depend only on the NLO process and the type of perturbation term.This means, for instance, that the ͓ 2 ␣͔ →ϱ III term in Eq. ͑16͒ can be found from the corresponding ͓ 2 ␣͔ ϭ0 III , ͓ 2 ␣͔ →ϱ I and ͓ 2 ␣͔ ϭ0 I terms.The next lowest order terms in ⌬ P ZPVA are two orders higher than those given in Eq. ͑11͒, i.e., ⌬ P ZPVA ϭ͓ P͔ I ϩ͓ P͔ III .͑18͒ Thus, the infinite frequency curvature results obtained by the finite field method are correct through the second ͑total͒ order of perturbation theory.It is unlikely that Eq. ͑18͒ will be utilized in the foreseeable future with standard quantum chemistry programs since the evaluation of ͓ P͔ 0,3 requires quintic force constants and ͓ P͔ 3,0 requires the fourth derivative of the electronic property.Without doing any further analysis we know that the additional terms in Eqs.͑15͒-͑17͒ would be of the same type but two orders higher than those already discussed.On the other hand, it may be feasible to evaluate the total ZPVA correction directly by techniques that sample the complete potential energy surface.Taking that idea one step further, one could determine the total electronic property as given by Eq. ͑2͒.Then the entire vibrational hyperpolarizability, including both the nuclear relaxation and curvature contributions, would be obtained at the same time rather than stepwise, as in the current procedure.

III. DISCUSSION
The starting point for a FF/curvature calculation is the ZPVA correction term given, in lowest order, by Eqs.͑11͒-͑13͒.It is important to note that only a subset of the mechanical and electrical anharmonicity constants, i.e., those with repeated indices, are needed for evaluation of ͓ P͔ 0,1 and ͓ P͔ 1,0 .Once the equilibrium geometry and normal coordinates have been determined, then the required anharmonic force constants can be obtained by taking the diagonal second derivatives of the energy gradients ‫ץ‬V/‫ץ‬Q a ϭϪF a ,

͑19͒
From the numerical point of view, it is clear that the computational effort that must be expended to obtain the set of cubic force constants, F abb , is similar to that involved in calculating the set of quadratic force constants, F ab .A similar conclusion applies to analytical differentiation assuming that an appropriate computer code is available.Essentially the same analysis also pertains to the anharmonic electronic property derivatives ‫ץ‬ 2 P/‫ץ‬Q a 2 .Beyond the simplifications discussed above we can also consider various approximations that may reduce the computational effort.One possibility is to carry out preliminary calculations at a lower level than desired in terms of electronic structure method and/or basis set.These preliminary calculations could be used to specify the normal modes and, perhaps, to eliminate some of these modes as being relatively unimportant.An investigation 36 of this approach for the case of nuclear relaxation has given promising initial results.
The finite field approach to dynamic vibrational hyperpolarizabilities has certain limitations that should be borne in mind.First of all, it is only valid in the infinite frequency limit.In that limit the electronic properties and their zeropoint vibrational averages vanish.In order to estimate the ͑zero-field͒ ⌬P ZPVA at an optical frequency, we can scale the static value in the same manner often used in treating electron correlation. 37An example is afforded by a recent treatment 38 of the mean dc-second harmonic generation in methane.The error due to the difference between the frequency-dependent and the static ⌬ P ZPVA was calculated to be 3.5% of the corresponding electronic property at ϭ0.06 a.u.This error grows to 10.1% at ϭ0.10 au, but is reduced by a factor of 2 when the static ⌬ P ZPVA is scaled.From the few other examples available [39][40][41][42] it appears that, if anything, methane corresponds to a worst case scenario.
The error in the infinite frequency approximation for nuclear relaxation has been examined by Bishop and Dalskov. 35For five small molecules they evaluated the mean dc-SHG, THG, and IDRI, as well as the isotropic and anisotropic Kerr effect at the He/Ne laser frequency ( ϭ0.072 a.u.) and compared them with the same properties at →ϱ.In those cases where nuclear relaxation is important, the maximum error ͑NH 3 ; anisotropic Kerr effect͒ due to the infinite frequency approximation ͑i.e., their approximation B͒ turns out to be less than 12% of the ϭ0.072 value.
Another aspect of the finite field approach, which is a limitation on the one hand and an advantage on the other, is the fact that one obtains the entire curvature ͑or nuclear re-laxation͒ effect without knowing the contribution of individual terms ͑cf.Ref. 16͒.For purposes of analysis, however, there are ͑at least͒ three different ways of dividing up the total contribution that could prove useful: ͑1͒ ͓ P͔ 0,1 ϩ͓ P͔ 1,0 ; ͑2͒ individual normal coordinates Q a in Eqs.͑12͒ and ͑13͒; and ͑3͒ individual normal ͑or internal͒ coordinates in the field-dependent geometry optimization.All three of these divisions can be carried out separately or in concert.
In implementing our FF/curvature procedure it should be borne in mind that all the quantities in Eqs.͑12͒, ͑13͒, and ͑19͒ are field dependent.That includes the normal coordinates, Q a , and the vibrational frequencies, a , as well as the forces, F a , and the electronic properties, P. For a given field one can use standard quantum chemistry programs to determine the field-dependent equilibrium geometry and all Q a , a and Pϭ␣ or ␤ ͑note that F a ϭ0͒.Then, with the field fixed, the required derivatives of P and F a can be calculated using finite displacements of the Q a .Finally, b 2 ZPVA , g 2 ZPVA , and g 3 ZPVA are found by varying the field and fitting the quantity in Eq. ͑14͒ to the analog of either Eq.͑6͒ or Eq.͑7͒.
In summary, the combination of this paper with BHK yields a simple, practical FF method for calculating vibrational NLO properties in the infinite frequency approximation.At the lowest level of treatment the results are complete through second-order perturbation theory ͑and include some of the higher-order terms͒, where we refer here to the total order in mechanical and electrical anharmonicity.An analysis of existing data suggests that the infinite frequency approximation will be adequate in most instances.

APPENDIX
In this Appendix we outline the derivation of Eq. ͑15͒; Eqs.͑16͒ and ͑17͒ may be obtained in an analogous fashion.Our starting point is the lowest order expression for the static ⌬␣ ZPVA ϭ⌬␣ ZPVA (0;0) as given by Eqs.͑11͒-͑13͒.In the notation of Ref. 17, where and q 2 ,␣␤ ϭa 12 ,␣␤ /2a 20 .͑A2͒ Instead of evaluating the potential energy derivatives in Eq. ͑A2͒ at Qϭ0, Fϭ0 as done previously, we regard them as a function of the electric field, F, and the field-free normal coordinates, Q.In order to carry out an expansion in F at the field-dependent equilibrium geometry R F , as desired, we follow the same two-step iterative procedure described in Refs.16 and 17: ͑1͒ the stationary condition (‫ץ‬V/‫ץ‬Q) R F ϭ0 is applied to find R F , and, then, ͑2͒ V is determined as a function of F at that geometry.In doing this it is important to correctly take into account the dependence of the vibrational force constants on F as discussed in the following article.Indeed, the same field-dependent unitary transformation that diagonalizes the harmonic force constant matrix must be applied to all the other a nm coefficients.Then, straightforward algebra leads to the first derivative, a 30 i jk a 30 i jl q 2 k,␣␤ q 1 l,␥ /A 20 i j ͬ ,