A systematic and feasible method for computing nuclear contributions to electrical properties of polyatomic molecules

An analytic method to evaluate nuclear contributions to electrical properties of polyatomic molecules is presented. Such contributions control changes induced by an electric field on equilibrium geometry ~nuclear relaxation contribution ! and vibrational motion ~vibrational contribution! of a molecular system. Expressions to compute the nuclear contributions have been derived from a power series expansion of the potential energy. These contributions to the electrical properties are given in terms of energy derivatives with respect to normal coordinates, electric field intensity or both. Only one calculation of such derivatives at the field-free equilibrium geometry is required. To show the useful efficiency of the analytical evaluation of electrical properties ~th so-called AEEP method !, results for calculations on water and pyridine at the SCF/TZ2P and the MP2/TZ2P levels of theory are reported. The results obtained are compared with previous theoretical calculations and with experimental values. © 1997 American Institute of Physics. @S0021-9606 ~97!02129-6#


I. INTRODUCTION
[3][4][5][6][7][8][9][10][11][12] Such properties give the response of a molecule under the effect of an electromagnetic radiation.Addressing only the most important effect of the electric field component, the induced dipole moment can be expanded in a Taylor series, x,y,z x,y,z where the linear response is provided by the polarizability ␣, and the nonlinear terms are the first and second hyperpolarizabilities ␤ and ␥, respectively.Dynamic properties are defined for time-oscillating fields and the static limit is achieved at the limiting case (→0), i.e., for uniform fields.
In this work, only the static limit of electrical properties ͑P, in general͒ is studied.
The validity of the Born-Oppenheimer approximation for computing nonlinear optical properties was assessed by Adamowicz and Bartlett. 13Within the Born-Oppenheimer approximation, electronic and nuclear motions are evaluated separately, so the total energy of a chemical system is yielded by the sum of potential ͑electronic and nuclear repul-sion͒, vibrational, and rotational energies ͑rotational energies are not considered in this paper͒.5][16][17][18][19][20][21][22][23][24][25][26] All these induced changes can be explained in terms of different contributions to the electrical properties, namely, electronic (P el ), nuclear relaxation (P nr ) and vibrational (P vib ).
Experimental evidence for nuclear contributions (P N ϭP nr ϩP vib ) to electrical properties can be found by comparing the Kerr effect with the electric field induced second harmonic generation ͑ESHG͒ data. 27The Shelton and Palubinskas' work shows that nuclear contributions ͑P nr and P vib ͒ are nearly frequency independent, and tend to a constant value in the high frequency limit for highly symmetric molecules like CH 4 , CF 4 , and SF 6 . 28For large polyethylene molecules, relationships between the bond order alternation ͑BOA͒ parameter provided by the nuclear relaxation contribution and the molecular ͑hyper͒polarizabilities have been established both experimentally 9 and theoretically. 29,30n general, theoretical evaluation of nuclear contributions to electrical properties has been done by using either finite field 13,[19][20][21][22][24][25][26] or perturbation theory [31][32][33][34][35][36][37][38][39][40][41] treatments. Both aproaches allow to use post Hartree-Fock levels of theory.The latter not only allows one to calculate static properties, but also is the only reliable approach to obtain frequency-dependent electrical properties. This echnique implying the summation over all electronic and vibrational states ͑SOS͒ method, is either restricted to small polyatomic molecules, 22,[31][32][33][34] or different kinds of truncations must be applied to the SOS expressions to be reliable for polyatomic molecules.[35][36][37][38][39][40][41] Furthermore, the finite field approach, that has been recently applied to evaluate electrical properties of large polymeric species, 24,42,43 can be used only to obtain components of nuclear contributions to static electrical properties which are parallel to the permanent dipole moment of the molecule.[19][20][21]26 Finally, Handy's group 22 have used finite differences of the vibrational energy to calculate the P v contribution for the HF and H 2 O molecules including anharmonicity corrections.
A third alternate way of to evaluating the total electrical properties (PϭP el ϩP N ) exists.From this alternate treat-a͒ Electronic mail: josepm@iqc.udg.esment, one can obtain the molecular properties directly from the wave function that describes a molecular system, and has been used for long time by different authors.Originally, Kern and co-workers [44][45][46] computed nuclear contributions to one-electron properties, like the dipole moment, by including only mechanical anharmonicity corrections to the SCF potential energy of the H 2 O molecule.Werner and Meyer 47 simplified and applied this method to compute nuclear contributions to dipole moments and polarizabilities of various small molecules.Very recently, Russell and Spackman 48 have computed the vibrational averaging of and ␣ for several small polyatomic molecules including both mechanical and electrical anharmonicity.In the early eighties, Pandey and Santry 49 calculated nuclear relaxation and vibrational contributions to , ␣, and ␤ of the CO, HCN, and H 2 O molecules.However, they did not include the mechanical anharmonicity in the potential energy expansion.Rinaldi et al. 50applied also this analytical approach to evaluate ␣ nr , at the semiempirical molecular orbital level.Recently, Castiglioni et al. 29 developed expressions for the nuclear relaxation contributions to the ␣, ␤, and ␥ at the harmonic level, and determined the P nr contributions from experimental ir and Raman bands.In a previous paper, and only for parallel components to the internuclear axis of diatomic molecules, the analytical evaluation of electrical properties was applied to the CO molecule. 23he goal of this paper is to devise a method to evaluate analytically nuclear contributions (P N ) to the static electrical properties, namely, nuclear relaxation P nr and vibrational P vib .Although some studies have dealt with this subject earlier, 29,30,[44][45][46][47][48][49][50][51][52] analytical and complete evaluation of all components of , ␣, ␤, and ␥ tensors for polyatomic molecules is still lacking.In this work, the potential energy of a chemical system is expanded as double power series in terms of both normal coordinates Qϭ(Q 1 ,Q 2 ,...,Q 3NϪ6 ) and field strength vector Fϭ(F x ,F y ,F z ).The method presented in this work is named AEEP ͑analytical evaluation of electrical properties͒ and has been coded and implemented successfully in our laboratory.The AEEP method as presented here shows the following features ͑i͒ it computes complete values of the total electrical properties, ͑ii͒ it calculates all components of the electrical properties tensors from just one computation of energy derivatives evaluated at zero-field equilibrium geometry; ͑iii͒ it incorporates easily the electron correlation and basis set effects; ͑iv͒ it obtains routinely nuclear contributions of large polyatomic molecules; and ͑v͒ it gives clear rules to know all terms needed to evaluate both nuclear relaxation and vibrational contributions.Moreover, one can analyze the results of the AEEP method and compare them with these obtained by more traditional perturbation treatment, 53 to which the AEEP procedure is related.
In Sec.II, we shall present the AEEP method: nuclear contributions are given by coefficients of the power series expansion of the potential energy that are directly related to energy derivatives.In Sec.III, details of the molecular orbital ab initio calculations carried out in this paper are given.In Sec.IV, the AEEP results for the water and pyridine molecules will be presented.Finally, in Sec.V, conclusions of this work will be summarized.

