Determination of vibrational polarizabilities and hyperpolarizabilities using field-induced coordinates

An analytical set of field-induced coordinates ~FICs! is defined. It is shown that, instead of 3 N 26 normal coordinates, a relatively small number of FICs is sufficient to describe the vibrational polarizability and hyperpolarizabilities due to nuclear relaxation. The fact that the number of FICs does not depend upon the size of the molecule leads to computational advantages. A method is provided for separating anharmonic contributions from harmonic contributions as well as effective mechanical from electrical anharmonicity. Hartree–Fock calculations on a dozen representative conjugated molecules illustrate the procedures and indicate that anharmonicity can be very important. Other potential applications including the determination of zero-point vibrational averaging corrections are noted. © 2000 American Institute of Physics. @S0021-9606 ~00!32137-7#


I. INTRODUCTION
The key role of vibrational contributions to the nonlinear optical ͑NLO͒ properties of conjugated polymers and organic molecules of practical consequence is now well established. 1,2At the microscopic level these properties are governed by the first hyperpolarizabilily tensor, ␤, and the second hyperpolarizability tensor, ␥.Typically, the longitudinal component of these tensors (␤ L ,␥ L ) will be dominant.Most treatments of the vibrational ␤ L and ␥ L have been carried out by the perturbation theory method of Bishop and Kirtman 3-5 ͑see, for example, Refs.6-8͒.Even at the simplest ͑double harmonic͒ level of approximation, however, ab initio calculations for large molecules have been limited to the Hartree-Fock level because a complete force constant determination is required.When electrical and mechanical anharmonicity is included the computational difficulties are exacerbated.
A few years ago Bishop, Hasan, and Kirtman ͑BHK͒ ͑Ref.9͒ formulated an alternative finite field method that does not require explicit calculation of the vibrational force constants.It utilizes an ''infinite optical frequency'' approximation which implicity includes low order anharmonicity terms.The BHK procedure gives the so-called nuclear relaxation contribution to the vibrational polarizability and hyperpolarizabilites. Except for zero-point vibrational averaging ͑ZPVA͒ this contribution contains the leading vibrational term of each type ͑see later͒.Tests of the infinite optical frequency approximation 10,11 have shown that it yields satis-factory results, at least for small molecules.The generalization of BHK for an arbitrary optical frequency has recently been presented 12 but it is more difficult to apply.
The first step of BHK involves determining the relaxation of the equilibrium configuration in response to a finite external static field.This is followed by calculation of the change thereby induced in the electronic dipole moment, linear polarizability, and first hyperpolarizability.Finally, the desired vibrational properties are obtained from these fieldinduced property changes by numerical differentiation with respect to the field.The first successful implementation of this procedure, which requires a careful treatment of the Eckart conditions, was reported a short while ago. 13Since then it has been extended to the static linear polarizability of polymers 14 and, very recently, 15 applied to obtain the vibrational ␥ L for eight different homologous series of conjugated oligomers, each containing up to twelve heavy atoms along the backbone.In addition, it has been shown 16 that an exactly analogous procedure can be used to determine nuclear relaxation corrections to the ZPVA.Despite these successes further improvements are desirable.One difficulty is that repeated geometry optimizations are necessary in order to determine the numerical derivatives and, for sufficiently accurate results, we find 15 that very tight thresholds on the residual forces have to be employed.As a consequence it has not yet been feasible to study the role of electron correlation for systems of interest in materials science.8][19] Furthermore, in contrast with perturbation theory, BHK does not yield the contribution of individual coordinates nor does it allow a separation of electrical and mechanical anharmonicity effects.
In this paper we present a combination of the perturbation theory and BHK ͑i.e., nuclear relaxation͒ approaches that paves the way for calculations which include electron correlation and, at the same time, allow one to specify the key vibrational coordinates and to separate electrical from mechanical anharmonicity effects.Our new approach is based on the determination of a small number of vibrational coordinates that exactly reproduce the BHK results for the complete vibrational space.There are quite a few calculations in the literature which show that, for longitudinally extended conjugated molecules or oligomers, 1,7,20 the double harmonic vibrational polarizabilities and hyperpolarizabilities are often dominated ͑at least, at the Hartree-Fock level͒ by contributions from a limited number of normal modes.However, no prescription has been given for determining the precise form of these modes or for predicting when other modes may become important.On the other hand, in our initial implementation of BHK ͑Ref.13͒ we learned that one could numerically generate a pair of field-induced coordinates ͑FICs͒ that were sufficient to accurately determine the vibrational nuclear relaxation ͑infinite frequency approxima-tion͒ longitudinal polarizabilities and hyperpolarizabilites for a small set of prototype -conjugated molecules.This discovery provided the impetus for the development of exact field-induced coordinates which provide the desired prescription.
The FICs are analytically defined in Sec.II A. These coordinates are associated with the displacement of the equilibrium geometry due to a static applied field.Corresponding to the displacements that are first-order, second-order, etc. in the field there are first-order, second-order, etc. FICs.In each order, beyond first, one can also define a pure harmonic coordinate by discarding the anharmonic component.For the longitudinal nuclear relaxation ͑infinite optical frequency͒ ͑hyper͒polarizabilities either one or, at most, two of these FICs are sufficient to reproduce exactly the results of a complete 3N-6 normal coordinate calculation.After proving these relations for the longitudinal component we generalize to the other diagonal, as well as off-diagonal, components of the polarizability and hyperpolarizability tensors.More FICs are needed, of course, to generate the entire tensor but the key point is that the number of coordinates remains the same regardless of the size of the system.Furthermore, our results are not limited to conjugated systems; they are valid for any molecule.
In Sec.II B the FIC formulas for evaluating the vibrational nuclear relaxation properties are presented along with an analysis of various computational methods-both ''exact'' and approximate-that can be employed to take advantage of the reduced number of active coordinates.The results of Hartree-Fock calculations on 12 representative conjugated molecules are reported in Sec.III.These serve to ͑1͒ evaluate some of the approximation methods; ͑2͒ determine the form of the FICs; ͑3͒ illustrate their use for interpretive purposes; and, by way of passing, ͑4͒ assess the potential importance of anharmonic effects.Finally, in the last section we discuss electron correlation calculations and other future applications including an extension to the zero-point vibrational averaging correction ͑as well as the effect of nuclear relaxation on the latter͒ and to determining deviations from the infinite optical frequency approximation.

