Variational Calculation of Static and Dynamic Vibrational Nonlinear Optical Properties

The vibrational configuration interaction method used to obtain static vibrational ͑hyper͒polarizabilities is extended to dynamic nonlinear optical properties in the infinite optical frequency approximation. Illustrative calculations are carried out on H 2 O and NH 3. The former molecule is weakly anharmonic while the latter contains a strongly anharmonic umbrella mode. The effect on vibrational ͑hyper͒polarizabilities due to various truncations of the potential energy and property surfaces involved in the calculation are examined.


I. INTRODUCTION
Nonlinear optical properties ͑NLOP͒ play a key role in a number of fields including materials science, communications, and medicine. 1For that reason there is much practical, and fundamental, interest in the ab initio prediction of these properties, [2][3][4][5] which characterize the response of a material to an external electrical field.At the molecular level, NLOP are governed by the first and second hyperpolarizabilities ͑␤ and ␥͒.Since the nuclei of molecules, as well as the electrons, respond to the field, there will be contributions due to both types of motion.Although a complete treatment would require a fully nonadiabatic calculation, 6 in this paper we invoke the Born-Oppenheimer approximation ͑BOA͒.A nonadiabatic treatment is currently feasible only for diatomic molecules and a comparison with the BOA for NLOP of LiH and LiD did not show significant differences. 7In fact, in the past it has often been assumed that the vibrational contribution as well could be neglected in comparison to the electronic term.5][26] The BK method is based on a double-harmonic zeroth-order approximation which takes into account the harmonic vibrational potential energy surface ͑PES͒ and electrical field interaction terms that are linear in the normal modes.In order to include all significant mechanical and electrical anharmonicity contributions, it has been found necessary to include fourth-and higher-order derivatives of the vibrational PES, as well as the electrical properties, with respect to the normal coordinates.This rapidly becomes impractical as the size of the molecule increases.
8][29][30] This finite field ͑FF͒ "nuclear relaxation" approach yields the vibrational contributions to static and to electro-optical hyperpolarizabilities, which are defined as changes in the dynamic polarizability and first hyperpolarizability induced by a static electric field.Three well-known examples are the dc-Pockels effect ͑dc-P͒, the Kerr effect, and electric field induced second harmonic generation ͑ESHG͒. 1 Although it is not required, the FF nuclear relaxation method is usually applied in conjunction with the infinite optical frequency ͑IF͒ approximation.The latter assumes that the optical frequencies are infinite in comparison to all vibrational frequencies.2][33] On the other hand, this approxima-tion is not valid for optical frequencies in or near the IR region. 31s might be expected, the BK perturbation treatment is unreliable for very anharmonic systems such as weakly bound dimers. 34,357][38] For these cases a variational treatment of the vibrational effects is desirable.As the first step, a variational procedure for solving the vibrational Schrödinger equation is required.Much progress in that direction has been made in recent times.The procedures developed are analogous in many respects to the variational methods used in electronic structure theory.Thus, the simplest approach is the vibrational self-consistent field ͑VSCF͒ method [39][40][41][42] in which each normal mode vibrates in the average field generated by all the other modes.The remaining correlation corrections due to mode-mode anharmonicity can, then, be added on by the vibrational analog of Møller-Plesset perturbation theory, [43][44][45] configuration interaction, [46][47][48] or coupled cluster techniques 49 ͑see Ref. 50 for a recent review͒.In this paper we use vibrational configuration interaction ͑VCI͒.
Given the VSCF and VCI wave functions, one can evaluate the static vibrational hyperpolarizabilities in several different ways.Very recently, Christiansen and co-workers showed how response theory methodology can be used to calculate static and frequency-dependent vibrational linear polarizabilities. 16,17,51Their approach can be extended to hyperpolarizabilities as well.In an earlier work by Torrent-Sucarrat et al. 18 static vibrational hyperpolarizabilities were determined variationally using the FF nuclear relaxation method.Here we extend the latter approach to obtain electrooptical hyperpolarizabilities in the IF approximation.
In Sec.II we summarize the relevant background theory, and then, in Sec.III the computational details are described.This is followed in Sec.IV by the presentation of sample calculations for water and ammonia, with special attention paid to the role of the ammonia umbrella motion.Finally, our conclusions and some future plans are given in Sec.V.

