Finite field treatment of vibrational polarizabilities and hyperpolarizabilities : On the role of the Eckart conditions , their implementation , and their use in characterizing key vibrations

In the finite field ~FF! treatment of vibrational polarizabilities and hyperpolarizabilities, the field-free Eckart conditions must be enforced in order to prevent molecular reorientation during geometry optimization. These conditions are implemented for the first time. Our procedure facilities identification of field-induced internal coordinates that make the major contribution to the vibrational properties. Using only two of these coordinates, quantitative accuracy for nuclear relaxation polarizabilities and hyperpolarizabilities is achieved in p-conjugated systems. From these two coordinates a single most efficient natural conjugation coordinate ~NCC! can be extracted. The limitations of this one coordinate approach are discussed. It is shown that the Eckart conditions can lead to an isotope effect that is comparable to the isotope effect on zero-point vibrational averaging, but with a different mass-dependence. © 1999 American Institute of Physics. @S0021-9606 ~99!01225-8#


I. INTRODUCTION
Since the laser was invented in 1960, interest in nonlinear optics ͑NLO͒ has grown continuously 1-3 and now ranges from fundamental studies of the interactions of light with matter to applications such as laser frequency conversion and optical switching.In the presence of ͑spatially uniform͒ static and dynamic electric fields, both the electronic and nuclear motions of a molecule will be perturbed.][10][11] Until very recently, 12 the only method available for calculating the vibrational contribution at an arbitrary optical frequency was the perturbation treatment of Bishop and Kirtman ͑BK͒, [13][14][15][16] which is based on the general sum-over-states formulas 17 for dynamic hyperpolarizabilities.In the BK procedure, the electronic and vibrational terms are separated by applying a canonical or clamped nucleus approximation, wherein the electronic and nuclear motions are considered sequentially, rather than simultaneously. 18The resulting expressions, given in order of electrical and mechanical anharmonicity, contain first-and higher-order derivatives of static clamped nucleus electrical properties with respect to normal coordinates, as well as harmonic and anharmonic vibrational force constants.The determination of these parameters, particularly the higher property derivatives and anharmonic force constants, can be computationally demanding.
0][21][22][23] In the FF approach, one first optimizes the molecular geometry in the presence of a static electric field and, then, calculates the electronic dipole moment ͑or energy͒, linear polarizability, and first hyperpolarizability in that field.A detailed analysis shows that the field-dependent change in the clamped nucleus electronic properties yields what is called 7 the nuclear relaxation ͑NR͒ contribution to the vibrational polarizability and hyperpolarizabilities.The electronic dipole moment ( e ) gives the static value of the nuclear relaxation term, whereas the values derived from the linear polarizability (␣ e ) and first hyperpolarizability (␤ e ) correspond to the infinite optical frequency limit 19 ͑also known as the enhanced approximation 24 ͒.This limit, it turns out, gives results similar to the exact frequency-dependent value at typical laser optical frequencies, at least in those cases that have been tested. 25here is a second source of vibrational polarizability and hyperpolarizabilities that arises because of the static fieldinduced change in the vibrational potential.It is known as the curvature contribution and can be obtained from the zeropoint vibrational averaging ͑ZPVA͒ correction due to motion about the field-dependent equilibrium geometry.The analysis 19 of the field-dependent ZPVA is exactly analogous to that of the clamped nucleus electronic property itself.As might be expected, the contributions to the ZPVA that arise from geometry relaxation are of higher order than the NR contributions associated with the clamped nucleus electronic properties.
Very recently the FF scheme has been generalized 12 to yield the NR and curvature contributions over the entire frequency range.This formulation, like its predecessors, begins with a geometry optimization in the presence of a static electric field.The difference occurs in a later step where one must evaluate dynamic, rather than static, properties at the field-equilibrated geometry.Calculations based on this new technique have yet to be carried out.Indeed, even the earlier and simpler FF approach has been applied only in a few special cases.The main reason is that, in general, the molecule will rotate during the field-dependent geometry optimization.This causes the applied field to be no longer in the intended direction with respect to the molecular axes.The sole exception 7,18,26,27 occurs when the permanent ͑if any͒ and field-induced dipole moments are in the same direction as the field.
It is known in principle 19 that the above difficulty can be overcome by proper application of the Eckart conditions. 28,29o our knowledge, however, this has not previously been accomplished.In Sec.II we present a simple, yet successful, implementation of the Eckart conditions for geometry optimization in the presence of a static electric field.
Our computer code permits us to determine the contribution of any coordinate, or set of coordinates, to the vibrational polarizability and hyperpolarizability.In Sec.III we present a systematic method for determining the key internal vibrational coordinate͑s͒ and apply this approach to three prototypical -conjugated molecules.Although it has been widely speculated that the vibrational contribution to the electrical properties of such molecules is dominated by a single internal coordinate-sometimes referred to as the effective conjugation coordinate ͑ECC͒ 30 -there are also indications 31,32 to the contrary.Our treatment of the three -conjugated molecules shows that, in fact, two internal coordinates are required to quantitatively reproduce the electrical properties.Furthermore, we are able to extract a single most efficient natural conjugation coordinate ͑NCC͒ from these calculations.The NCC has its limitations but is shown to be preferable to the ECC.
The Eckart conditions explicitly involve atomic masses and thereby give rise to orientational mass-dependent effects on the vibrational polarizability and hyperpolarizability that have not been generally recognized in either the BK or FF treatments.We analyze these effects in Sec.IV using illustrative calculations for three isotopes of water (H 2 O, D 2 O, and HDO͒.Finally, Sec.V contains a summary of this paper, as well as some of our plans for future investigations based on the FF approach with the Eckart conditions included.

II. IMPLEMENTATION OF ECKART CONDITIONS
As discussed in Sec.I, the first step in calculating vibrational contributions to polarizabilities and hyperpolarizabilities by the FF method lies in determining field-dependent equilibrium geometries.This has limited the treatment since, except in a few special cases, the molecule will rotate with respect to the finite field during geometry optimization.In order to prevent the molecule from rotating, one must apply the Eckart conditions, 28,29 ⌺ i m i ͑ a i ϫr i ͒ϭ0.͑1͒ Here, r i is the vector from the center of mass to atom i, a i is the equilibrium value of that vector, and m i is the atomic mass.͑For convenience the translational condition is not explicitly displayed since it is not involved in the following analysis.However, it is rigorously enforced in our treatment.͒A key point to recognize is that, according to the theory underlying the FF method, it is the field-free Eckart conditions that are required even when an external field is imposed.Thus, the field-dependent geometry optimization may be carried out by the following procedure: ͑i͒ Using any standard quantum chemistry program a complete optimization of the field-free geometry is performed; ͑ii͒ At this geometry, three rotation ͑two for linear mol-ecules͒ and three translational Eckart coordinates are constructed; ͑iii͒ A set of variable coordinates is read as input.These coordinates may be of any type ͑internal, normal etc.͒ and any number from one to 3N-6(5); ͑iv͒ The input coordinates expressed in terms of massweighted Cartesian displacements are orthogonalized to the field-free translations and rotations; and ͑v͒ A field-dependent geometry optimization is performed in this subspace of orthogonalized coordinates.
In our implementation, step ͑v͒ is carried out by means of the Broyden, Fletcher, Goldfarb, and Shanno ͑BFGS͒ method. 33he required field-dependent energy and gradients can be obtained from any standard quantum chemistry program.We developed our code in conjunction with GAUSSIAN 94. 34 Our program was tested by determining fully optimized ͑3N-6 or 3N-5 coordinates͒ field-dependent geometries, and then calculating NR contributions to the vibrational polarizability and hyperpolarizabilities by means of the FF procedure of Bishop, Hasan, and Kirtman ͑BHK͒. 19As noted above, exactly the same results for the vibrational properties can be obtained analytically.However, the latter computations are much more demanding because, in contrast with the FF approach, they require explicit input values of: ͑1͒ quadratic, cubic, and quartic vibrational force constants; ͑2͒ first derivatives ͑with respect to nuclear coordinates͒ of the electronic dipole moment ( e ), polarizability (␣ e ), and first hyperpolarizability (␤ e ); ͑3͒ second derivatives of e and ␣ e ; and ͑4͒ third derivatives of e .The simplicity of the FF/ relaxation procedure, however, does not come without a price.In order to ͑implicitly͒ calculate the first, second, and third derivatives of ␤ e , ␣ e , and e , respectively, with sufficient accuracy, one needs to use tight thresholds in the field-dependent geometry optimization.Experience has taught us that the residual atomic forces must be no greater than 10 Ϫ6 hartree/bohr or hartree/rad and that the thresholds on the energy (10 Ϫ12 a.u.) and density matrix (max ϭ10 Ϫ10 a.u.; rmsϭ10 Ϫ12 a.u.) must be decreased accordingly.Nonetheless, except for small molecules, the FF method is far more efficient than the perturbation theory procedure for determining the complete NR contribution to the vibrational properties, particularly if the properties are longitudinally dominant.

III. CHARACTERIZATION OF MOST IMPORTANT VIBRATIONAL MOTIONS/FIELD-INDUCED COORDINATES
The results obtained from a FF calculation with any complete set of coordinates can also be determined by perturbation theory.In the perturbation treatment of ␣ NR , ␤ NR , and ␥ NR one computes these quantities as a sum over normal modes ͑SOM͒.At the double harmonic level of approximation there is no coupling between modes, and the contribution of each vibration can be immediately determined.When anharmonicity is taken into account, the assignment to individual modes can be accomplished 35 by equipartitioning.This provides a straightforward method of characterizing the vibrational coordinates that make the largest contributions to the NR polarizability and hyperpolarizabilities.There are, however, two important drawbacks to the SOM approach.One of these is the difficulty of doing a complete perturbation calculation for medium and large systems; in fact, no such calculations have been reported with anharmonicity included.Equally significant is the nonintuitive and nontransferable nature of the normal coordinates.Both of these drawbacks can be overcome through use of the FF method.
It is often assumed for -conjugated systems-which are of special interest in NLO applications-that there is a single coordinate, 30,36 known as the effective conjugation coordinate ͑ECC͒ or bond length alternation ͑BLA͒ coordinate, which makes a dominant contribution to the NR properties.To our knowledge, however, this hypothesis has not been assessed and, moreover, the ECC is usually not precisely defined.On the other hand, the FF method leads logically to field-induced coordinates ͑FICs͒, from which a most efficient natural conjugation coordinate ͑NCC͒ can be rigorously extracted and tested.In our first application to three prototype -conjugated systems ͑see Fig. 1͒-hexatriene, NO 2 Ϫ͑CHϭCH͒ 2 ϪNH 2 , and p-nitroaniline ͑PNA͒-we find that the NCC gives semiquantitative results for the properties of both donor/acceptor ͑D/A͒ molecules but, beyond that, a second FIC is required.Using two FICs in the subsequent restricted field-dependent geometry optimizations yields almost exact agreement in all cases with the values calculated using a complete set of 3N-6 coordinates.
The FICs are generated simply by taking the difference between the optimized geometry with and without an external field ͑F͒ being present.For the -conjugated systems in this study, F is chosen to lie in the longitudinal direction ͓parallel to the dipole moment for PNA and the donoracceptor ͑D/A͒ butadiene; along the line that bisects the central and end double bonds for hexatriene͔ since ␣, ␤, and ␥ are all longitudinally dominant.The special significance of this type of coordinate is suggested by physical as well as mathematical considerations.From a mathematical perspective the zero field and field-dependent geometries determine a linear synchronous transit ͑LST͒ path 37 between the two structures.By moving along this path one can obtain a first approximation to the optimum geometry at any field from zero to ͉F͉ and even somewhat larger.For PNA and NO 2 Ϫ͑CHϭCH͒ 2 ϪNH 2 , specifically, there is also an argument which suggests that this path has special physical significance.It is based on the assumption that the polarizability and hyperpolarizabilities arise from electric dipole transitions between a primarily covalent ground state and one ͑or more͒ charge-transfer excited states.9][40] Although this valence-bond charge-transfer model has serious deficiencies, 10,41 as currently formulated the underlying concept is basically valid.The choice that remains in defining a particular FIC is the magnitude and sign of the electric field.It seems clear that the field should be in the direction that induces the charge-transfer structure͑s͒, which is also the alternative that leads to the larger geometry change.Based on both the mathematical and physical perspectives, one is tempted to select the strongest field used in the FF calculation; this yields the strongest field-induced coordinate ͑SFIC͒.Our results for the NR polarizabilities and hyperpolarizabilities shown in Table I were obtained by restricted one coordinate field-dependent geometry optimizations.They may be contrasted with the reference values given in Table II, which were determined using a complete set of 3N-6 coordinates.
Typically, in FF calculations we employ fields of magnitude 2 k F where kϭ0,1,2,... and Fϭ0.0004 a.u.The factor 2 k is introduced in order to apply the Romberg fitting procedure 42,43 which removes higher-order hyperpolarizability contaminants.In this case a maximum value of kϭ5 was employed; thus, the strongest field is 0.0128 a.u.
The geometry optimizations were done at the RHF/6-31G level and the electronic properties needed for the FF treatment were calculated by the CPHF/6-31G method.For the two D/A molecules the one vibrational degree of freedom SFIC model yields small errors ͑see Table I for the % error in parentheses͒ in the static ␣ NR (0;0), the static ␥ NR (0;0,0,0), and the optical kerr effect ͑OKE͒ ␥ NR (Ϫ;,0,0).The error in the remaining properties varies from about 8% to 16%.In the case of hexatriene, the deviations from the complete 3N-6 coordinate treatment are more serious.This is consistent with the fact that, from the physical point of view, the SFIC model is less appropriate for hexatriene than for the two D/A molecules.
In order to improve on the SFIC model it was decided to add a second coordinate derived from a weak field geometry optimization.This weak field induced coordinate ͑WFIC͒ was obtained in the same manner as the SFIC except that we used the weakest field ͑0.0004 a.u.͒ from our FF calculations.The argument here is simply that the two limiting cases might be expected to provide the necessary flexibility over the entire range of fields.Indeed, in a single coordinate treatment, the WFIC gives results ͑see Table III͒ that are to a large extent, complementary to those found for the SFIC.In particular, for the D/A molecules the static ␤ NR (0;0,0); the dc-Pockels ͑dc-P͒ effect ␤ NR (Ϫ;,0), and the dc-induced second harmonic generation ͑dc-SHG͒ ␥ NR (Ϫ2;,,0) are given quite well, whereas the opposite is true in Table I.It turns out that the two-coordinate WFIC/SFIC model is remarkably accurate for all three molecules, with a maximum error of 1.6% in every instance.Although this is quite satisfactory, one final improvement of the two-coordinate model was undertaken.
Since there is some ambiguity with respect to the choice of the weakest and strongest fields, the sensitivity of our results was tested by examining the additional pairs 0.0008/ 0.0064 and 0.0016/0.0032a.u.Table IV summarizes the values obtained for the latter.Again, the agreement with the complete 3N-6 coordinate calculation is excellent.It is even better than for the WFIC/SFIC model and the same is true for the 0.0008/0.0064a.u.pair.In order to avoid any ambiguity, then, we suggest that the 2-FIC treatment be employed with coordinates generated by the fields F ϭ0.0016/0.0032a.u., which will typically be in the middle of the range used in the FF procedure.
Our results imply that, at each field used in the FF calculations, one can find a linear combination of two midrange FICs (Fϭ0.0016/0.0032a.u.) that accurately reproduces the fully optimized geometry.This linear combination varies, of course, from one FF to another.A simple average over all fields yields an overall single most effective coordinate, which we call the natural conjugation coordinate ͑NCC͒.When the NCC is used in a one coordinate FF treatment, the polarizabilities and hyperpolarizabilities shown in Table V are obtained.The results are quite reasonable, with the most notable exceptions being ␥ NR (0;0,0,0) and ␥ NR (Ϫ;,0,0) for hexatriene; in these cases the FF procedure is nonconvergent.The convergence problem can be overcome by symmetrizing the hexatriene NCC on the grounds that ϩF and ϪF are degenerate for that molecule.This leads to the rows contained in square brackets in Table V.Unfortunately, when that is done the values of the polarizability and dc-SHG are seriously degraded.One other possibility considered for the NCC was to use a weighted average in which F 2 is the weighting factor.Although the errors in ␥ NR (0;0,0,0) and ␥ NR (;,0,0) for the D/A molecules are, thereby, significantly reduced, the errors in all the other properties increase.Thus, we define the NCC as being the simple average and at the same time note the limitations pointed out above.
Finally, it is of interest to compare the rigorously computed NCC with the a priori ECC.As observed earlier, there is no strict definition of the ECC.Even for hexatriene, one can vary the CC bond lengths while keeping all angles fixed or move the CH groups (CH 2 at the chain ends͒ as a rigid unit in the longitudinal direction.We make the latter choice, depicted in Fig. 1͑a͒, since it turns out to give better results.
For the two D/A molecules the ECC in Figs.1͑b͒ and 1͑c͒ was drawn by alternating the longitudinal displacements as one proceeds from one atom to the next along the conjugated backbone ͑which includes the oxygen atoms of the NO 2 group͒.Again, the hydrogens move rigidly along with the heavy atom to which they are attached.When these coordinates are used in a one-dimensional FF treatment, the NR polarizabilities and hyperpolarizabilities reported in Unlike the ECC, the NCC displayed in Fig. 2 has a transverse as well as a longitudinal component ͑both shown separately using the same scale͒.That is why, for example, the latter gives a good value for the linear polarizability of hexatriene ͓cf.Ref. 35͔.The longitudinal component of the NCC is similar to the ECC in that both exhibit alternating displacements along the conjugated backbone.On the other hand, in the NCC the displacements are far from uniform in length and the attached hydrogen atoms are, for the most part, either stationary or nearly stationary.We have not attempted a further analysis of the NCC since that is beyond the scope of this paper, nor have we carried out any onedimensional calculations using only the longitudinal component of the NCC.

IV. ANALYSIS OF ISOTOPE EFFECT FOR NR VIBRATIONAL POLARIZABILITIES AND HYPERPOLARIZABILITIES
Although the perturbation theory calculation of ␣ NR , ␤ NR , and ␥ NR is equivalent to a FF treatment, the Eckart conditions appear in a more subtle way since the field does not occur explicitly and infinitesimal displacements are presumed.For example, consider the perturbation theory formula for the static ␣ NR given by 44 where R i is a ͑nonredundant͒ internal coordinate, F is the force constant matrix in terms of the R i , and , are Cartesian directions.At first glance it might appear that the expression on the rhs of Eq. ͑2͒ is independent of atomic masses.However, one should recall that the dipole partial derivatives must be evaluated with the coordinates corresponding to infinitesimal rotations ͑and translations͒ held equal to zero, i.e., the Eckart conditions ͓cf.Eq. ͑1͔͒ must be satisfied.As a result, although it is often not recognized, ␣ NR ͑0;0͒ will be mass-dependent; in particular, it will exhibit an isotope effect.In this section we analyze some of the factors that govern this mass-dependence not only for ␣ NR but for  the first hyperpolarizability (␤ NR ) and second hyperpolarizability (␥ NR ) as well.We also carry out illustrative calculations for three isotopes of water.
The expression on the rhs of Eq. ͑2͒ has been studied by vibrational spectroscopists [45][46][47] in connection with infrared isotope intensity rules.It has been shown that under certain circumstances some, or all, of the dipole first derivatives will be mass-independent.Three such instances, basically determined by symmetry, are: ͑1͒ R i is a symmetry coordinate that cannot give rise to an infinitesimal molecular rotation; ͑2͒ the molecule does not have a permanent dipole moment ͑ordinarily due to inversion symmetry͒; ͑3͒ the dipole derivative is in the direction of the permanent moment.Apart from these cases, the remaining dipole derivatives will generally be mass-dependent.
Next, we extend the above analysis to hyperpolarizabilities treated at the double harmonic level.This level is sufficient to yield the complete NR contribution to several NLO processes in the infinite frequency approximation 5,19,21 -namely ␤ NR (Ϫ;,0), ␥ NR (Ϫ2;,,0), and the intensity-dependent refractive index ͑IDRI͒ ␥ NR (Ϫ;,,Ϫ).For these properties the derivatives of ␣ e and ␤ e are needed in addition to those of e .Ordinarily, the derivatives of ␣ e and ␤ e will also be mass-dependent.
For ␣ e only exception ͑1͒ above applies; for ␤ e exception ͑2͒ is also valid, since inversion symmetry causes both e and ␤ e to vanish.There is no simple analog of ͑3͒ in the case of ␣ e and ␤ e since they are tensor, rather than vector, properties.
In order to obtain the NR contribution to other vibrational hyperpolarizabilities that are often of interest, including the static ␤ NR and ␥ NR as well as the OKE, higher-order electronic property derivatives are necessary.As a general rule, these derivatives will also be mass-dependent due to the Eckart conditions unless the symmetry coordinate cannot lead to a molecular rotation.
Thus far we have focused on NR to the exclusion of contributions to ␣ v , ␤ v and ␥ v that are related to the ZPVA of the electronic properties.Besides the ZPVA per se, the latter category includes higher-order ''curvature'' terms derivable from the ZPVA. 20͑The ZPVA is the lowest-order curvature term and is usually considered to be an electronic property but here, for convenience, we discuss it as part of the vibrational contribution.͒In contrast with NR, the isotope effect due to the ZPVA has been fairly well studied, [48][49][50] at least for the polarizability of small molecules.A massdependence occurs in this case not only through the electronic property derivatives but, more importantly, through averaging over the vibrational coordinates.For the total isotope shift both NR and ZPVA-related contributions must be taken into account.In this paper our interest lies particularly in the previously ignored NR isotope shift and its analysis in terms of the Eckart conditions.However, in the example ͑below͒ we also compute the ZPVA terms in order to compare them with NR.
For illustrative, and further analytical purposes, calculations were undertaken on the three ͑common͒ isotopes of water-H 2 O, D 2 O, and HDO.Our results turn out to be revealing as to the role played by electronic vs inertial symmetry in determining ␣ NR , ␤ NR , and ␥ NR .They also shed light on the relationship between the mass-dependence due to NR and that due to ZPVA.The H 2 O and D 2 O isotopes have just one internal coordinate ͑asymmetric stretch͒ that has the same symmetry as a rotation ͑which, in this case, is within the molecular plane͒.From the lhs of Eq. ͑1͒ one can determine 29 the mass-dependent angle of rotation generated by any set of displacement having that symmetry.In order to satisfy the Eckart conditions the molecule must be rotated to restore the original orientation, and that leads to a rotation of e ͑as well as the ␣ e and ␤ e tensors͒.Therefore, the electronic properties with respect to space-fixed axes will be isotope-dependent.HDO has the same C 2V electronic sym- metry as the other two isotopes ͑the electronic properties of all three are identical͒ but its inertial symmetry is reduced to C s .As a result, the two symmetric internal coordinates can also give rise to a mass-dependent rotation within the molecular plane.
Our NR results for H 2 O, D 2 O, and HDO are reported in Table VII.All calculations were carried out in the polarized ͓5s3 p2d/3s2 p͔ Sadlej basis 51 using six d functions.The permanent dipole moment lies along the z axis and yz is the molecular plane. 52Completely optimized field-dependent geometries were determined at the Hartree-Fock level using the Eckart program.Then the electrical properties were evaluated at the optimum geometry by the coupled perturbed Hartree-Fock ͑CPHF͒ method.For H 2 O our values agree very closely with those calculated by Bishop and Dalskov 25 using the same basis set; for D 2 O and HDO there are no previous results.
The largest components are those in the yz plane.Consider, first, the diagonal planar components of the static polarizability given by Eq. ͑2͒.If the magnitude of e is , then rotation in the molecular plane by an angle will generate the dipole components y e ϭ sin and z e ϭ cos .Since is proportional to the internal coordinate displacement, ␦R, the power series expansion of sin will give a leading term proportional to ␦R and a remainder of order (␦R) 3 .Similarly, z has the form ͓1ϩ(␦R) 2 ͔.Thus, ‫ץ‬ y e /‫ץ‬R 0 and ␣ yy NR (0;0) exhibit an isotope effect, while ‫ץ‬ z e /‫ץ‬R vanishes and there is no mass-dependence for ␣ zz NR (0;0).Analogous arguments can be made to show that ‫␣ץ‬ zz e /‫ץ‬R and ‫␤ץ‬ zzz e /‫ץ‬R are zero, which accounts for the iso-topic invariance of the dc-P effect and dc-SHG.Note that the above analysis applies to HDO even though not required by the inertial C s symmetry.
It is evident from Table VII that the reduced inertial symmetry in HDO leads to an isotope effect ͑vs H 2 O and D 2 O) for the diagonal z component when anharmonic terms are included ͓see ␤ NR (0;0,0), ␥ NR (0;0,0,0), and ␥ NR (Ϫ,,0,0)].For HDO, unlike H 2 O and D 2 O, symmetric motions can generate a rotation and, indeed, that is the origin of the mass-dependent anharmonicity contributions in this isotope.Thus, we see that the effect upon ␣ NR , ␤ NR , and ␥ NR of applying the Eckart conditions depends upon both the electronic and inertial asymmetry.In general, even in the absence of a significant electronic effect, substitutions that cause a rotation of the inertial axes will alter these properties.Finally, the above arguments can easily be extended to explain all the other isotope variations in Table VII.
In order to compare with the ZPVA contributions and, ultimately, with experiment, the various components given in Table VII were combined to form the mean properties ͑indicated by an overhead bar͒, the polarizability anisotropy ͑⌬␣͒, and the Kerr anisotropies ␤ K and ␥ K .For the polarizability we use the definitions of Ref. 50; for the hyperpolarizabilities, Ref. 5. Our results for the separate NR and ZPVA contributions are reported in Table VIII.For the ZPVA of H 2 O and D 2 O these values agree well with previous work. 49,53,54It is clear that the NR term behaves very differently from the corresponding ZPVA term as far as isotopic substitution is concerned.For NR, the HDO value is either the largest or the smallest of all three isotopes.On the other hand, the ZPVA  VIII, the range of isotope shifts in ␣ ͑0; 0͒ and ␤ ͑0; 0, 0͒ will be somewhat larger ͑by a factor of 1.5-2.5͒for ZPVA than for NR.This is just a preliminary assessment because electron correlation has not been taken into account and the effects of the higher-order curvature ͑ZPVA-related͒ terms remain to be examined.Our main purpose has been to determine the relative differences in magnitude and isotopic pattern between the ZPVA and NR vibrational contributions in a typical case.Recently, 55 gas-phase EFISH measurements have been reported which yield the ␤ ¯(Ϫ2;,) and ␥ ¯(Ϫ2;,,0) values for H 2 O and D 2 O.The NR contribution to both of these quantities is very small at optical frequencies.In order to emphasize the NR isotope effect it would be preferable to have measurements of the static properties, particularly ␣ and ␤, since it is difficult to calculate ␥ ZPVA .From Table VIII it is clear that, in these instances, the total isotope shift will be quite different from that predicted for ZPVA alone.The theoretical NR contribution to ␤͑0; 0, 0͒ can also be checked, at least approximately, by means of infrared and Raman intensity measurements. 56