A. Definition and properties
In the presence of a uniform static external electric field the equilibrium geometry of a molecule will relax to a new field-dependent configuration.This is due to the fact that the electrostatic interaction with the field depends linearly upon the field-free normal coordinates.The new minimum in the potential energy may be obtained as usual by applying the stationary condition.One may solve the resulting relation order-by-order in the field 21,22 for the value of the ith fieldfree normal coordinate at the field-dependent equilibrium geometry, where, and the sums over 3NϪ6 normal coordinates reduce to 3N Ϫ5 for linear molecules.In Eq. ͑2͒, n is the total order of differentiation with respect to normal coordinates, while m is the total order with respect to the fields.Derivatives with mϾ0 and nϾ1 are electrical anharmonicity parameters.For mϭ0 those derivatives with nϾ2 are the usual mechanical anharmonicity constants.The quantity q 1 i,a in Eq. ͑1͒ represents the linear field-induced change along the field-free normal coordinate Q i .Thus, we define the first-order FIC in the direction a by The second-order FICs are defined analogously,

͑5͒
Here the second term on the rhs of Eq. ͑5͒ contains electrical anharmonicity parameters, while the third term involves mechanical anharmonicity.Using just the first term of Eq. ͑5͒, which is the pure harmonic component, we may also define 2,har Note that 1,har a ϭ 1 a .Although it would be straightforward to extend these definitions to third-and higher-order, that is not necessary for our purposes.Recently, we have become aware of a related treatment of infrared intensity-carrying modes by Torii and co-workers. 23It turns out that those modes are a special case of the first-order FICs defined in Eq. ͑4͒, where the vibrational force constants, i.e., a 20 ii , are the same for all i.This is the appropriate choice for the infrared intensity problem.
We are now in a position to use the first-order FIC in the longitudinal direction ͑L͒ to obtain an analytical expression for the nuclear relaxation contribution to the electro-optic Pockels effect ͑EOPE͒, ␤ L nr (Ϫ;,0) →ϱ .This may be done following a procedure 21 based on BHK.The first step is to define ͕ i ͖ as a set of vibrational coordinates with 1 equal to l L and ͕ 2 , 3 ,..., 3NϪ6 ͖ orthogonal to each other ͑for convenience͒ and to 1 L , i.e., where M is an orthogonal matrix.From Eqs. ͑4͒ and ͑7͒ M 1 j ϭϪq 1 J,L .Obviously, the value of i at the fielddependent equilibrium geometry is given by Then, the field-dependent linear polarizability may be written as a power series in the field FϭF L , in which R 0 and R F are, respectively, the field free and fielddependent equilibrium geometry.In Eq. ͑9͒ the second term, involving the normal coordinate displacements Q i F , gives rise to the nuclear relaxation contribution to the Pockels effect in the infinite optical frequency approximation, i.e., ␤ L nr (Ϫ;,0) →ϱ .Using the chain rule to express ‫␣ץ‬ L /‫ץ‬Q i in terms of ‫␣ץ‬ L ‫ץ/‬ i , the fact that ‫ץ‬Q i F /‫ץ‬F L ϭM 1i , and Eq.͑8͒ for i F , we have This demonstrates that in the infinite frequency approximation a single FIC, i.e., 1 ϭ 1 L , yields exactly the same nuclear relaxation contribution as the complete set of 3N-6 normal coordinates.
Following an analogous treatment, but replacing ␣ L (R F ,FϭF L ) in Eq. ͑9͒ by either L (R F ,FϭF L ) or ␤ L (R F ,FϭF L ), it is easy to demonstrate that ␣ L nr (0;0), and the nuclear relaxation contribution to dc-second harmonic generation ͑dc-SHG͒, ␥ L nr (Ϫ2;,,0) →ϱ , can also be calculated ͑see Table I͒ using only 1 L .On the other hand, TABLE I. Analytical formulas for the polarizabilities and hyperpolarizabilities in terms of field-induced coordinates ͑FICs͒.The upper index on each sum is either 1 or 2, depending upon the particular property and component ͑see Table II͒.For those properties where more than one coordinate is required, we have assumed that a prior linear combination has been made to satisfy the condition ‫ץ‬ 2 V/‫ץ‬ i ‫ץ‬ j ϭ0 for i j. , ii ͬ ␤ L nr (0;0,0) and the nuclear relaxation contribution to the electro-optic Kerr effect ͑EOKE͒, ␥ L nr (Ϫ;,0,0) →ϱ , are determined 9 by the second derivatives of L (R F ,FϭF L ) and ␣ L (R F ,FϭF L ), respectively.Nonetheless, for ␤ L nr (0;0,0), it is easy to show that the 1 L FIC, by itself, gives the full 3NϪ6 normal coordinate result.In this derivation, and those that follow, we may start with the 3NϪ6 normal coordinate perturbation expressions of Bishop and Kirtman [3][4][5] and our Eqs.͑1͒-͑3͒ yields If the chain rule is employed to write ‫ץ‬ 2 L /‫ץ‬Q i ‫ץ‬Q j and ‫ץ‬ 3 V/‫ץ‬Q i ‫ץ‬Q j ‫ץ‬Q k , as well as ‫␣ץ‬ L /‫ץ‬Q i , in terms of the coordinates, i , then Eq. ͑12͒ becomes .