A. FF nuclear relaxation approach
Although it has been shown that the BK perturbation theory and FF nuclear relaxation formulas have an exact correspondence, 28 the two methods divide the total vibrational contribution to the ͑hyper͒polarizability differently.In the BK method the total value is the sum of the zero-point vibrational averaging correction to the electronic property ͑P zpva ͒ plus a remaining so-called pure vibrational term ͑P ͒.On the contrary, in the FF nuclear relaxation approach the vibrational ͑hyper͒polarizability is the sum of a nuclear relaxation term ͑P nr ͒ and a curvature contribution ͑P c ͒. P nr is the change in the pure electronic ͑hyper͒polarizability caused by field-induced nuclear relaxation; P c is the sum of P zpva plus P c−zpva , where the latter is the change in P zpva caused by field-induced nuclear relaxation.Thus, P nr + P c−zpva = P . 28In order to obtain the total vibrational correction one must add P zpva to either side of this relation.
P nr contains the lowest-order anharmonicity terms of P , whereas P c−zpva contains all the remaining higher-order terms.Simple analytical formulas are available for P nr . 28owever, the formulation of P c−zpva is far more complex.In this paper P c−zpva is obtained from P zpva by the numerical FF nuclear relaxation method of Kirtman, Luis, and Bishop ͑KLB͒. 30For satisfactory performance it is necessary to determine P zpva accurately.In this regard we use VSCF and VCI to calculate anharmonic wave functions and vibrationally averaged ͑hyper͒polarizabilities 17 ͑see Secs.II B, II C, and III͒.
According to KLB, P c−zpva may be obtained, in the IF approximation, from the coefficients in a power series expansion of ⌬P zpva = P zpva ͑R F , F͒ − P zpva ͑R 0 ,0͒ as a function of the static electric field F. We use R F here to indicate the relaxed equilibrium geometry in the presence of the field.The relevant property expansions are given by and where The subscripts ␣, ␤, ␥, and ␦ refer to the Cartesian axes and → ϱ denotes the IF limit.The accuracy of P c−zpva depends only on the accuracy of the P zpva values.

B. VSCF and VCI methods
In this paper we focus on the role of vibrations, while ignoring rotations and rotation-vibration coupling.We also neglect the mass-dependent terms in the effective vibrational potential of the full Watson Hamiltionian. 52,53Under these approximations, the vibrational Hamiltonian is given by where M is the number of vibrational modes ͑i.e., 3N −5 or 3N −6͒ and Q i is a normal coordinate.We expand V ˆ͑Q 1 , Q 2 , ... ,Q M ͒ as a Taylor series in the normal coordinates.Then, the vibrational Hamiltonian can be expressed as a linear combination of products of one-mode operators, where the one-mode operator h ˆi k ͑Q i ͒ may be one of three types: ‫ץ‬ ˆ2 / ‫ץ‬Q i 2 , Q ˆi n , or just the unit operator.The first M terms in the sum over k on the right hand side of Eq. ͑11͒ correspond to the kinetic energy operators whereas the remaining terms correspond to V ˆ͑Q 1 , Q 2 , ... ,Q M ͒.
In the VSCF methodology the ground state vibrational wave function 0 This wave function is exact in the case of the harmonic approximation for V͑Q 1 , Q 2 , ... ,Q M ͒ or similarly separable potentials.When interaction between modes is included, one can obtain optimal modals by solving the Hartree-type equation where the effective potential is Equations ͑13͒ and ͑14͒ must be solved self-consistently since the evaluation of V i 0 ͑Q i ͒ depends on the modals.In the implementation used in this paper each modal is expanded in a set of harmonic oscillator eigenfunctions for that mode. 54he VSCF vibrational energy E VSCF is just the expectation value of the Hamiltonian, The energy is given as a sum of products of one-mode integrals, which are easily evaluated if the modals are expanded in terms of harmonic oscillator functions.
As explained in the Introduction, the VSCF method incorporates anharmonic coupling between different modes only in an average sense.In this paper dynamical correlation due to mode-mode anharmonicity is incorporated using VCI.The VCI wave function is a linear combination of the ground state VSCF wave function and excited state wave functions generated by replacing occupied modals of the ground state with virtual ͑i.e., unoccupied͒ modals,