V. CONCLUSIONS AND FUTURE WORK
In the FF treatment of vibrational polarizabilities and hyperpolarizabilities, a key step is determining the optimized molecular geometry in the presence of an external field.The major stumbling block in such a calculation is that the molecule will generally change its orientation with respect to the field during optimization.In order to avoid this difficulty we have developed a convenient implementation of the field-free Eckart conditions on rotational motion.These conditions introduce a mass-dependence into the derivatives of e ,␣ e ,␤ e with respect to internal coordinates.One consequence is an isotope effect for the NR polarizabilities and hyperpolarizabilities.Using the three common isotopes of water as an example, we have analyzed the role of electronic and inertial asymmetry on NR polarizabilities and hyperpolarizabilities.Preliminary calculations show that the NR isotope effect in this example is comparable in magnitude to that due to ZPVA, but the isotopic pattern is quite different.Although both isotope effects are small for water, they may be larger in other cases.
Apart from making FF calculations feasible, the most important contribution of our treatment is to provide a convenient method for studying the role of individual internal coordinates or sets of internal coordinates.We have shown that the fundamental idea undergirding the FF method can be used as the basis for a systematic approach to determining the most important vibrational motions.In this approach, field-induced coordinates ͑FICs͒ are derived directly from field-dependent geometry optimizations.A specific treatment for -conjugated molecules, where there is a natural direction for the applied field, has been given here.Calculations on three prototype -conjugated molecules reveal that quantitative accuracy can be achieved by means of a 2-FIC treatment.The particular values of the field chosen to generate the two-coordinate basis do not seem to be important, but to avoid ambiguity we suggest the mid-range values F ϭ0.0016/0.0032a.u. in all cases.
From the 2-FIC treatment, one can extract a single welldefined natural conjugation coordinate ͑NCC͒.This NCC, which for these prototypical conjugated molecules is more accurate than the a priori ECC, provides semiquantitative NR properties where applicable, and also yields a physically intuitive picture.
We have initiated further studies based on the 2-FIC method.One such investigation relates to a recent paper 10 demonstrating that the two-state valence-bond chargetransfer ͑VB-CT͒ model is not valid, in its current formulation, for NR polarizabilities and hyperpolarizabilities.This conclusion was based on trial calculations for a set of ten molecules that cover the entire range of VB-CT mixing.The very same set of molecules should pose a strenuous test of our new 2-FIC treatment.We hope that this study will eventually lead to a new VB-CT model that will be both practical and generally applicable.
From the viewpoint of methods, as opposed to applica- tions, we are working on a more detailed analysis of the 2-FIC treatment.The possibility that a two-coordinate FF procedure can be combine with perturbation theory in a way that achieves ''the best of both worlds'' is also under investigation.