͑13͒
Although the second and third terms on the rhs of Eq. ͑13͒ involve electrical and mechanical anharmonicity only a single anharmonicity parameter occurs in either case.
A procedure parallel to the one above can be employed to obtain a simple FIC equation for ␥ L nr (Ϫ;,0,0) →ϱ .This time we start with Eq. ͑12͒ of Ref. 21, It is convenient in this case to define the FICs so that 1 ϭ 1 L , 2 ϭ 2,har L and the set ͕ 3 , 4 ,..., 3NϪ6 ͖ is orthogo- nal to 1 L and 2,har L ϭ 2,har LL .In that event ͓cf.Eq.
Then, after applying the chain rule to convert from ͕Q i ͖ to ͕ i ͖, taking advantage of orthonormality and where, now, 1 ϭ 1 L and 2 ϭ 2,har L are sufficient to obtain the exact ␥ L nr (Ϫ;,0,0) →ϱ .Alternatively, by combining the terms in Eq. ͑14͒ that contain q 2 one can obtain the expression with the FICs defined so that 1 ϭ 1 L and 2 ϭ 2 L .Thus, the EOKE can also be written entirely in terms of 1 L and 2 L .For ␥ L nr (0;0,0,0) one can follow a similar procedure to show that this quantity depends only on the 1 L and 2 L coordinates.We start with Eq. ͑7͒ of Ref. 21, 3NϪ6 24 ͩ a 31 i jk,L q 1 i,L q 1 j,L q 1 k,L 3NϪ6 24 ͩ a 40 i jkl q 1 i,L q 1 j,L q 1 k,L q 1 l,L ϩ 3a 30 i jk a 21 kl,L a 20 24 ͩ 9a 30 i jk a 30 klm 4a 20 kk q 1 i,L q 1 j,L q 1 l,L q 1 m,L ͪ .