͑16͒
Here is a compound index specifying an excitation out of the VSCF state in a specific number of modes to a specific set of levels for these modes resulting in .Optimum values for the coefficients C 0 and C are obtained by applying the variational principle.The full VCI wave function, which is exact for a given vibrational Hamiltonian and harmonic oscillator basis set, is obtained when all possible excitations are included.The molecules considered in this work contain six or less degrees of freedom making full vibrational configuration interaction ͑FVCI͒ quite feasible.For somewhat larger systems other wave function parametrizations may be more efficient than the VCI, for example, vibrational coupled cluster methods ͑VCC͒. 49n this paper we use the direct VCI implementation of Christiansen, 49,54 wherein the VCI Hamiltonian is never explicitly constructed.The VCI states are generated using an iterative method, in which trial vectors span a reduced space and the VCI eigenvalue problem is solved in that space.New trial vectors based on the solution vector within a given reduced space are created using a Davidson update. 55

C. Calculation of vibrational average properties
The zero-point vibrational averaging ͑ZPVA͒ correction to the electronic property P is just the difference between the average value of the property in the ground vibrational state and the value at the equilibrium geometry P e , As in the case of V ˆ͑Q 1 , Q 2 , ... ,Q M ͒, we expand P ˆ0 e ͑Q 1 , Q 2 , ... ,Q M ͒ as a Taylor series in the normal modes and write the result as a linear combination of products of one-mode operators, Since the one-mode operators are of the same type as that occurring in the Hamiltonian, the zpva calculation uses the same algorithms as the energy evaluation.Thus, the VSCF wave function algorithms, including the direct VCI procedure described in Ref. 17, can be used straightforwardly.

III. COMPUTATIONAL DETAILS
In order to illustrate the method described above the diagonal vibrational ͑hyper͒polarizability component parallel to the dipole moment has been computed for water and ammonia.][59][60] Vibrational ͑hyper͒polarizabilities were determined by the FF nuclear relaxation method using the fields of ±0.0002, ±0.0004, ±0.0008, ±0.0016, ±0.0032, ±0.0064, ±0.0128, ±0.0256, and ±0.0512 a.u.The most accurate derivatives were selected from a Romberg method triangle. 62Fielddependent geometry optimizations were done with our own program which rigorously enforces the Eckart conditions. 63t each optimized field-dependent geometry we obtained Taylor series expansions in the normal modes through sixth order for the potential energy, fifth order for the dipole moment, and fourth order for the static linear polarizability and first hyperpolarizability.At the field-free equilibrium geometry we also obtained the fourth-order Taylor series expansion of the static second hyperpolarizability.Third-and higher-order potential energy derivatives were calculated from analytical quadratic force constants by numerical differentiation.Similarly, the higher-order dipole moment derivatives were calculated from the analytical first derivative and the ͑hyper͒polarizability derivatives from the ͑hyper͒polarizability itself.
In the VSCF and VCI calculations we used a basis set for each mode consisting of the v =0,1,2, ... ,6 harmonic oscillator functions.Our VCI calculations used the full space of all possible 1,2,…,6 modal excitations, and hence, they are designated FVCI.We have checked that the change in the VSCF and VCI ground state energies is negligible ͑i.e., Ͻ10 −7 a.u.͒ when the harmonic oscillator basis set is increased to v =8.