II. ANALYTICAL EVALUATION OF ELECTRICAL PROPERTIES (AEEP) METHOD: NUCLEAR CONTRIBUTIONS TO STATIC ELECTRICAL PROPERTIES
The potential energy of a polyatomic molecule placed under the effect of an electric field is function of both normal coordinates Qϭ(Q 1 ,Q 2 ,...,Q 3NϪ6 ) and field strength vector Fϭ(F x ,F y ,F z ).Thus the potential energy hypersurface of a chemical system can be expressed as a double ͑Q and F͒ power series expansion, x,y,z x,y,z a nm i 1 ...i n ,... j 1 ...

͑2͒
where n refers to displacements along the 3NϪ6 normal coordinates Q i , m refers to changes of the field strength F j , and N is the number of atoms of the molecule.Differentiation of equation two with respect to either normal coordinate displacements, field strength, or both leads to relationships between coefficients of the power series expansion and potential energy derivatives .

͑3͒
All these derivatives of the potential energy hypersurface are evaluated at the equilibrium geometry (Qϭ0) and for the zero field case (Fϭ0).Then, electronic contributions to the electrical properties are defined by

͑5͒
where indexes xyz and i jk refer to field strength vector com-ponents and normal coordinates, respectively.Molecular property derivatives can also be obtained from the expansion coefficients, e.g., The zero subscript stands for all derivatives involved in Eqs.
͑3͒ and ͑6͒ that are evaluated at the field free equilibrium geometry.
In this section, terms up to nϩmр4 have been considered in the power series expansion of the potential energy ͓Eq.͑2͔͒.Such an expansion implies calculation up to fourthorder derivatives, thus requiring the state-of-the-art techniques in the evaluation of energy derivatives. 54,55This expansion also carries inclusion of first-and second-order mechanical anharmonicity ͑a 30 and a 40 terms͒, first-and second-order electrical anharmonicity for dipole moment (a 21 and a 31 terms͒, first-order electrical anharmonicity for polarizability ͑a 22 terms͒, and harmonic approximation for first hyperpolarizability ͑a 13 terms͒.Under these restrictions, the expansion used for the double power series of the potential energy is given by where abcd and i jkl run over field strength vector components (F x ,F y ,F z ) and normal coordinates (Q 1 ,Q 2 ,..., Q 3NϪ6 ), respectively.Using this truncation for the potential energy, it will be shown that a complete evaluation of nuclear relaxation contribution to dipole moment nr , polarizability ␣ nr , and hyperpolarizabilities ␤ nr and ␥ nr , and vibrational contribution to dipole moment vib and polarizability ␣ vib can be obtained for the static limit of these properties.The ␤ vib and ␥ vib are fifthand sixth-order properties, respectively, that could be exactly evaluated including the fifth-and sixth-order terms in the power series expansion of the potential energy surface.Thus using the fourth-order truncation for the potential energy ͓Eq.͑7͔͒ only partial values of ␤ vib and ␥ vib can be calculated.

A. Nuclear relaxation contributions
Recently, Castiglioni et al. 29,30 have pointed out that the nuclear relaxation contributions to electrical properties are due to the change of the equilibrium geometry induced by the applied field.Nuclear displacements of a molecule, caused by an electric field, can be easily obtained from normal coordinates displacements applying the stationary point condition to Eq. ͑7͒.
Iterative solution of the resulting system of nonindependent equations is detailed in Appendix A, and leads to the equilibrium field-dependent normal coordinates given by

͑8͒
where runs over the 3NϪ6 normal coordinates.
It is worth noting that this procedure allows for the direct determination of the field-dependence behavior of the nuclear coordinates.For instance, when finite field techniques 13,[19][20][21][22][24][25][26] are used molecular internal coordinates are reoptimized for each field strength. In thepresent study, only one calculation at the zero-field optimised geometry is required, followed by trivial application of Eq. ͑8͒.
Evaluating the first derivative of Eq. ͑8͒ with respect to electric field strength at zero field case, the first order nuclear relaxation q 1 ,a is obtained.
Moreover, the first-order nuclear relaxation has been previously defined for several authors as the ratio between the dipole moment derivative and the force constant, although only for diatomic molecules. 56It must be noted that both dipole derivatives and force constants can be obtained experimentally from the ir bands.[59] q 2 ,ab ϭ a 12 ,ab 2a 20 , q 3 ,abc ϭ a 13 ,abc 2a 20 .͑10͒ Substitution of the equilibrium field-dependent normal coordinates ͓Eq.͑8͔͒ into Eq.͑7͒ will lead to a field-dependent potential energy at the equilibrium geometry V eq (F), which will include both the electronic and nuclear relaxation contributions, given by

͑11͒
Comparison between this equation and the Taylor series ͓Eq.͑1͔͒ and subtraction of the purely electronic contributions (a 01 , a 02 , a 03 and a 04 terms͒ leads to definition of the nuclear relaxation contributions to molecular electrical properties for all components of their tensors.This definition of the P nr contributions to the static electrical properties shows the ad-ditive character between P el and P nr .At the equilibrium geometry and for the zero-field case, each component of the nuclear relaxation dipole moment has been shown to be null Each component of the nuclear relaxation polarizability is given by It can be seen that components of the ␣ nr tensor are a function of only harmonic terms ͑second-order terms͒.Each component of the nuclear relaxation first hyperpolarizability is given by ͑14͒ which are functions of harmonic ͑a 20 , a 11 , and a 12 ͒ and first order anharmonic ͑a 30 and a 21 ͒ terms ͑second and third order terms͒.This expression of ␤ abc nr cannot be simplified to because a 12 and q 1 in the first term, and a 21 and q 1 in the second are not equivalent.Therefore, all permutations of the field indexes abc must be considered to evaluate the ␤ nr values correctly.These permutations are not present in ␣ nr because a 11 and q 1 are equivalent in Eq. ͑13͒ ͓see Eq. ͑9͔͒.
Finally, each component of the nuclear relaxation contribution to the second hyperpolarizability is given by

3NϪ6
ͩ a 40 i jkl q 1 i,a q 1 j,b q 1 k,c q 1 l,d ϩ 3a 30 i jk a 21 kl,a a 20

͑16͒
where the permutation of all field indexes abcd has not been specified, like in the ␤ nr expression ͓Eq.͑14͔͒ for simplicity reasons, and is given in detail in Appendix B.

B. Vibrational contributions
To obtain the vibrational contribution to the molecular properties, derivatives of the vibrational energy with respect to the field strength at the corresponding order must be calculated.We assumed here that the vibrational energy is given by the harmonic zero-point energy ͑for considering excited vibrational states see Ref. 19͒ On the contrary, anharmonicity has been included through the potential energy expansion ͓Eq.͑7͔͒.The validity of this approximation for the polarizability of the HF molecule in the vibrational ground state has been already tested. 51,52nder these assumptions, the vibrational contribution to static electrical properties is function of the field-dependent force constants (k i (F)ϭ i ) and their derivatives with respect to the field strength.Then, in atomic units the tensor components of vib and ␣ vib are given by

͑18͒
where m i are the reduced masses associated to each normal coordinate, and i ϭ2 i , where i are the harmonic vibrational frequencies.Both the field-dependent force constants and their derivatives are evaluated at the zero-field equilibrium geometry.By double differentiating the power series expansion ͓Eq.͑7͔͒, the field-dependent force constants at the equilibrium geometry are given by where Q i eq are the equilibrium field-dependent normal coordinates presented in Eq. ͑8͒, and runs over all normal coordinates.From Eqs. ͑8͒ and ͑19͒ field-dependent force constant derivatives with respect to electric field are a function of Q i eq derivatives and the a 21 , a 22 , a 30 , a 31 , and a 40 coefficients of the power series expansion of the potential energy.All these coefficients are included in the potential energy expansion only when anharmonic corrections are considered.Then, in the harmonic model for the potential energy, P vib contributions to the molecular properties are zero everywhere.Therefore, vibrational contributions to the electrical properties are a consequence of the anharmonicity of the potential energy hypersurface.This aspect was originally pointed out by Kern and Matcha 44 in the late sixties.Thus final expressions of the vibrational contributions to dipole moment and polarizability tensor components are given by

͑21͒
where all permutations between the nonequivalent coefficients in each term have already been taken into account.
It has been shown that complete evaluation of P vib must include all a nm terms for which the following conditions hold ͑i͒ nϩmр1ϩ2, ͑ii͒ mр1, and ͑iii͒ n 0, where 1 equals 1 for vib , 2 for ␣ vib , and so on.For this reason, when the potential energy expansion includes all terms up to fourthorder ͑with all a nm with nϩmр4͒ like Eq. ͑7͒ only vib and ␣ vib can be obtained.Then, inclusion of higher-order terms in the field-dependent potential energy expansion ͓Eq.͑7͔͒ will not incorporate additional terms in either vib or ␣ vib , because i and their derivatives are evaluated for the zerofield case at the equilibrium geometry.
To obtain expressions for ␤ vib and ␥ vib , the same procedure must be followed including fifth-and sixth-order terms in the potential energy expansion ͓Eq.͑7͔͒, respectively.Otherwise, the calculated values of ␤ vib and ␥ vib will be only approximate.

III. METHODOLOGY
Calculations have been carried out at the ab initio SCF-MO level of theory.Electron correlation has been introduced using second-order Moller-Plesset perturbation theory ͑MP2͒.The Dunning-Huzinaga 60 basis set has been used in this work; for water and pyridine, the (9s5p/4s)/ ͓5s3p/3s͔ contraction has been used, and two sets of polarization functions have been added ͑TZ2P basis set͒ with exponents 0.75 and 0.15 for C, 0.80 and 0.15 for N, 0.85 and 0.15 for 0, and 0.75 and 0.15 for H.The C 2 axis and the molecular plane of these molecules have been assigned to the z and yz Cartesian coordinates of the molecular system of reference, respectively.
Nuclear contributions to electrical properties are given in terms of energy derivatives with respect to either normal coordinates, field strength, or both.To show the usefulness of the AEEP method, data presented in this work use only second-and third-energy derivatives that can be routinely calculated by the GAUSSIAN-94 61 series of programs.At both SCF and MP2, derivatives through second-order have been evaluated analytically, whereas third-order derivatives have been obtained by numerical differences.At the SCF level, P el has also been obtained analytically. 54,55Nuclear contributions ͑P nr and P vib ͒ to the static electrical properties have been calculated by the AEEP 62 code developed in our laboratory, following straight forward expressions presented in Sec.II.Because of only second-and third-energy derivatives having been calculated, the values of ␣ vib and ␤ vib presented in this paper might be improved by including higher-order energy derivatives.Furthermore, ␣ nr , ␤ nr , and vib values reported for water and pyridine can be improved only by including a more accurate evaluation of electron correlation or by using more flexible basis sets to determine the wave function of the molecular system.
The present procedure may be compared to the finite differences approaches 13,[19][20][21][22][24][25][26] that involve several geometry optimizations. In thepresent study, only one typical energyϩthird derivatives calculation at the field-free geometry equilibrium geometry is required, together with straight forward application of formulas of Sec.II.Thus the AEEP method is simple, easy to use and applicable to any molecule, with the practical bottleneck lain in the third-energy derivatives calculation.

IV. RESULTS
In this section the AEEP method is used to compute analytically nuclear contributions to static electrical properties of water and pyridine.The water molecule, which is one of the simplest nonlinear polyatomic molecules, has already been chosen by different authors 22,[44][45][46]49 to evaluate nuclear contributions to the static electrical properties using different theoretical approaches. Calulated AEEP electrical properties of water and pyridine will be compared with previous theoretical results, and available experimental data.For these two molecules, the agreement between the AEEP results and the experimental data will be shown to be excellent.
In Table I, the nonzero tensor elements of water obtained by the AEEP method are collected.These data have been computed at the respective minima on the SCF/TZ2P and MP2/TZ2P potential energy surfaces.From values presented in Table I, the total dipole moment (ϭ z ), the mean polarizability ͓͗␣͘ϭ(␣ xx ϩ␣ yy ϩ␣ zz )/3͔ and the vector component of the hyperpolarizability tensor parallel to the permanent dipole moment (␤ ʈ ϭ␤ z ϭ3(␤ xxz ϩ␤ yyz ϩ␤ zzz )/5) are 0.773, 8.51, and Ϫ3.47 a.u. at the SCF level, and 0.737, 9.72, and Ϫ9.59 a.u. at the MP2 level, respectively.
The experimental total dipole moment and mean polarizability of water are 0.724, and 9.92 a.u., 63 respectively.While the and ␣ SCF values show relative errors of 7% and Ϫ14%, respectively, these errors are reduced to 1.7% and Ϫ2.0% by inclusion of electron correlation at MP2.Such an efficiency of the MP2 level of theory to obtain accurate values of electrical properties has already been shown by different authors. 26,64,65For the water molecule, the agreement between the P N contributions calculated at the MP2 and SCF levels is excellent, hence electron correlation affects mainly to P el contributions.This fact must be study in more detail before being generalized.Thus a good strategy to reproduce experimental values might be to compute P N at SCF/TZ2P and only consider electron correlation to calculate P el . 26Using this strategy, , ͗␣͘, and ␤ ʈ turn out to be 0.738, 9.77, and Ϫ9.22 a.u., respectively.The relative errors with respect to the experimental dipole moment and polarizability are 1.9% and Ϫ1.5%, respectively.In consequence, the nuclear contributions to the electrical properties seem to be accurately calculated at the SCF/TZ2P level.
Nuclear contributions to electrical properties have also been calculated at the SCF/6-311ϩϩG(3d f ,2pd) level to check the reliability of the TZ2P basis set when predicting P N .Using this large and flexible basis set the vib , ͗␣ nr ͘, ͗␣ vib ͘ and ␤ z vib values are unchanged with respect to the TZ2P ones, and only ␤ z nr has changed to 7.44 a.u.As we will see in following paragraphs, the basis set effect on ␤ nr , which seems to be mainly related to the polarizability derivative, must be studied in detail.
Comparison of the P N AEEP values with earlier theoretical data must be done carefully in order to compare equivalent data, even when they are named in different ways.Both the finite field values [19][20][21]26 and the AEEP data of parallel components of the molecular properties ͑ z , ␣ zz , and ␤ zzz ͒ presented in this work can be compared directly, because they represent two different approaches to obtain the P nr and P vib contributions to electrical properties. Ken and co-workers [44][45][46] calculated the zero-point vibrational averaging contribution to the dipole moment ZPVA .This contribution was a function of the anharmonicity of the potential energy surface and dipole moment derivatives.Then, this calculated ZPVA can be compared with the vib value as given by Eq. ͑20͒.Werner and Meyer 47 used also Kern's expressions to compute ZPVA and ␣ ZPVA . Ten, their reported values can be compared with the AEEP vib and ␣ vib values, yet ␣ ZPVA only includes the two first terms of ␣ vib ͓Eq.͑21͔͒.The recent work by Russell and Spackman 48 includes terms up to fourth-order, so their ␣ ZPVA values can be compared with AEEP's ␣ vib .Moreover, Pandey and Santry 49 defined vibrational corrections to dipole moment v , polarizabilities ␣ v , and hyperpolarizabilities ␤ v .They went one step further by partially including the nuclear relaxation contribution, that was recognized to be 1 order of magnitude larger than the vibrational averaging values obtained by Kern et al. [44][45][46] Pandey and Santry 49 presented the harmonic part of P nr and the field derivatives of the force constants part of P vib together, the latter given by a 2m , and the former by the a 1m coefficients ͓Eq.͑20͔͒.These authors also obtained a zero value for nr . Martı´and Bishop 53 showed that the perturbative nuclear contributions P v ϩP ZPVA values can only be compared with the P nr ϩP vib values given either by the AEEP method or by the finite field treatment, even though the perturbation method expressions of P v ϩP ZPVA include only terms up to third-order ͑a mn with mϩnр3͒ together with first and second derivatives of the molecular property with respect to normal coordinates ͑a 1m and a 2m coefficients͒.The latter had already been included through P ZPVA , originally by Kern and co-workers.Moreover, Cohen et al. 22 calculated ␣ v , ␤ v , and ␥ v , from finite differences of the vibrational energy.This approach deserves especial attention because they introduced both mechanical and electrical anharmonicity in the potential energy expression, and also introduced anharmonicity in the vibrational energy expression ͓Eq.͑8͒ of Ref. 22͑a͔͒.Although they tried to use Bishop and Kirtman's notation arising from perturbation theory ͑P v for the pure vibrational contribution, and P ZPVA for the vibrational averaging͒, they actually computed the so-called nuclear relaxation and vibrational contributions to electrical properties ͑see Ref. 53͒.
As far as nuclear relaxation contribution is concerned, since both the Cohen method and our AEEP method introduce anharmonicity up to fourth-order in the potential energy, derivatives of V F 0 ͓Eq.͑8͒ of Ref. 22͑a͔͒ are completely comparable to AEEP P nr .For instance, their second derivative of V F 0 (␣ v ͓0200͔) is equivalent to our ␣ nr .As Cohen et al. show ͓and also the present Eq.͑13͔͒ the nuclear relaxation contribution to polarizability is not affected by anharmonicity at all.While the AEEP method uses the harmonic aproximation to calculate the vibrational energy, Cohen et al. introduced anharmonicity and the E 0 Dunham coefficient.Hence, derivatives of vibrational energy and the AEEP P vib values are only partially comparable.Thus their approach is indeed somewhat more accurate than ours, although it is more expensive and the amount of the computational burden is much larger, because geometries must be reoptimized a larger number of times, and many dipole-moment higher derivatives must be calculated at every strength field ͑e.g., for ␤ vib of a C 1 -symmetry molecule they need ten calculations involving four field strengths͒.On the contrary, the AEEP method needs only one calculation of energy derivatives at the field-free equilibrium geometry.
In Table II, the nuclear contributions to electrical properties of H 2 O reported by different authors are gathered and compared to the AEEP results.The ␣ zz nr , ␤ zzz nr , and z vib values obtained by the finite-field approach 26 and the AEEP method are nearly identical, because all terms of Eqs.͑13͒, ͑14͒, and ͑20͒ have been included in this work.Only the AEEP ␣ zz vib and ␤ zzz vib values are lower than the finite-field ones, because the AEEP data have been calculated here including only second and third order energy derivatives ͓Eq.͑21͔͒.These results are strongly encouraging, and show that the error due to the lack of electron correlation in P el is much more important than the error due to neglecting of fourthand fifth-order terms in ␣ vib and ␤ vib .By using a different contraction scheme of the Dunning-Huzinaga basis set, Kern et al. [44][45][46] obtained a vibrational contribution to larger than the AEEP one by a factor of 3. The P ZPVA values reported by both Werner and Meyer, 47 and Russell and Spackman, 48 for dipole moment and polarizability are in good agreement with the AEEP values.Both works calculated as well the diagonal tensor elements ␣ xx vib and ␣ yy vib to be 0.13 and 0.33 a.u., 47 and 0.0887 and 0.3982 a.u., 48 respectively.These values are also in excellent agreement with data reported in this work.Although Pandey and Santry 49 included partially the P nr contribution and used a different contraction scheme of the Dunning-Huzinaga basis set, agreement between their results and our AEEP values is satisfactory, especially for ␤.Data from the perturbative sum over states presented in the last column of Table II, are clearly lower than other theoretical predictions; further, they have been obtained using a completely different basis set.
One can see in Table II that, for the parallel component of the nuclear contribution to polarizability, the difference between results in Ref. 22 and the AEEP values is lower than 0.03 a.u.The origin of this small difference can be found in the second derivative of the vibrational energy, since no meaningful differences are found between ␣ zz v ͓0200͔ ͑0.876 a.u.͒ and ␣ zz nr ͑0.87 a.u.͒ values.On the contrary, for the first hyperpolarizability, ␤ v differs from ␤ nr by a factor of 2. The difference in sign in the nuclear contribution to ␤ is due to the different criteria used in the dipole moment definition.In this case, the main source of divergence must be found between the third derivative of V F 0 (Ϫ12.287a.u.) and ␤ zz nr ͑5.92 a.u.͒.Analyzing these two values in more detail, we find that this difference is caused by the first derivative of the polarizability with respect to normal coordinates, which is very sensitive to the basis set used in the calculation.This basis set effect on the computation of ␤ N , which has been previously shown in this work for the 6-311 ϩϩG(3d f ,2pd) basis, might be the origin of the difference among the three different results, namely ͑1͒ those found by Pandey et al. 49 and us, ͑2͒ those reported by Cohen et al. 22͑a͒ and ͑3͒ those found by Bishop and co-workers. 33s a test of the behavior of the AEEP method on medium-sized polyatomic molecules, we report calculations on pyridine.Table III shows the nonzero tensor elements of the electronic and nuclear contributions to the molecular properties of pyridine, predicted by the AEEP method and calculated at the SCF/TZ2P level, and the nonzero tensor elements of el , ␣ el and ␣ nr , computed at the MP2/TZ2P level.Comparison of the last two columns of Table III shows that the agreement between the calculated AEEP and ͗␣͘ and the experimental values 66 is excellent.Relative errors are lesser than 1% and Ϫ4% for SCF dipole moment and main polarizability, respectively.This fact supports the usefulness of the AEEP method to compute electrical properties for large polyatomic molecules.For pyridine, as for water, the SCF and MP2 ␣ nr values are very similar.However, when the SCF ␣ nr value is replaced by the MP2 result, the main polarizability value is closer to experiment.Finally, we have applied a strategy similar that used in water, and have computed the P el contributions by MP2 trying to reproduce the experimental results.These recomputed z , ͗␣͘ and ␤ l turn out to be 0.8926, 64.54, and 1.70 a.u.The relative error of mean polarizability is thus reduced to less than 1%, showing again that the main source of error in the SCF electrical properties arises from the P el contribution.The relative weights of the P N contributions in the total z , ͗␣͘ and ␤ l are 2%, 4.5% and more than 50%, respectively.These weights are of the same order of magnitude of the values previously reported for water.This result shows that nuclear contributions to the electrical properties are also essential to yield theoretical results that reproduce the experimental electrical properties of large polyatomic molecules like pyridine, and that theoretical values can be routinely compute by the AEEP method as presented in this work.

V. CONCLUSIONS
We have presented a method to analytically compute the nuclear contributions (P nr ϩP vib ) to electrical properties.Expressions to compute both nuclear relaxation and vibrational contributions have been deduced from a power series expansion of the potential energy.Such contributions to the electrical properties are given in terms of energy derivatives with respect to normal coordinates, field strength or both.The accuracy of the AEEP values is only determined by the quality of the wave function used to describe the molecular system.The AEEP method is quite simple, despite some amount of algebraic development to obtain equilibrium field-dependent normal coordinates in Appendix A. For end- users, the AEEP method requires just one calculation of the energy and its derivatives at the field-free optimized geometry, followed by trivial application of the formulas in Sec.II.
The AEEP method predicts the order of the necessary energy derivatives required for a complete computation of a specific nuclear contribution ͑e.g., for nuclear relaxation to ␣ only second derivatives are necessary͒.To our knowledge, this interesting advantage is exclusive of AEEP method and allows important savings in computational time.
For water and pyridine, the reported SCF nuclear contributions combined with MP2 electronic contributions, allow to reproduce well the experimental data, showing error lesser than 2%.The calculated nuclear contributions of water and pyridine represent a meaningful amount of these molecular properties.The fourth-order terms included by Cohen et al., and by Russell and Spackman allow to quantitatively obtain the ␣ vib values; on the contrary, when these terms are neglected, like in the presented AEEP values, only qualitatively ␣ vib values are obtained.The analytical evaluation of electrical properties is the only systematic and computational feasible method for dealing with large polyatomic molecules using standard quantum mechanical programs.The presented ␣ vib and ␤ vib AEEP values of electrical properties might be improved by inclusion of fourth-order derivatives.Further work in this subject is being carried out in our laboratory.

ACKNOWLEDGMENTS
This work has been funded through the Spanish DGI-CYT Project No. PB92-0333.One of us ͑J.M.L.͒ acknowledges the financial help provided by the Generalitat de Catalunya through the CIRIT Project No. FI/95-5101.We also thank Mr. Josep Martı ´for helpful discussions.

APPENDIX A
Application of stationary conditions to Eq. ͑7͒ allows to obtain the nuclear displacements of a molecule induced by the applied field, which are given by ‫ץ‬V͑Q,F͒ ‫ץ‬Q ϭ2a 20 where runs over the 3NϪ6 normal coordinates.This nuclear displacement can be calculated from the first term, x,y,z a 13 ,abc F a F b F cͬ , ͑A2͒ where each normal coordinate Q depends on all normal coordinates and the field strength vector.Thus to obtain the field-dependent normal coordinates Q , a nonlinear system of 3NϪ6 equations with 3NϪ6 nonindependent variables must be solved.This system of equations can be solved iteratively by the following procedure: ͑a͒ Use the set of independent field-dependent normal coordinates given by that is obtained by applying the stationary condition to the harmonic expansion of the potential energy.͑b͒ Substitute this set of independent normal coordinates Q H into Eq.͑A2͒ to get a first set of field-dependent normal coordinates Q Ј .͑c͒ Substitute the new set of normal coordinates Q Ј into Eq.͑A2͒ to obtain a second set of field-dependent normal coordinates Q Љ .
͑d͒ Repeat step ͑c͒ until the desired convergence is achieved.
In step ͑b͒ the converged first-order field-dependent normal coordinates given by Q eq ͑ F x ,F y ,F z ͒ϭϪ ͚ a x,y,z q 1 ,a F a ͑A4͒ are obtained, whereas in step ͑c͒ the converged second-order field-dependent normal coordinates are obtained, x,y,z ͫ q 2 ,ab Ϫ ͚ i 3NϪ6 a 21 i,a a 20 q 1 i,b ϩ ͚ i, j 3NϪ6 3a 30 i j 2a 20 q 1 i,a q 1 j,b ͬ F a F b ͑A5͒ and so on.The desired convergence is limited by the expansion of the potential energy used.In the present case, because all terms up to fourth-order have been included in the expansion of the potential energy, the field-dependent normal coordinates converged can be obtained only up to thirdorder.However, to reach the reported expressions of the nuclear contributions to the electrical properties only the field-dependent normal coordinates up to third order are needed.

APPENDIX B
The complete expression of the nuclear relaxation contribution to the second hyperpolarizability, with all explicit permutations is ␥ abcd nr ϭ6 ͚ i 3NϪ6 ͑ a 13 i,abc q 1 i,d ϩa 13 i,abd q 1 i,c ϩa 13 i,acd q 1 i,b ϩa 13 i,bcd q 1 i,a ͒ϩ4 ͚ i 3NϪ6 ͑a 12 i,ab q 2 i,cd ϩa 12 i,ac q 2 i,bd ϩa 12 i,ad q 2 i,bc ͒ Ϫ4 ͚ i, j 3NϪ6 ͑a 22 i j,ab q 1 i,c q 1 j,d ϩa 22 i j,ac q 1 i,b q 1 j,d ϩa 22 i j,ad q 1 i,b q 1 j,c ϩa 22 i j,bc q 1 i,a q 1 j,d ϩa 22 i j,bd q 1 i,a q 1 j,c ϩa 22 i j,cd q 1 i,a q 1 j,b ͒ Ϫ4 ͚ i, j 3NϪ6 ͑a 21 i j,a q 1 i,b q 2 j,cd ϩa 21 i j,b q 1 i,a q 2 j,cd ϩa 21 i j,a q 1 i,c q 2 j,bd ϩa 21 i j,c q 1 i,a q 2 j,bd ϩa 21 i j,a q 1 i,d q 2 j,bc ϩa 21 i j,d q 1 i,a q 2 j,bc ϩa 21 i j,b q 1 i,c q 2 j,ad ϩa 21 i j,c q 1 i,b q 2 j,ad ϩa 21 i j,b q 1 i,d q 2 j,ac ϩa 21 i j,d q 1 i,b q 2 j,ac ϩa 21 i j,c q 1 i,d q 2 j,ab ϩa 21 i j,d q 1 i,c q 2 j,ab ͒ ϩ6 ͚ i, j,k 3NϪ6 ͑a 31 i jk,a q 1 i,b q 1 j,c q 1 k,d ϩa 31 i jk,b q 1 i,a q 1 j,c q 1 k,d ϩa 31 i jk,c q 1 i,a q 1 j,b q 1 k,d ϩa 31 i jk,d q 1 i,a q 1 i,b q 1 k,c ͒ ϩ12 ͚ i, j,k 3NϪ6 ͑ a 30 i jk q 1 i,a q 1 j,b q 2 k,cd ϩa 30 i jk q 1 i,a q 1 j,c q 2 k,bd ϩa 30 i jk q 1 i,a q 1 j,d q 2 k,bc ϩ3 more terms ͑ a 40 i jkl q 1 i,a q 1 j,b q 1 k,c q 1 l,d ͒Ϫ6 ͚ i, j,k,l 3NϪ6 ͩ a 30 i jk a 21 kl,a a 20 k q 1 i,c q 1 j,d q 1 l,b ϩ a 30 i jk a 21 kl,b a 20 k q 1 i,c q 1 j,d q 1 l,a ϩ a 30 i jk a 21 kl,a a 20 k q 1 i,b q 1 j,d q 1 l,c ϩ9 more terms••• ͪ ϩ9 ͚ i, j,k,l,m 3NϪ6 ͩ a 30 i jk a 30 klm a 20 k q 1 i,a q 1 j,b q 1 l,c q 1 m,d ϩ a 30 i jk a 30 klm a 20 k q 1 i,a q 1 j,b q 1 l,c q 1 m,d ϩ4 more terms••• ͪ .͑B1͒

3 ,
abc F a F b F cͬ ͑A3͒
a See Reference 63.
d g h k

TABLE III .
Nonzero tensor elements of the electronic and nuclear contributions to the electrical properties of pyridine calculated at SCF/TZ2P and MP2/TZ2P.All values are given in atomic units ͑a.u.͒.