͑20͒
The terms involving a 20 in the denominator can be eliminated by using ͓cf.Eq. ͑1͔͒ 3NϪ6 24 ͩ a 31 i jk,L q 1 i,L q 1 j,L q 1 k,L Ϫ 3 4 This time we employ the vibrational coordinates ͕ i ͖,

͑23͒
Clearly, the only coordinates that appear in Eq. ͑23͒ are 1 ϭ 1 L and 2 ϭ 2 L .Although no static fields are involved in the intensitydependent refractive index ͑IDRI͒ effect, ␥ L nr (Ϫ;, Ϫ,) →ϱ , the vibrational contribution to this property does not vanish in the infinite frequency approximation.To obtain the FIC formula we utilize Eq. ͑14a͒ of Ref. 21, For this case only one FIC, 1 ϭ 2,har L is required.As usual we define ͕ 2 , 3 ,..., 3NϪ6 ͖ so that this set is orthogonal to 2,har

L
. Following a procedure completely parallel to the previous cases Eq. ͑24͒ can be transformed to Thus, ␥ L nr (Ϫ;,Ϫ,) →ϱ can be calculated exactly using only the 1 ϭ 2,har L coordinate.

B. Working analytical formulas for properties and computational considerations
Now that the FICs needed for each property have been defined, we are ready to develop compact analytical working expressions for these properties.This is done by expanding the potential energy in terms of the required FICs.It is convenient to diagonalize the Hessian in the reduced FIC basis.Then, one can use exactly the same formulation as in the usual scheme based on normal coordinates. 21The resulting equations are presented in Table I in a form that extends the above treatment to all components of the polarizability and hyperpolarizability tensors.For each diagonal component of the property one needs the FICs in the corresponding direction; for off-diagonal components the three independent offdiagonal components of the second-order FICs may be needed as well.Thus, from Table II we see that 15 FICs are required to obtain all components of all properties.That number is reduced to 9 if only the diagonal components are desired or 6 if the IDRI is excluded.Still further reductions occur for various subsets of the properties and/or components.In order to determine the longitudinal Pockels effect, for example, just a single FIC is required.For average values, however, each component must be computed separately.
A key point is that the number of FICs involved is independent of the size of the system.Note that the same fundamental quantities appear in Table I as in Eqs.͑2͒ and ͑3͒ except that the derivatives are taken with respect to FICs as opposed to normal coordinates.In fact, we do not use the FICs themselves but, rather, the linear combinations that diagonalize the Hessian.
An immediate question that arises is how to take advantage of the reduction in the number of coordinates upon transforming from normal modes to FICs.From a computational perspective two distinct phases are involved in our treatment.The first is a determination of the FICs.Both 1 a and 2 ab can be obtained either by evaluating appropriate derivatives of the potential or from field-dependent geometry optimizations.The first-order FICs depend upon the dipole derivatives ‫‪Q‬ץ/ץ‬ i ϭϪ‫ץ‬ 2 V/‫ץ‬F‫ץ‬Q i and the vibrational force constants ͑i.e., the Hessian͒.2,har ab is similar except that it involves the derivative of ␣ rather than ; whereas 2 ab depends, in addition, upon anharmonicity parameters determined by one further derivative with respect to the normal coordinates.Even at the Hartree-Fock level, to our knowledge there are no commonly available programs for analytically evaluating all these anharmonicity parameters.On the other hand, 2 ab is readily obtained from finite field geometry optimizations-as is 1 a -and this procedure also avoids calculating the full 3NϪ6ϫ3NϪ6 Hessian.Since only loworder FICs are required the fields should be kept small, which is a feature we did not utilize in earlier work. 13Small fields usually, though not always, make it easier ͑fewer steps͒ to attain the tight convergence in geometry evaluation that was previously 15 found necessary to obtain accurate results ͑primarily for static ␥͒ using the BHK treatment.Another possibility ͑not pursued here͒ is to forego tight convergence and, instead, employ one or two additional FICs based on semiempirical calculations and/or intuitive considerations.Unfortunately, 2,har ab , which is needed only for the IDRI, cannot be obtained from field-dependent geometry optimizations.It can, of course, be found simply by evaluating the polarizability derivatives ‫␣ץ‬ ab /‫ץ‬Q i and the harmonic force constants.If one wishes to avoid ab initio computation of the Hessian a reasonable approach is to employ semiempirical force constants and normal coordinates.In order to improve on the accuracy in this case one could augment the basis with ͑a small number of͒ FICs as mentioned above.
The second phase of our treatment is to evaluate the nuclear relaxation polarizabilities and hyperpolarizabilities in terms of the FICs.Again, this may be done either by explicit evaluation of the derivatives involved or by the BHK finite field procedure.The former yields each of the individual contributions shown in Table I, which is desirable for interpretive purposes.It is also the only way to determine the IDRI and, for smaller molecules, especially if one is interested in all tensor components, it may be the computationally more efficient method for all properties.As noted above there are no commonly available programs that permit analytical evaluation of the anharmonicity parameters in Table I.For large molecules the number of parameters is dramatically lessened by replacing normal coordinates with FICs, which makes their numerical computation feasible.Indeed, even at levels where the analytical Hessian is not available, the use of the FICs with the BHK method could be advantageous.