IV. SAMPLE CALCULATIONS
In this section we present an analysis of the calculated vibrational static and electro-optical ͑hyper͒polarizabilities of water and ammonia.Topics of discussion include the breakdown into contributions from P nr , P zpva , and P c−zpva , as well as the effect of anharmonic mode-mode coupling treated at different truncation levels in the normal mode expansion of the potential energy and property surfaces.The significance of the variational approach is that, for the first time, vibrational electro-optical properties are calculated beyond the lowest order of perturbation theory.From the perspective of this work water is particularly harmonic and serves as a first test to our new approach.Ammonia, on the other hand, possesses a very anharmonic umbrella mode, which is of particular interest for us and more challenging to the methodology.Note that our goal here is to validate the methodology rather than to calculate benchmark values for comparison with experiment.
The P zpva values in Tables I-III were computed using different truncations of the Taylor series expansions in terms of normal modes.These truncations are denoted as XmMtT, 16 where X may be the potential energy ͑V͒ or the electronic contribution to the static electrical property ͑ , ␣ , ␤ , ␥͒, m is the maximum number of coupled modes, and t is the highest power in the expansion ͑M and T are included for cosmetic reasons only͒.As explained in Sec.II A P c−zpva is obtained from finite field derivatives of P zpva .The c − zpva truncations given in Tables I-III are those used in evaluating the zpva from which they are computed.As usual P c−zpva contains terms induced by nuclear relaxation that involve higher-order derivatives than the P zpva from which it is obtained.In fact, the highest order will increase by the order of the field derivative performed in the c − zpva evaluation, e.g., if ␥ zzzz c−zvpa ͑0;0,0,0͒ is calculated from zpva , then the maximum order will increase by 3. Finally, when referring to our "best" P zpva and P c−zpva values we mean the ones computed using the most flexible property and potential energy treatments in the FVCI calculation.
It is important to distinguish between orders in the normal coordinate expansions and orders of perturbation theory ͑PT͒.For example, the 2M2T / V3M4T treatment in Table I includes all first-order PT contributions to zpva plus some ͑but not all͒ third-and higher-order PT contributions.Naturally, the subsequent c − zpva calculation will yield results that are different from those obtained by a KLB nuclear relaxation treatment 30 based on using all the zpva terms through a given order of PT.
Although water is quite harmonic, the vibrational contributions cannot be neglected if one wants to perform an accurate evaluation of the NLO properties.While the vibra-tional contribution to is only about 0.1%, it reaches 57% and 22% for the static first and second hyperpolarizabilities, respectively.Ordinarily the values of P nr and P c−zpva decrease when the number of static fields involved in the property decreases.In fact, when no static field is present the vibrational contribution to NLO properties vanishes within the IF approximation ͓␥͑− ; ,− , ͒ →ϱ is an exception because there is a cancellation of + and − that simulates a static field͔.Examples of this behavior may be found in Tables II and III.For our best c − zpva values, the ratio ͑␥ zzzz nr ͑0;0,0,0͒ + ␥ zzzz c−zvpa ͑0;0,0,0͒͒ / ␥ zzzz e ͑0;0,0,0͒ is 0.23, whereas ͑␥ zzzz nr ͑− ; ,0,0͒ →ϱ + ␥ zzzz c−zpva ͑− ; ,0,0͒ →ϱ ͒ / ␥ zzzz e ͑0;0,0,0͒ is only 0.04 and ͑␥ zzzz nr ͑−2 ; , ,0͒ →ϱ + ␥ zzzz c−zpva ͑−2 ; , ,0͒ →ϱ ͒/␥ zzzz e ͑0;0,0,0͒ is further reduced to 0.01.The same rule is also followed by the static and dc-P first hyperpolarizabilities although ␤ zzz c−zpva ͑− ; ,0͒ →ϱ is, surprisingly, larger than its static counterpart.This unusual behavior of the first hyperpolarizability can probably be attributed to the fact that its magnitude is exceptionally small for a strongly dipolar molecule.
As an indicator of convergence in the magnitude of the vibrational contributions with respect to the order of anharmonicity we have, in the past, 18 used the ratios P zpva / P e and P c−zpva / P nr .Using our best P zpva and P c−zpva values for water we find that P zpva / P e Ͻ 0.08 and P c−zpva / P nr Ͻ 0.21 in all cases.Since these ratios are much less than unity, we con- clude that the vibrational anharmonicity terms are converging rapidly, in agreement with our previous finding for the static properties. 18e now proceed to analyze the convergence with respect to the maximum number of coupled modes and the truncation level of the Taylor series expansion for P zpva .Consider, first, the convergence behavior of the property for a fixed V3M4T potential energy expansion.The dipole moment, linear polarizability, and second hyperpolarizability calculated at the P2M2T level are very close to those obtained using our highest level treatment.For the first hyperpolarizability the behavior is less clear, but again this could be related to the very small magnitude of this property in comparison with the large dipole moment.One might wonder whether the differences between P1M1T and P2M2T, which are substantial, are due to the truncation of the mode-mode coupling or the property expansion or both.We have found that the P1M2T results ͑P2M1T is impossible͒ are identical to P2M2T.This indicates the nonessential role of two-mode coupling terms.In fact, the observed behavior can be rationalized on the basis of PT.In the PT formulas for P zpva ͑see, for instance, Eqs.͑11͒-͑13͒ of Ref. 30͒ there are two first-order terms.One involves only mechanical anharmonicity and, thereby, requires only first derivatives of the property ͑P1M1T͒.The other requires second derivatives but is limited to diagonal terms, i.e., P1M2T.Thus, in first-order PT, P1M2T is necessary and sufficient to treat electrical anharmonicity.
Next, we discuss the convergence with respect to the potential energy surface for a given ͑in particular, the best͒ property expansion.The V3M4T potential energy expansion was improved by adding two-mode fifth-and sixth-order terms ͑V3M4T + V2M6T͒ except for ‫ץ‬ 6 V / ‫ץ‬Q i 3 ‫ץ‬Q j 3 .The threemode ‫ץ‬ 5 V / ‫ץ‬Q i 3 ‫ץ‬Q j ‫ץ‬Q k terms were included as well since they could be obtained without any extra computational cost.Except, as usual, for ␤ zzz c−zpva ͑0;0,0͒ the changes in P zpva and P c−zpva compared to V3M4T are negligible ͑Ͻ0.3% of the total property͒.Although not conclusive, this does indicate that the values are converged at the V3M4T truncation level.
In agreement with Kongsted and Christiansen, 17 the VSCF and FVCI static zpva contributions to the four electrical properties are very similar to one another.For c − zpva the values of ␣ and ␥ are also virtually identical.At this point it should not be surprising that the agreement for ␤ is not as ␥ zzzz nr ͑−2 ; ,0,0͒ →ϱ −11 The 3M5T derivatives ‫ץ‬ 5 V / ‫ץ‬Q i 3 ‫ץ‬Q j ‫ץ‬Q k are also included. b The 2M6T derivatives ‫ץ‬ 6 V / ‫ץ‬Q i 3 ‫ץ‬Q j 3 are not included.
The 2M4T derivatives ‫ץ‬ 4 e / ‫ץ‬Q i 3 ‫ץ‬Q j are also included. d The 2M5T derivatives ‫ץ‬ 5 e / ‫ץ‬Q i 4 ‫ץ‬Q j are also included.