FIG. 1 .
FIG. 1.The a priori effective conjugation coordinate ͑ECC͒.All arrows are drawn to the same scale.

FIG. 2 .
FIG. 2. The natural conjugation coordinate ͑NCC͒.For each molecule the transverse and longitudinal components are shown separately.All arrows are drawn to the same scale except as indicated.

TABLE I .
Nuclear relaxation ͑NR͒ polarizabilities and hyperpolarizabilities ͑in a.u.͒ of prototype -conjugated molecules calculated using the one vibrational degree of freedom SFIC model.
aThe percentage error compared to the complete 3N-6 coordinate treatment of TableIIis given in parentheses.b Infinite frequency approximation.

TABLE II .
Reference values ͑in a.u.͒ for nuclear relaxation ͑NR͒ polarizabilities and hyperpolarizabilities of prototype -conjugated molecules obtained using a complete set of 3N-6 coordinates.

TABLE III .
Nuclear relaxation ͑NR͒ polarizabilities and hyperpolarizabilites ͑in a.u.͒ of prototype -conjugated molecules calculated using the one vibrational degree of freedom WFIC model.
aThe percentage error compared to the complete 3N-6 coordinate treatment of TableIIis given in parentheses.b Infinite frequency approximation.

TABLE IV .
Nuclear relaxation ͑NR͒ polarizabilities and hyperpolarizabilities ͑in a.u.͒ of prototype -conjugated molecules calculated for the 2-FIC treatment with Fϭ0.0016/0.0032a.u.The percentage error compared to the complete 3N-6 coordinate treatment of Table II is given in parentheses.
a b Infinite frequency approximation.