III. TEST CALCULATIONS
Ab initio RHF/6-31G calculations were carried out for testing purposes on hexatriene, hexasilane, and a representative set ͑see Fig. 1͒ of ten push-pull -conjugated molecules taken from the work of Bishop et al. 24 Although average values of the various properties were calculated to confirm the efficiency of the FIC procedure, our focus here will be on the longitudinal component of the hyperpolarizability tensors, which is dominant for these molecules.We took the longitudinal direction to be along the principal axis associated with the largest rotational constant.The first-order FIC, 1 L , was determined both analytically and by finite field geometry optimization.Since 2,har L cannot be obtained by the latter procedure we used the analytical coordinate ͑see below however͒.On the other hand, it is excessively timeconsuming to compute 2 L analytically so that coordinate was found by the finite field geometry optimization technique.Given the FICs we want to evaluate the properties using the derivative expressions of Table I.For our molecules it is feasible at the Hartree-Fock/6-31G level to use the GAUSS- IAN98 suite of programs 25   employed, together with derivatives obtained as above, to compute the properties listed in Table I ͑except for the static ␥ which requires 2 L ͒.For comparison the same properties were calculated using all 3NϪ6 normal coordinates.Because the normal coordinate treatment is so time-consuming we did not evaluate the derivatives in this case but, instead, applied the equivalent finite field procedure which is embodied in our recently developed Eckart program. 13As expected the 3N-6 normal coordinate results were reproduced by the FICs for all molecules and for all properties within 0.3% numerical round-off error.
Next we obtained 1 L and 2 L by finite field geometry optimization with Fϭ0.0,Ϯ0.0004a.u. and repeated the FIC calculations with these numerical coordinates.Remarkable agreement with the analytical coordinate values-within 1.5% in every case-was found.Thus, the field-dependent geometry optimization method for determining the first-and second-order FICs appears to be very robust.The IDRI was not included in the data set because it is determined by 2,har L rather than 1 L and 2 L .For the molecules I, II, III, IX, and XII it turns out that the 1 L , 2 L pair gives good values of ␥ L nr (Ϫ;,Ϫ,) →ϱ , but that is not true for the other donor/acceptor molecules or for hexasilane.It is interesting to note that the donor/aceptor molecules showing a small effect due to the anharmonicity in 2 L are covalent, whereas the remaining molecules are either of medium or large ionicity. 24The possible relationship between the importance of anharmonicity and the degree of covalency is a subject that deserves further investigation.We tried to improve the nuclear relaxation values of the IDRI simply by using an additional FIC generated from a geometry optimization carried out at a relatively large field, i.e., Fϭ0.0064 a.u.However, the results obtained showed no improvement.This indicates that it is probably preferable to approximate 2,har L directly, as discussed in Sec.II B, rather than by purifying 2 L of the anharmonic contribution.
The harmonic FICs 1 L and 2 L are displayed in Figs.2-7.For the two centrosymmetric C 2h molecules, namely, Si 6 H 14 ͑XI͒ ͑Fig.6͒ and C 6 H 8 ͑XII͒ ͑Fig.7͒, these FICs are symmetry coordinates-1 L has the same symmetry (b u ) as L while 2,har L has the same symmetry (a g ) as ␣ L .As expected 1 L bears a similarity to the normal mode that makes the major contribution to ͓ 2 ͔ 0 , i.e., the TAM-2 mode in C 6 H 8 ͑Ref.26͒ and the H-wagging mode in Si 6 H 14 . 7,8Likewise, 2,har L is clearly related to the normal coordinates that FIG. 1. Formula/structure of molecules studied in this paper.͑For molecules I-X and XII bonds of length 1.38-1.40Å are considered to be intermediate between single and double bonds and are indicated by dashed lines.͒dominate ͓␣ 2 ͔ 0 , in this case the LAM-1 mode for Si 6 H 14 and a carbon skeletal motion for C 6 H 8 .The latter is obtained from two ECC-type normal modes where the skeletal carbon motion is combined in-phase and out-of-phase with the H-wag.It should be noted that the natural conjugation coor-dinate ͑NCC͒ defined in our earlier paper 13 is not a single FIC but rather a composite of ͑primarily͒ 1 L and 2 L .All the remaining compounds are -conjugated donoracceptor ͑D/A͒ molecules.We have chosen I ͑Fig.͑Fig.4͒ as representative of those that have a polyenic or Zwitterionic structure whereas VI ͑Fig.3͒ and X ͑Fig.5͒ are cyaninelike.In the former case both 1 L and 2,har L exhibit substantial bond length alternation ͑BLA͒ character-the aromatic↔quinoid motion in p-nitroaniline ͑IX͒ is an example-as well as H atom displacements.For the molecules with a cyanine structure other longitudinal, as well as transverse, displacements are more important, and there is also a non-negligible out-of-plane component ͑not shown in Figs. 3 and 5͒.These results confirm our previous conclusion 13 that no two-state valence bond-charge transfer model based on a single BLA coordinate can successfully describe vibrational NLO properties of D/A -conjugated organic chains.
Finally, we turn to the separation of anharmonic from harmonic contributions as well as electrical anharmonicity from mechanical anharmonicity.There are three nuclear relaxation properties that are influenced by anharmonicity, namely, ␤ L nr (0;0,0), ␥ L nr (Ϫ;,0,0) →ϱ , and ␥ L nr (0;0,0,0).Tables III and IV provide a breakdown of the values of these properties into the various terms as classified by perturbation theory.27 In the case of ␤ L nr (0;0,0), for example, it is readily seen that the first term on the rhs of the expression in Table I is the ͑double harmonic͒ ͓␣͔ 0,0 term in the notation of Refs.3-5.The second term on the rhs, which is first-order in electrical anharmonicity, is ͓ 3 ͔ 1,0 , whereas the third term is first-order in mechanical anharmonicity and denoted by ͓ 3 ͔ 1,0 .For ␥ L nr (Ϫ;,0,0) →ϱ the expression in Table I takes the same form whether one uses 2 L or 2,har L along with 1 L .However, in order to assign a specific order of perturbation theory to each term it is necessary to employ 2,har L .Then, the five successive terms on the rhs are ͓␤͔ 0,0 ,͓␣ 2 ͔ 0,0 , the two components of ͓ 2 ␣͔ 1,0 , and ͓ 2 ␣͔ 0,1 .For ␥ L nr (0;0,0,0) a second calculation is required where 2,har L is used in the formula of Table I rather than 2 L .This gives the correct value for the first five terms in addition to one component of ͓ 4 ͔ 2,0 ͑i.e., the sixth term͒ and one component of ͓ 4 ͔ 0,2 ͑i.e., the eighth term͒.Combining this information with the 2 L result; which gives the total ͓ 4 ͔ contribution and with Eq. ͑22͒, we obtain the combinations ͓ 4 ͔ 2,0 ϩ1/2͓ 4 ͔ 1,1 and ͓ 4 ͔ 0,2 ϩ1/2͓ 4 ͔ 1,1 .It is not possible by means of FICs to make a separation of all three second-order components, nor is it possible with any coordinates to make a clean separation of electrical and mechanical anharmonicity.On the other hand, the second-order combinations given above would seem to be the most reasonable choice for splitting the total ͓ 4 ͔ contribution into effective mechanical and electrical anharmonicity components.
The purpose of our calculations, which are reported in Tables III and IV, is to show how the FICs can be utilized for interpretive purposes.In contrast with our previous findings 13 for planar -conjugated oligomers, we note that anharmonicity plays a major role in determining all three properties.For the EOKE this is especially true when the ͓␣ 2 ͔ 0,0 and ͓␤͔ 0,0 terms nearly cancel as in compounds V and XI.Typically, for ␥ L nr (0;0,0,0) the anharmonicity is dominant due to the ͓ 4 ͔ electrical anharmonicity term.͑Note that terms of the same type, e.g., ͓␣ 2 ͔ 0,0 , have different values in Tables III and IV because they occur with different coefficients.͒The above observations, though indica-tive, must be qualified by the fact that we have used a small basis and omitted electron correlation.