B. Ammonia
In Tables IV-VI we present the vibrational contributions to the dipole moment and to the diagonal component of the ͑hyper͒polarizabilities along the C 3 symmetry axis of ammonia.The corresponding pure electronic contributions are also presented, as well as the vibrational contributions obtained with the umbrella normal mode frozen.
Vibrations play a more important role in the NLO properties of ammonia than water.For the static first and second hyperpolarizabilities our best calculations yield vibrational contributions that are, respectively, 5.7 and 2.5 times larger than the corresponding pure electronic value.These contributions will certainly be important at frequencies in the infrared range. 31On the other hand, the vibrational electrooptical hyperpolarizabilities are reduced by at least an order of magnitude in the infinite optical frequency limit.
The data presented in Tables IV-VI show that the umbrella mode has a major effect in most cases.For instance, our best values for the vibrational static second hyperpolar-izability with and without the umbrella mode are 17 454 and 351 a.u., respectively.The dominance of the umbrella mode is not always as extreme for the other vibrational properties, although the reduction in magnitude due to freezing this mode remains greater than 76% ͑except for the dipole moment and ESHG͒.As far as the various vibrational terms listed in the tables are concerned, nuclear relaxation ͑nr͒ is usually the most important and, in that case, the other 3N − 7 modes never contribute more than 14% ͑4% for static properties͒ of the total.The contribution of the c − zpva and ͑usually smaller͒ zpva terms is more variable.For the static c − zpva, in particular, the role of the umbrella mode is once more overwhelming ͑i.e., always larger than 98%͒.
Again the PT treatment can help analyze the above results.When the number of static fields associated with the property decreases, higher-order terms are eliminated in the formula for P nr ͑or P c−zpva ͒.For example, while the PT expression for ␤ zzz nr ͑0;0,0͒ ͑two static fields͒ contains zerothand first-order terms, the expression for ␤ zzz nr ͑− ; ,0͒ →ϱ ͑one static field͒ depends only on the zeroth-order term.Thus, the reduction in the role of the umbrella mode in going from static to electro-optical P nr ͑and P c−zpva ͒ confirms that the importance of this mode increases with the order of anharmonicity included, as might have been expected.
For our best NH 3 values, the convergence indicator The 3M5T derivatives ‫ץ‬ 5 V / ‫ץ‬Q i 3 ‫ץ‬Q j ‫ץ‬Q k are also included.
The 2M6T derivatives ‫ץ‬ 6 V / ‫ץ‬Q i 3 ‫ץ‬Q j 3 are not included.
P zpva / P e is always less than 0.03.However, P c−zpva / P nr is much larger and, in the case of the static ␥, is greater than unity.When the umbrella mode is frozen the latter ratio is reduced to 0.33.This tells us that, for the static ␥, the vibrational anharmonicity of the umbrella mode must be treated more accurately by a method not involving an expansion in normal modes.
In order to investigate the convergence behavior with regard to the property and potential energy expansions in the FVCI approach we consider, first, the various property expansions for a fixed V4M4T potential energy expansion.For P zpva the differences between P2M2T and our most flexible property expansion are negligible.For P c−zpva the differences increase, due entirely to the umbrella mode, but are still never larger than 2.7% of the total property value.Thus, we conclude that the property expansions are converged at the P2M2T truncation level.
Our next step is to compare the V4M4T and V4M4T + V2M6T potential energy expansions, using in both cases our best property expansion.Now the differences are much larger, reaching 26.0% and 63.3% of the total property value ͑including P e ͒ for the static first and second hyperpolarizabilities, respectively.In these two cases convergence has not been achieved.Again, it is obvious by inspection that the problem is due to the umbrella mode and the solution lies in a treatment of this mode not involving a power series expansion of the potential energy associated with this vibration.
The poor convergence is confined to the static hyperpolarizabilities.For the dipole moment, static polarizability, and infinite frequency electro-optical hyperpolarizabilites the differences are always lower than 2.2% ͑even without taking into account the fact that the frequency-dependent P e will be larger than the static value͒.Not surprisingly, when the umbrella mode is frozen the difference is much smaller yetless than 0.03% ͑including static hyperpolarizabilities͒.
The effect of vibrational correlation varies widely depending on the particular contribution and property.It is negligible ͑Ͻ0.3% of the VSCF total property value͒ for the zpva contribution to all static properties.However, vibrational correlation becomes more important for the c − zpva contribution, in which case it increases the magnitude of the VSCF total property value ͑best calculation͒ by anywhere from 0.7% to 17.3% ͑excluding ESHG which is negligibly affected͒.Finally, and most importantly, when the umbrella mode is deleted, the difference between VSCF and FVCI is negligible.