TABLE V .
Nuclear relaxation ͑NR͒ polarizabilities and hyperpolarizabilities ͑in a.u.͒ of prototype -conjugated molecules calculated using the one vibrational degree of freedom NCC model.
aThe percentage error compared to the complete 3N-6 coordinate treatment of TableIIis given in parentheses.b Infinite frequency approximation.

TABLE VI .
Nuclear relaxation ͑NR͒ polarizabilities and hyperpolarizabilities ͑in a.u.͒ of prototype -conjugated molecules calculated from the one coordinate a priori ECC model.
aThe percentage error compared to the complete 3N-6 coordinate treatment of TableIIis given in parentheses.b Infinite frequency approximation.

TABLE VII .
Nuclear relaxation ͑NR͒ polarizabilities and hyperpolarizabilities ͑in a.u.͒ for H 2 O, D 2 O, and HDO.The frequency-dependent properties are obtained in the infinite frequency approximation.for HDO is almost exactly the average of the corresponding H 2 O and D 2 O values.Furthermore, while the magnitude of the ZPVA term decreases by & ͑Refs.49,50͒ from H 2 O to D 2 O, the NR contribution typically increases.According to the numerical values in Table

TABLE VIII .
Nuclear relaxation ͑NR͒ and zero-point vibrational averaging ͑ZPVA͒ contributions ͑in a.u.͒ to polarizabilities and hyperpolarizabilities of H 2 O, D 2 O, and HDO.The frequency-dependent properties are obtained in the infinite frequency approximation.See Ref. 50 for definition of polarizabilities; Ref. 5 for definition of hyperpolarizabilities. a