IV. CONCLUSIONS AND FUTURE WORK
By transforming from normal coordinates to fieldinduced coordinates ͑FICs͒ we have shown that the vibrational degrees of freedom required to completely describe nuclear relaxation polarizabilities and hyperpolarizabilities is reduced from 3N-6 to a relatively small number which does not depend upon the size of the molecule.A summary of the FICs required for each property is given in Table II.In terms of normal coordinates the number of elements of the Hessian that must be determined, even for the lowest ͑double har-monic͒ level of calculation, scales as (3NϪ6) 2 .If anharmonic effects are included-and our calculations indicate they may be quite important in donor/acceptor systems for EOKE, as well as the static ␤ and ␥-then the number of anharmonicity parameters scales as (3NϪ6) 3 and (3NϪ6) 4 for the first and second hyperpolarizabilities, respectively.Various schemes for taking computational advantage of the zeroth-order scaling in terms of FICs, based in large part on field-dependent geometry optimizations, have been dis-   cussed.Since these optimizations are carried out in a reduced coordinate space they can be done much more efficiently than in the space of 3NϪ6 normal coordinates.It has also been shown that the vibrational contributions can be separated into different types, as in perturbation theory, for interpretive purposes.When more than one FIC is employed a breakdown into the individual degrees of freedom and their interactions could prove useful.With further study we are hopeful that a practical intuition as to the nature of the FICs will develop.There are several possibilities for applying FICs to obtain vibrational hyperpolarizabilities beyond the →ϱ nuclear relaxation approximation.One of these is to account for finite optical frequencies either approximately, by using the →ϱ FICs in the exact nuclear relaxation perturbation theory formulas, or by generating exact finite frequency FICs.Another is to determine the zero-point vibrational average ͑ZPVA͒ of the hyperpolarizability.This contribution is difficult to compute and, for that reason, it is not yet known how important it is for large organic molecules of interest in nonlinear optics.The ZPVA of any property can be divided into two terms 2 -one due to mechanical anharmonicity and the other due to electrical anharmonicity.In a forthcoming paper 28 we show that both can be written compactly in terms of the zero-point vibrational energy.The contribution due to mechanical anharmonicity, then, depends only on a single FIC, whereas the electrical anharmonicity term does not depend explicitly on vibrational coordinates.As in the case of nuclear relaxation, the use of FICs to determine the ZPVA provides an important simplification.