V. CONCLUSIONS
It is now well established that molecular vibrations often make an important, or even major, contribution to the total The 3M5T derivatives ‫ץ‬ 5 V / ‫ץ‬Q i 3 ‫ץ‬Q j ‫ץ‬Q k are also included. b The 2M6T derivatives ‫ץ‬ 6 V / ‫ץ‬Q i 3 ‫ץ‬Q j 3 are not included.
The 2M4T derivatives ‫ץ‬ 4 e / ‫ץ‬Q i 3 ‫ץ‬Q j are also included.
The 2M5T derivatives ‫ץ‬ 5 e / ‫ץ‬Q i 4 ‫ץ‬Q j are also included.
͑hyper͒polarizability.This vibrational contribution can depend significantly upon the anharmonic, as well as the harmonic, parameters that govern the shape of the potential energy and electric dipole property surfaces.Certain low-order ͑in terms of perturbation theory͒ anharmonic terms are included in a so-called nr treatment.This is sufficient in some cases.In general, however, it is necessary to account for higher-order contributions.
The entire remaining effect of anharmonicity ͑after nr͒ on static vibrational ͑hyper͒polarizabilities can be obtained, in principle, by calculating the zero-point vibrational energy as a function of a uniform external static electric field.This quantity can be computed in lowest order and, with some difficulty, by perturbation theory.However, for molecules with strongly anharmonic modes, the lowest-order approximation is inadequate.Previously, it has been found more satisfactory to use vibrational CI together with a normal coordinate expansion of the potential energy surface.In principle, an analogous procedure can be followed to calculate dynamic vibrational electro-optical ͑hyper͒polarizabilities in the infinite optical frequency approximation.For the latter it is necessary to obtain field-dependent zero-point vibrationally averaged ͑hyper͒polarizabilities as well as vibrational potential energies.In this paper we have shown the feasibility of carrying out such calculations, at least for small molecules ͑at present͒, using essentially a full VCI in the direct method of Ref. 17.
Two specific examples were considered, namely, H 2 O and NH 3 .These calculations were limited to diagonal ͑hyper͒polarizabilities along the symmetry axis.The former molecule provides a test of the methodology for weak anharmonicity, whereas the latter contains the strongly anharmonic umbrella mode.For comparison purposes, we also obtained the corresponding static vibrational and pure electronic properties.The effect of various truncations of the potential energy and property surfaces were examined.In the case of water it was found ͑with one exception͒ that, as far as the property surface is concerned, the effect of mode-mode coupling is negligible and it is sufficient to terminate the normal coordinate expansion at second order.On the other hand, for The 3M5T derivatives ‫ץ‬ 5 V / ‫ץ‬Q i 3 ‫ץ‬Q j ‫ץ‬Q k are also included.
The 2M6T derivatives ‫ץ‬ 6 V / ‫ץ‬Q i 3 ‫ץ‬Q j 3 are not included.
The 2M4T derivatives ‫ץ‬ 4 e / ‫ץ‬Q i 3 ‫ץ‬Q j are also included.
The 2M5T derivatives ‫ץ‬ 5 e / ‫ץ‬Q i 4 ‫ץ‬Q j are also included.
the potential energy surface a fourth-order normal coordinate expansion is required.The one exception to these observations is the static first hyperpolarizability, where large percentage changes are found because the property value is unusually small.In the case of ammonia the vibrational contributions are much more important.Calculations done with the umbrella mode frozen indicate that this mode plays the primary role except in some instances where the vibrational contribution to the property is very small.The umbrella vibration does couple strongly to other modes so that one needs to include two-mode coupling in the property surface, although second order in the normal mode expansion remains sufficient.A significant difference from water is that the potential energy surface obtained using a fourth-order normal mode expansion and up to four coupled modes is no longer adequate for the static hyperpolarizabilities.On the other hand, this level of truncation is satisfactory for the electro-optical properties.
There remains a need for extension of the computational methodology presented here in a number of directions including separate numerical treatment of highly anharmonic modes and improved efficiency for application to larger molecules with more vibrational degrees of freedom.Efforts along both of these lines are currently in progress.

TABLE I .
Electronic and vibrational contributions to and to the diagonal component of ␣ in the direction of , calculated for H 2 O by VSCF and FVCI.Results are reported for different truncation levels, with respect to the number of coupled modes and the order of the Taylor series expansion in normal modes ͑see text for notation͒, for both the potential and the relevant property.All quantities are in a.u.

TABLE II .
Electronic and vibrational contributions to the diagonal component of ␤ in the direction of calculated for H 2 O by VSCF and FVCI.Results are reported for different truncation levels, with respect to the number of coupled modes and the order of the Taylor series expansion in normal modes ͑see text for notation͒, for both the potential and the relevant property.All quantities are in a.u.

TABLE III .
Electronic and vibrational contributions to the diagonal component of ␥ in the direction of , calculated for H 2 O by VSCF and FVCI.Results are reported for different truncation levels, with respect to the number of couple modes and the order of the Taylor series expansion in normal modes ͑see text for notation͒, for both the potential and the relevant property.All quantities are in a.u.

TABLE IV .
Electronic and vibrational contributions to and to the diagonal component of ␣ in the direction of , calculated for NH 3 by VSCF and FVCI.Results are reported for different truncation levels, with respect to the number of coupled modes and the order of the Taylor series expansion in normal modes ͑see text for notation͒, for both the potential and the relevant property.Results with ͑3N −6͒ and without ͑3N −7͒ umbrella mode degree of freedom are present.All quantities are in a.u.
bThe 2M5T derivatives ‫ץ‬ 5 e / ‫ץ‬Q i 4 ‫ץ‬Q j are also included.c

TABLE V .
Electronic and vibrational contributions to the diagonal component of ␤ in the direction of calculated for NH 3 by VSCF and FVCI.Results are reported for different truncation levels, with respect to the number of coupled modes and the order of the Taylor series expansion in normal modes ͑see text for notation͒, for both the potential and the relevant property.Results with ͑3N −6͒ and without ͑3N −7͒ umbrella mode degree of freedom are present.All quantities are in a.u.

TABLE VI .
Electronic and vibrational contributions to the diagonal component of ␥ in the direction of calculated for NH 3 by VSCF and FVCI.Results are reported for different truncation levels, with respect to the number of coupled modes and the order of the Taylor series expansion in normal modes ͑see text for notation͒, for both the potential and the relevant property.Results with ͑3N −6͒ and without ͑3N −7͒ umbrella mode degree of freedom are present.All quantities are in a.u.