FIG. 3 .FIG. 4 .
FIG. 2. Atomic displacements for the FICs 1 L ͑top͒ and 2,har L ͑bottom͒ of molecule I from Fig. 1.The length of the arrow is proportional to the displacement.

FIG. 5 .FIG. 6 .
FIG. 5. Atomic displacements for the FICs 1 L ͑top͒ and 2,har L ͑bottom͒ of molecule X from Fig. 1.The length of the arrow is proportional to the displacement.
or, equivalently, with those of Ref. 21.Using Eq. ͑6͒ of Ref. 21, i.e., where 1 ϭ 1 L , 2 ϭ 2 L , and ͕ 3 , 4 ,..., 3NϪ6 ͖ are orthogonal to 1 L and 2 L .In that event ‫ץ‬Q i F /‫ץ‬F L ϭM 1i , ‫ץ‬ 2 Q i F /‫ץ‬F L 2 ϭM 2i and it is straightforward to verify that Eqs.͑15͒-͑17͒ remain valid ͑although now 2 and M 2i have a different meaning͒ if we replace ‫ץ(‬ 2 i F /‫ץ‬F L 2 ) har by ‫ץ‬ 2 i F /‫ץ‬F L 2 .Then, applying the chain rule to transform a 13 i,LLL ,a 12 i,LL ,... from ͕Q i ͖ to ͕ i ͖ one obtains to obtain analytical results for a 20 , a 01 , a 11 , a 02 , a 12 , and a 03 .Then, numerical differentiation of a 20 , a 11 , a 12 , and a 03 with respect to the FICs yields a 30 , a 21 , a 22 , and a 13 , respectively.a 40 and a 31 were computed by double numerical differentiation of a 20 and a 11 .
As a first step the analytical FICs 1 L and 2,har L were

TABLE II .
FICs required to calculate diagonal and off-diagonal elements of the property tensor.The directions a, b, c, and d are along the x, y, and z Cartesian axes.

TABLE III .
Breakdown of RHF/6-31G ␤ L →ϱ for the molecules of Figs.2-7into the harmonic and anharmonic terms defined in Ref.3.All the values are given in a.u.