Microscopic and macroscopic polarization within a combined quantum mechanics and molecular mechanics model

A polarizable quantum mechanics and molecular mechanics model has been extended to account for the difference between the macroscopic electric field and the actual electric field felt by the solute molecule. This enables the calculation of effective microscopic properties which can be related to macroscopic susceptibilities directly comparable with experimental results. By seperating the discrete local field into two distinct contribution we define two different microscopic properties, the so-called solute and effective properties. The solute properties account for the pure solvent effects, i.e., effects even when the macroscopic electric field is zero, and the effective properties account for both the pure solvent effects and the effect from the induced dipoles in the solvent due to the macroscopic electric field. We present results for the linear and nonlinear polarizabilities of water and acetonitrile both in the gas phase and in the liquid phase. For all the properties we find that the pure solvent effect increases the properties whereas the induced electric field decreases the properties. Furthermore, we present results for the refractive index, third-harmonic generation ~THG!, and electric field induced second-harmonic generation ~EFISH! for liquid water and acetonitrile. We find in general good agreement between the calculated and experimental results for the refractive index and the THG susceptibility. For the EFISH susceptibility, however, the difference between experiment and theory is larger since the orientational effect arising from the static electric field is not accurately described. © 2005 American Institute of Physics. @DOI: 10.1063/1.1831271 #


I. INTRODUCTION
4][5][6][7] The main reason for this is that TD-DFT provides a level of accuracy which in most cases is sufficient at a lower computational requirement than other methods.[10][11][12][13][14][15][16][17] In recent years applications of TD-DFT to calculate properties of molecules in solution has also been presented. 18 -27Among the methods for treating solvent effect on molecules are the combined quantum mechanical and molecular mechanics ͑QM/MM͒ models. 28 -39The QM/MM model has recently been extended to study molecular response properties in solution within TD-DFT. 18 -20,22-24An example of such a QM/MM method is the discrete solvent reaction field ͑DRF͒ model which we have recently developed. 18 -20The model combines TD-DFT ͑QM͒ description of the solute molecule with a classical ͑MM͒ description of the discrete solvent molecules.The latter are represented using distributed atomic charges and polarizabilities.
An important feature of the model is the inclusion of polarizabilities in the MM part which allows for the solvent molecules to be polarized by the solute and by interactions with other solvent molecules.Van Duijnen et al. 50applied basically the same polarizable QM/MM formalism, implemented in a conventional wave function package, to calculate the static ͑hyper͒polarizability of acetone in ten different solvent.Also Kongsted et al. 41 stressed the importance of inclusion of polarizabilities in the calculations of response properties.Another important feature of the model is the inclusion of a short range damping of the interactions.This has been included in two ways.The first is the use of the modi-fied dipole interaction model of Thole 42 which avoids the''polarization catastrophe'' by introducing smeared out dipoles which mimic the overlapping of the charge distributions at short distances.The second way is that at short distances between the QM and the MM part the QM/MM interactions are damped to account for the short range repulsion in an approximate way.This is done by replacing the point charge by a Gaussian charge distribution with a unit width and point dipoles smeared out in a similar manner. 18,19,43n this work we will extend the QM/MM formalism to also include the so-called local field factors, i.e., the difference between the macroscopic electric field and the actual electric field felt by the solute.This will enable the calculation of effective microscopic properties which can be related to the macroscopic susceptibilities.The macroscopic susceptibilities can then be compared directly with experimental results.There exist in the literature some other approaches to calculate effective properties and relate these to the macroscopic properties, 40,44 -50 which differ from this work in the way the solvent was represented.First we will describe the theoretical framework and then we will present numerical results for liquid water and liquid acetonitrile.These liquids are chosen since there exist several theoretical and experimental studies of the microscopic and macroscopic properties.

II. THEORY
The basic concept of nonlinear optics is the expansion of the total macroscopic polarization in a material in powers of the macroscopic electric field where the expansion coefficients define the macroscopic ͑nonlinear͒ susceptibilities.Similarly, the total microscopic polarization ͑dipole moment͒ is expanded in terms of the total microscopic electric field with expansion coefficients defining the microscopic ͑nonlin-ear͒ polarizabilities.However, in the literature several different conventions exist for describing nonlinear optical properties 51 which differ in the numerical coefficients used.Therefore, in order to compare values obtained in different conventions it is important to correct for the differences in the numerical factors used.This has been clarified by Willetts et al. 51 but remains a problem, because often it is not stated explicitly which convention is used.In this work we will use a perturbation series expansion for the macroscopic polarization, which is often used for experimental properties, and a Taylor series expansion for the microscopic polarization, which frequently is used for theoretical properties.

A. The macroscopic polarization
The macroscopic polarization of a material in the presence of a macroscopic electric field F mac is expressed as a power series in the field strength as 52,53 where P 0 is the permanent polarization, (1) the linear optical susceptibility, (2) the second-order nonlinear optical susceptibility, and (3) the third-order nonlinear optical sus-ceptibility.The subscripts I,J,K,L, . . .denotes space-fixed axes and the Einstein summation convention is used for repeated subscripts.If we consider the macroscopic field to be a superposition of a static and an optical component, the macroscopic polarization can be expressed as 52,53 The Fourier amplitudes of the polarization are then given in terms of the frequency-dependent susceptibilities as 52,53 P I s ϭ␦ s ,0

͑4͒
where the output frequency is given as the sum of input frequencies s ϭ ͚ a a .The numerical coefficients K (Ϫ s ; a ,...) arise from the Fourier expansion of the electric field and polarization and ensures that all susceptibilities of the same order have the same static limit.A tabulation of the coefficients can be found in Ref. 54 and 55.The frequency-dependent susceptibilities can then be found from Eq. ͑4͒ by differentiation which gives the linear susceptibility the second-order nonlinear susceptibility and the third-order nonlinear susceptibility .

B. The microscopic polarization
Similarly the microscopic polarization ͑dipole moment͒ can be expanded in terms oscillating at different frequencies as 44,51 The microscopic dipole moment is then usually given by a Taylor expansion as 44,51

͑9͒
where ␣ 0 is the permanent electric dipole moment, ␣ ␣␤ (Ϫ s ; s ) is the polarizabillity, ␤ ␣␤␥ (Ϫ s ; a , b ) is the first hyperpolarizability, and, ␥ ␣␤␥␦ (Ϫ s ; a , b , c ) is the second hyperpolarizability.The numerical coefficients K(Ϫ s ; a ,...) are the same as for the macroscopic polarization and again this ensures that all ͑hyper͒polarizabilities of the same order have the same static limit.The subscripts ␣,␤,␥,... denote molecule-fixed axes and again the Einstein summation convention is used for repeated subscripts.The microscopic polarization is expanded in terms of the actual total electric field F b ,␥ tot felt by the molecule.In the condensed phase the actual electric field felt by the molecule is different from the macroscopic electric field.Therefore, in order to express the macroscopic properties in terms of the microscopic properties we need to relate the actual electric field at a molecule to the macroscopic electric field.

C. The local electric field
The concept of relating the actual electric field, often called the internal or local field, to the macroscopic field dates back to the work of Lorentz. 56,57Lorentz 56 derived a simple relation between the internal electric field, the macroscopic electric field, and the macroscopic polarization of the system, and due to its simplicity Lorentz local field theory is still used. 52,53,57,58The central idea is that only close to the molecule we need to consider explicitly the field from nearby molecules, so the total system is separated into a macroscopic region far from the molecule and a microscopic region close to the molecule.The molecules in the macroscopic region can then be described by the average macroscopic properties.Therefore, inside a macroscopically small, but microscopically large, virtual cavity V we subtract the contribution from the macroscopic electric field and replace it by the correct discrete local field, where F s ,␣ pol is the macroscopic electric field in the cavity V and F s ,␣ disc (⍀) is the discrete electric field in the cavity V which depends on the local configuration ⍀ of the molecules inside the cavity.Since we are not allowing the macroscopic region to adjust to the presence of the cavity the polarization remains homogeneous. 57This approach neglects that a static electric field tends to orient molecules with a permanent dipole 57,59 and therefore a correction due to Onsager 59 is often used for a static electric fields.Lorentz showed 56 that for a cubic arrangement of identical particles the discrete field was zero.This is also true on average for a completely random distribution where there is no correlation between the induced dipoles and the position of the molecules. 57For a spherical cavity the macroscopic field is simply given in terms of the macroscopic polarization 56,57 and the total electric field can be written as where we have split the discrete electric field F disc into two different contributions, F ind and F perm .The first term arises from the interactions of the macroscopic electric field with the other molecules in the cavity, i.e., accounts for the induced polarization of the surrounding molecules due to the electric field.The second term accounts for the interactions between the molecules when there is no electric field present, i.e., arises from the permanent charge distribution of the surrounding molecules.However, depending on the theoretical model used for describing the microscopic region, this splitting of the discrete electric field is not always possible nor necessary.The last two terms depend strongly on the local configuration of the molecules in the cavity and are inherently microscopic in nature and it is therefore better to treat these fields explicitly within the microscopic model used.

D. The effective molecular properties
Instead of expanding the induced dipole moment in terms of the total field, Eq. ͑9͒, we expand it in terms of an effective macroscopic electric field This expansion in terms of the effective field defines the so-called effective properties 44 as

͑13͒
These effective properties give an induced dipole moment due to the effective macroscopic electric field which is identical to the induced dipole moment in Eq. ͑9͒ and are the properties which we will relate to the experimental susceptibilities.This means that the microscopic contributions to the total field are incorporated into the effective properties.These effective properties could be compared with experi-mental results corrected for differences between the total field and the macroscopic electric field by the Lorentz/ Onsager local field method. 44

E. The solute molecular properties
Since we have separated the discrete field into the two contributions mentioned above we can also choose to expand the induced dipole moment in terms of the field arising directly from the macroscopic electric field, where the field arising from the interactions between the molecules when there is no macroscopic field is incorporated into the properties.This gives an expansion which defines the so-called solute properties 44 as

͑15͒
These solute properties relate to the macroscopic properties corrected for the field from the dipoles of all other molecules induced by the macroscopic field in addition to the Lorentz local field.This corresponds to a thought experiment where the macroscopic field is allowed to propagate inside the cavity without being modified by interactions with the molecules.

F. Relating the macroscopic and the microscopic polarization
The macroscopic polarization is related to the average microscopic dipole moment per molecule by 52,58

͑16͒
where N d is the number density and the brackets, ͗͘, denote orientational averaging, and relate the molecule-fixed axes to the space-fixed axes. 58,60Inserting the expansion of the dipole moment in terms of the effective macroscopic field, Eq. ͑13͒, we can express the macroscopic polarization in terms of the effective ͑hyper͒polarizabilities as,

͑17͒
We see that the averaging is done on the product of the ͑hyper͒polarizabilities and the effective fields.This is exactly the reason why the total electric field was split into an effec-tive macroscopic part and a microscopic part which was incorporated into the ͑hyper͒polarizabilities in Eq. ͑13͒ by expanding the dipole moment in terms of the effective field.
Since the effective field is macroscopic we can take it outside the averaging and express the macroscopic polarization in terms of orientational averages of the effective ͑hyper͒polarizabilities as

G. Local field factors
Comparing Eq. ͑18͒ with Eq. ͑4͒ we see that we have to consider derivatives of the effective electric field with respect to the macroscopic electric field.This is usually done by introducing the so-called local field factors which relate the macroscopic field to the effective field.
Using the definition of the effective electric field, Eq. ͑12͒, we obtain a local field factor where ⑀ (1) ( a ) is the optical dielectric constant at frequency a .This is the Lorentz form of the local field factors.However, as mentioned above this does not account for the fact that for a static macroscopic field the molecules tend to orient.Onsager 59 analyzed this problem and suggested the following form for the local field factor, to be used for static electric fields, where ⑀ (1) (0) is the dielectric constant and n (1) (0) is the refractive index at zero frequency.We see that the Onsager field factor reduces to the Lorentz factor for optical fields by using the relation n (1) ( a ) 2 ϭ⑀ (1) ( a ).

H. Orientational averaging
In order to relate the molecule-fixed axes to the spacefixed axes we need also to consider molecular rotations ͑or orientations͒.The orientational averaging and thermal averaging of the dipole moment will be done using classical theory and is given by 60 where ⌽ ␣I is the cosine of the angle between the molecular axis ␣ and the laboratory axis I.The angular dependent part of the energy in the presence of the electric field is given by We note that it is only the solute properties which are responsible for the change in the energy due to the electric field.If we expand the exponential and only terms of the order (kT) Ϫ1 are retained we get By combining the definitions of the susceptibilities in Eqs.͑5͒, ͑6͒, and ͑7͒ with the expression for the macroscopic polarization in terms of the effective ͑hyper͒polarizabilities we can obtain a link between the macroscopic and the microscopic properties.

I. Refractive index
The macroscopic quantity determining the refractive index is the linear optical polarization due to an optical electric field.By inserting the definition of the effective field, Eq. ͑12͒, into Eq.18 and using the definition of the linear susceptibilitiy, Eq. ͑5͒, we obtain

͑25͒
It is noted that there is no contribution from the rotation of the dipole moment since the optical field is considered to be oscillating faster than the permanent dipole moments can be oriented.The linear susceptibility can then be written in terms of the mean effective polarizability by rewriting Eq. ͑25͒ as which is the standard expression for the susceptibility corrected for the Lorentz local field, 52,61 although using the effective rather than the gas phase polarizability.The susceptibility is related to the refractive index or the optical dielectric constant of the system as which is the familiar Lorentz-Lorenz or Clausius-Mossotti equation, 56,57,61 again with the effective polarizability instead of the gas phase polarizability.

J. Dielectric constant
In a dielectric constant measurement the polarization due to a static electric field is measured and the corresponding susceptibility is given by The first term is the rotational contribution arising from the permanent dipole moment and is given by The second term is the isotropic orientation average of the polarizability, often referred to as the mean polarizability, given by 60 and denoted ␣ ¯.By combining these two terms the equation for the linear susceptibility can be rewritten as

͑31͒
The dielectric constant is then related to the linear susceptibility through the usual equation ⑀ (1) ͑ 0 ͒ϭ1ϩ4 ZZ (1) ͑ 0;0͒.͑32͒ It should be noted that the susceptibility in Eq. ͑31͒ depends on the dielectric constant through the Onsager local field factor, L 0 , and is therefore not a defining equation.However, this will allow us to make an estimate of the dielectric constant.

K. Third-harmonics generation
The first nonlinear susceptibility we will consider is the third-order nonlinear susceptibility which arises from three optical electric fields corresponding to the THG experiments.The THG susceptibility is then obtained by inserting Eq. ͑12͒ into Eq.͑18͒ and using the definition of the susceptibilitiy, Eq. ͑7͒.This gives the susceptibility as where we see that we have a contribution both from the linear susceptibility and from the third-order nonlinear susceptibility.The isotropic orientation average of the second hyperpolarizability, often referred to as the mean or parallel second hyperpolarizability, is given by 60

͑34͒
We can then rewrite the THG susceptibility as where we have used the relation between the linear susceptibility and the dielectric constant in Eq. ͑27͒.Using Eq. ͑26͒ we can express the term with the effective polarizability in terms of the dielectric constant at 3. This allows us to express the third-order nonlinear susceptibility as which is the form for the nonlinear susceptibility well known from standard Lorentz local field theory with nϩ1 local field corrections, where n is the number of applied fields.

L. Electric field induced second-harmonic generation
The second nonlinear susceptibility we will consider is the third-order nonlinear susceptibility arising from two optical and one static electric field and corresponds to the electric field induced second-harmonic generation ͑EFISH͒ experiments.The EFISH susceptibility is then obtained by inserting Eq. ͑12͒ into Eq.͑18͒ and using the definition of the susceptibilitiy, Eq. ͑7͒.This gives the susceptibility as We see that the EFISH susceptibility consists of three terms and compared with the THG expression there is a term depending on the effective first hyperpolarizability.The second term is a rotational contribution analogous with the dielectric constant and is

͑38͒
where the mean hyperpolarizability ␤ ¯ʈ in the direction of the dipole moment, here the z axis, is introduced This allows us to rewrite the EFISH susceptibility as

M. The discrete solvent reaction field model
In the QM/MM method the solvent molecules ͑MM͒ are treated with a classical force field and the interactions between the solute and solvent are described with an effective operator.In the QM/MM method the total ͑effective͒ Hamiltonian for the system is written as 28 -39 H ˆϭH ˆQM ϩH ˆQM/MM ϩH ˆMM , ͑41͒ where H ˆQM is the quantum mechanical Hamiltonian for the solute, H ˆQM/MM describes the interactions between solute and solvent, and H ˆMM describes the solvent-solvent interactions.
We have recently developed such a method for studying solvent effect on molecular properties which we denoted the DRF ͑Refs.18 -20͒ where the QM part is treated using DFT.
Within the discrete solvent reaction field model the QM/MM operator at a point r i is given by where the first term el is the electrostatic operator and describes the Coulombic interaction between the QM system and the permanent charge distribution of the solvent molecules.The second term pol , is the polarization operator and describes the many-body polarization of the solvent molecules, i.e., the change in the charge distribution of the solvent molecules due to interaction with the QM part and other solvent molecules.The charge distribution of the solvent is represented by atomic point charges and the many-body polarization term is represented by induced atomic dipoles at the solvent molecules.
For a collection of atomic polarizabilities in an electric field, assuming linear response, the induced atomic dipole at site s is given by where ␣ a,␣␤ is a component of the atomic polarizability tensor at site s, which for an isotropic atom gives ␣ s,␣␤ ϭ␦ ␣␤ ␣ s , and is the screened dipole interaction tensor. 42,62,18Here we neglect the frequency dependence of the classical part, i.e., the atomic polarizability is frequency independent, but the model can easily be extended to include also this effect. 43,63 s,␤ init () is the initial electric field at site s and is in this work extended to also include the macroscopic electric field.The initial field then consist of four terms where F t,␤ QM,el () is the field arising from the frequencydependent electronic charge distribution of the QM part, F t,␤ QM,nuc the field from the QM nuclei, F t,␤ M M,q the field from the point charges at the solvent molecules, and F t,␤ mac () the macroscopic electric field.The inclusion of the macroscopic electric field in Eq. ͑44͒ describes the induced dipole moments in the solvent due to macroscopic electric field.This was the reason for splitting the discrete electric field in Eq. ͑11͒ into a part induced by the macroscopic electric field and a part existing even without a macroscopic field.Therefore, if the macroscopic electric field is included in Eq. ͑44͒ we will be calculating the effective properties and if it is excluded we calculate the solute properties.
We can now calculate and distinguish between different effects in going from microscopic properties to macroscopic properties.Since the permanent discrete electric field in Eq. ͑11͒ is always present we will associate this with a pure solvent effect, i.e., the solute properties include this solvent effect.The effective properties then includes the effects of a microscopic induced field in the solvent due to the macroscopic electric field.Finally, by combining the effective properties with the macroscopic local field factors ͑Lorentz/ Osager͒ we can obtain the macroscopic susceptibilities.

III. COMPUTATION DETAILS
7][68] The details of the implementation are described in Refs.18 -20 and is in this work only extended to include the macroscopic electric field in Eq. ͑44͒.
The basis set used for water consists of a large eventempered basis set of Slater-type orbitals with orbital exponent ϭ␣␤ i , iϭ1,...,n ͑details given in Ref. 19͒.For acetonitrile the standard TZP basis set was augmented with firstorder and second-order field induced polarization functions taken from Ref. 69.All the calculations were done using the BP-GRAC ͑gradient-regulated asymptotic connection BP͒ potentials. 13,14The BP-GRAC potential sets the highest ordered molecular orbital ͑HOMO͒ level at the first ionization potential ͑IP͒ and therefore requires the IP as input.For water IPϭ0.45 a.u. and for acetonitrile IPϭ0.46143 a.u. was used.The values for IP were obtained from calculations using the SAOP potential for which it has been shown that the HOMO level corresponds well with the experimental IP. 70he molecular dynamics ͑MD͒ simulations were performed with the discrete reaction field polarizable force field 35,62,[71][72][73][74] using the DRF90 program. 71For both water and acetonitrile a MD-simulation of 50 ps with a timestep of 1 fs was performed at 298K, using a Nose-Hoover thermostat 75 ͑with ϭ1 ps) to keep the temperature constant and a soft wall force potential 71 to keep the particles inside the simulation box.After every 0.5 ps, the configuration of solvent molecules was kept and the QM/MM calculations were performed.In the simulation of water, 256 molecules were placed in a spherical box of 23.12 bohrs; for acetonitrile, 128 molecules were placed in a spherical box of 26.15 bohrs.The sizes of the simulation boxes were chosen so that the simulated macroscopic densities correspond to the experimental values of 0.998 and 0.786 kg/l, respectively.In the simulations, the molecules were treated as rigid bodies using quaternions. 76Standard atomic polarizabilities 62 for all atoms were used in the simulations.For acetonitrile MDC-d charges were used 77 that were obtained from DFT calculations in a TZ2P basis set with the ADF program, while for water charges were fitted to reproduce the experimental dipole moment.For the van der Waals interaction a 6-12 Lennard-Jones potential were used for water and the standard DRF90 potential for acetonitrile.For water the Lennard-Jones parameters were taken from Ref. 78 and adjusted to match the point charges and atomic polarizabilities used in this work.The new parameters obtained are Rϭ1.7385Å and ⑀ϭ0.2900 kcal/mol located on the oxygen atom.By inspecting radial distribution functions, it was checked that the solvent shells around the central molecule were correctly represented.
The atomic parameters, i.e., point charges and atomic polarizabilities, needed for the solvent molecules in the QM/MM calculations are: For water the point charges are q H ϭ0.3295 and q O ϭϪ0.6590 a.u. and the atomic polarizabilities are ␣ H ϭ0.0690 and ␣ O ϭ9.3005 a.u.For acetonitrile the point charges are q C1 ϭ0.288 340 a.u, q C2 ϭϪ0.009 643 a.u, q H ϭ0.017 028 and q N ϭϪ0.329781, where C2 is the carbon atom attached to the nitrogen.The atomic polarizabilities are ␣ C ϭ8.6959, ␣ H ϭ2.8382, and ␣ N ϭ3.5042 a.u.

IV. RESULTS
In the following we will present microscopic and macroscopic properties for the two liquids water and acetonitrile.The solute and the effective properties will be presented as averaged over the 101 different solvent configurations.The standard deviation will also be displayed to indicate the average fluctuation in the properties due to the different solvent configuration.All microscopic properties will be given in atomic units ͑a.u.͒ whereas the macroscopic susceptibilities will be presented in cgs units ͑esu͒.

A. Gas phase results
In Table I we present DFT results for , ␣ ¯(Ϫ;), ␤ ¯ʈ(Ϫ2;,), and ␥ ¯ʈ(Ϫ2;,,0) for water and acetonitrile in the gas phase.The results have been calculated at ϭ0.0428 a.u.(ϭ1064 nm) for water and ϭ0.0885 ( ϭ514.5 nm) for acetronitrile and are compared both with experimental and ab initio coupled cluster single doubles ͑CCSD͒ results.In general we find good agreement between the DFT results and the CCSD results for all properties.The largest difference of ϳ25% between the calculated values is in ␥ ¯ʈ and ␤ ¯for acetonitrile.The DFT results for ␤ ¯ʈ is lower than the CCSD results whereas for ␥ ¯ʈ the opposite is found.If we compare with results obtained from EFISH experiments we see that for water there is an excellent agreement between the calculated and the experimental results for all properties.In the case of acetonitrile we see that for ␤ ¯ʈ the DFT is in better agreement with the experiment whereas for ␥ ¯it is the CCSD results.However, as mentioned in the theory section the measured quantity in the EFISH experiment is ⌫ ¯ʈϭ z ␤ ¯ʈ/3k b T ϩ␥ ¯ʈ .Therefore, this value is also reported for the different methods in Table I.Again, we see that there is good agreement between theory and experiment for water.For acetonitrile the DFT and CCSD results are within 10% and 20%, respectively, of the experimental results.

B. Microscopic response properties
In Table II we present the dipole moment of acetonitrile and water in the gas phase and in the liquid phase.For both water and acetonitrile we see that there is a large enhancement of the dipole moment in going from the gas phase to the liquid phase.The enhancement for water is ϳ40% and for acetonitrile ϳ25%.From Table II we see that the dipole moment is completely determined by the z component both in the gas phase and in the liquid phase.Since the liquid phase dipole moment is obtained from an averaging over 101 configurations the fact that the other component of the dipole moment in the liquid phase is zero indicates that the averaging is close to isotropic.The standard deviations of the dipole moment is also presented in Table II and amounts to ϳ5% for both water and acetonitrile.In Fig. 1 we display the dipole moment of the individual configurations for water.As can be seen from the figure there is strong dependence on the configurations and the dipole moment oscillate between 2.2 and 2.8 debye.In this work there is no difference between the solute and the effective dipole moment because the MD simulations are done without the static electric field present.Therefore, the orientational effect due to the electric field is not accounted for in the MD simulations.In classical MD simulations the orientational effects on the dielectric constant can be obtained from the fluctuation in the dipole moments of the molecules. 79However, this approach will not work in the present QM/MM simulations since there is only one space-fixed molecule in the QM part.
The frequency-dependent polarizability components, mean value and anisotropy of water and acetonitrile both in the gas and liquid phase are presented in Table III.For the liquid phase we present both the solute and the effective properties.We also present both the average of the anisotropy ͗⌬␣͘ and the anisotropy of the average polarizability ⌬͗␣͘.First we note that the solvent effects are not very large both for water and acetonitrile.In both cases the solute properties are larger than the gas phase values.The mean value of the effective properties are very close to the gas phase values.However, for water the components of the polarizability tensor is different for the effective and the gas phase properties.The fact that the effective mean polarizability is close to the gas phase values shows that the discrete electric field in Eq. ͑11͒ is, to first order, close to zero.If we compare ͗⌬␣͘ with ⌬͗␣͘ we see that for water they are very different.The reason for this is that for water the anistropy is very small and therefore the off-diagonal tensor components becomes important.The off-diagonal elements are on average equal to zero due to the isotropic sampling and, therefore, the anisotropy of the averaged polarizability is small and close to the gas phase value.We also note that the fluctuations are slightly larger for the effective properties than for the solute properties.In the Figs.2͑a͒ and 2͑b͒ we display, respectively, the solute and effective mean polarizability of water for the individual configurations.The solute polarizability oscillates between 10 and 10.6 a.u.whereas the effective polarizability oscillations between 9.6 and 10.4 a.u.In the case of the solute polarizability all results are larger than the gas phase value whereas for the effective polarizability some of the  configurations give a polarizability smaller than in the gas phase.
In Table IV we present the frequency-dependent ͑SHG͒ first hyperpolarizability at the frequency, ϭ0.0428 a.u., for water and acetonitrile in the gas and liquid phase.For the liquid phase we present both the solute and the effective properties.First we note that for the first hyperpolarizability the solvent effects are very large.For water this leads to a change in sign for the mean first hyperpolarizability.For acetonitrile the solute and effective mean hyperpolarizability are, respectively, a factor of 5 and 6 larger than the gas phase value.In both cases we see that the fluctuations due to the different solvent configurations are very large.The first hyperpolarizability is therefore extremely sensitive to the local structure of the solvent.In Figs.3͑a͒ and 3͑b͒ we display, respectively, the solute and effective mean first hyperpolarizability of water for the individual configurations.For both the solute and the effective mean first hyperpolarizability the fluctuations are large: they oscillates between Ϫ5 and 23 a.u.again illustrating the strong sensitivity to the solvent configurations.
In Table V we present the frequency-dependent ͑EFISH͒ second hyperpolarizability for water and acetonitrile in the gas and liquid phase.For the liquid phase we present both the solute and the effective properties.For water we see that the solute second hyperpolarizability is slightly larger than the vacuum results, whereas the effective second hyperpolarizability is smaller.For acetonitrile both the solute and the effective second hyperpolarizability is larger than the vacuum results.For both water and acetonitrile the effective second hyperpolarizability is smaller than the solute in agreement with the trend found for the linear polarizability.In Figs.4͑a͒ and 4͑b͒ we display, respectively, the solute and effective mean second hyperpolarizability of water for the individual configurations.The solute mean second hyperpolarizability oscillates between 1800 and 2600 a.u.whereas the effective mean second hyperpolarizability oscillates between 1300 and 1900 a.u.The solute properties are in general above the gas phase value whereas the effective second hyperpolarizability is always lower than the gas phase value.
For all the properties we find that the pure solvent effect, i.e., the difference between the gas phase and the solute properties, increases the properties.On the other hand the induced electric field, i.e., the difference between the solute and effective properties, decreases the properties.Therefore, the electric field induced in the solvent due to the macroscopic electric field produces a screening of the electric field whereas the field from the charge distribution of the solvent molecules produces an enhancement of the electric field.The fact that the first hyperpolarizability is nearly unaffected by the induced electric field is likely to arise from the sensitivity to the short range screening. 20

C. Macroscopic response properties 1. Refractive index
We have calculated the refractive index of liquid water and acetonitrile using Eq.͑27͒ and the effective polarizability presented above.For water we used a number density of N d ϭ0.3338ϫ10 23 cm Ϫ1 and N d ϭ0.1153ϫ10 23 cm Ϫ1 was used for acetonitrile.The results for the refractive index of water are n()ϭ1.334and n(2)ϭ1.342,where ϭ0.0428 a.u., obtained from an effective polarizability of 9.95 and 10.17 a.u., respectively.The results are in good agreement with the experimentally determined refractive index at ϭ0.0428 of nϭ1.326 ͑Ref.80͒ and at ϭ0.077 of nϭ1.333. 81or acetonitrile the calculated refractive indices are n ϭ1.362, 1.377, and 1.434 at a frequency of ϭ0.000, 0.0885, and 0.1770 a.u., calculated from effective polariz-abilities of 31.02,32.16 and 36.4 a.u., respectively.The experimental results for the refractive index at ϭ0.0856 is n()ϭ1.347and n(2)ϭ1.384. 82Again, the calculated results are in agreement with the experiments although in this case the calculated results are somewhat larger than the experimental results.
It is well known that the Lorentz-Lorenz equation often ͑but not always͒ gives a good relation between the gas phase polarizability and the refractive index.As described in the theory section it is not the gas phase polarizability but rather the effective polarizabilities that should be used in the Lorentz-Lorenz equations.We have shown that by including both the solvent effects and the effect of the local field induced in the solvent due to the electric field in the calculation of the liquid phase polarizability we obtain a value of the effective polarizability close to the gas phase value.It is therefore not surprising that the refractive index are in agreement with the experimental results.

Local field factors
In order to calculate the nonlinear susceptibilities we need to consider the local field factors described in the theory section.From the refractive index calculated above we can obtain the optical local field factors given by Eq. ͑20͒.For water at ϭ0.0428 we obtain a local field factors of L ϭ1.26 and L 2 ϭ1.27.For acetonitrile the local field factor at ϭ0.000 a.u. is L ϭ1.28, and at ϭ0.0856 the local field factors are L ϭ1.30 and L 2 ϭ1.35.
Furthermore, to calculate the EFISH susceptibility in Eq. ͑40͒ we also need to consider the static local field factor given in Eq. ͑21͒.The static local field factor depends on the dielectric constant of the liquid.However, since we have not included the orientation effect due to the static electric field in the calculation it is not likely that we can calculate this quantity correctly.In fact the calculation of the dielectric constant, in particular for water is a highly complicated problem. 57Although the Onsager local field factor is introduced to account for some of the orientational effect it is not included in a consistent manner.
We can, however, calculate the Onsager local field factors from the experimental data.For water using the experimental dielectric constant of ⑀ 1 (0)ϭ78 and the refractive index nϭ1.326for water we obtain a static local field factor of L 0 ϭ1.86.We can use this local field factor to estimate the dielectric constant using the susceptibility obtained from Eq. ͑31͒.This give a dielectric constant of ⑀ (1) (0)ϳ43.This is significantly lower than the experimental value.For acetonitrile the experimental dielectric constant is 37.5 and the refractive index is 1.339.Using this we obtain a local field factor of L 0 ϭ1.85.Again using this local field factor we estimate the dielectric constant of acetonitrile to be ⑀ (1) (0) ϳ57.For acetonitrile the estimate is significantly larger than the experimental value.Since the refractive index is predicted correctly and the largest contribution to the dielectric constant is from the dipole moment it is the orientational term which is not described accurately.

THG susceptibility
Since we calculate the effective second hyperpolarizability by finite differentiations of the frequency-dependent polarizability we cannot obtain the THG second hyperpolarizability directly.The THG second hyperpolarizability can however be obtained by calculating the EFISH second hyperpolarizability at different frequency and then use dispersion formulas 2 to estimate the THG second hyperpolarizability.Here we will estimate the THG susceptibility directly from the EFISH second hyperpolarizability as The estimated THG susceptibility for water is then (3) ϭ1.07ϫ10 Ϫ14 esu at ϭ0.0428 a.u.The experimental result for water relative to the reference value of fused silica is 0.90ϫ SiO 2 (3) measured at the same frequency. 80In the original experimental work 80 a reference value for fused silica of SiO 2 (3) ϭ3.11ϫ10 Ϫ14 esu ͑Ref.83͒ was used, however, recently a value of 1.43ϫ10 Ϫ14 esu ͑Ref.84 and 85͒ has been measured and is believed to be more accurate.Adopting the latter reference value the THG susceptibility for water is (3) ϭ1.29ϫ10 Ϫ14 esu.The estimated THG susceptibility is somewhat lower but in good agreement with the experimental result.Part of the difference can be attributed to the lower dispersion arising from the EFISH second hyperpolarizability compared with the THG second hyperpolarizability.However, it is likely that the effective second hyperpolarizability is too small.The THG susceptibility of liquid water has also been calculated using different continuum and discrete local field models giving results in the range (3) ϳ1 -2ϫ10 Ϫ14 esu depending on the local field model used and in good agreement with the results presented here. 45he experimental result for acetonitrile measured at ϭ0.0239 a.u. is (3) ϭ2.54ϫ10 Ϫ14 esu using the old reference value of SiO 2 (3) ϭ2.79ϫ10 Ϫ14 esu. 83The new reference value at ϭ0.0239 a.u. is SiO 2 (3) ϭ1.15ϫ10 Ϫ14 esu. 84,85Correcting for the differences between the two reference values gives a THG susceptibilty of (3) ϭ1.05ϫ10 Ϫ14 esu.Since we have not calculated the effective second hyperpolarizability for acetonitrile at this frequency we will use the static result to estimate the THG susceptibility instead.The static effective second hyperpolarizability for acetonitrile is ␥ ϭ4891.1 a.u.Using this value the static THG susceptibility is (3) ϭ1.27ϫ10 Ϫ14 esu.The result is larger than but in agreement with the experimental result.

EFISH susceptibility
Finally, we will use the effective mean first and second hyperpolarizability to calculate the EFISH susceptibility given in Eq. ͑40͒.This gives for water an EFISH susceptibility of (3) ϭ4.1ϫ10 Ϫ14 esu at ϭ0.0428 a.u.For liquid water there has been one EFISH experiment at ϭ0.0428 a.u. 86In that work a value of (3) ϭ17.6 ϫ10 Ϫ14 esu was reported.The measurement was done in the X convention, 51 i.e., no numerical factors in the expansion of the polarization, and relative to a quartz crystal reference where a value of d 11 ϭ0.8ϫ10Ϫ9 esuϭ0.335pm/V was adopted.However, the currently accepted value for quartz is d 11 ϭ0.30pm/V. 87,88Correcting for the difference in the reference value and conventions gives (3) ϭ10.5ϫ10Ϫ14 esu (17.6ϫ2/3ϫ0.3/0.335),i.e., significantly larger than the calculated value.Since the larger contribution comes from the term depending on the dipole moment and first hyperpolarizability it is likely that the-not well describedorientational contribution is responsible for the difference.The agreement between theory and experiment for the THG susceptibility also support this conclusion.
For acetonitrile we calculate the EFISH susceptibility at ϭ0.0885 a.u. to be (3) ϭ35.7ϫ10Ϫ14 esu.The experimental result 89 at the same frequency is (3) ϭ13.6 ϫ10 Ϫ14 esu ͑corrected with 2/3 to convert it to the convention used here͒.The calculated value is much larger than the experimental results in agreement with the trend found for the dielectric constant.There has recently been a study of both the microscopic and macroscopic properties of acetonitrile using ab initio method combined with an Onsager model for the solvation where in general a good agreement was found between experiments and the calculated results. 46

V. CONCLUSIONS
We have in this work presented an extension of the QM/MM formalism to include the so-called local field factors, i.e., the difference between the macroscopic electric field and the actual electric field felt by the solute molecule.This enables the calculation of effective microscopic properties which can be related to the macroscopic susceptibilities directly comparable with experimental results.By separating the discrete local field into two distinct contribution we can define two different microscopic properties, the so-called solute and effective properties.In the solute properties the pure solvent effect, i.e., effects even when the macroscopic electric field is zero, are accounted for and in the effective both the pure solvent effect and the effect from the induced dipoles in the solvent is accounted for.We have presented results for linear and nonlinear polarizabilities of water and acetonitrile both in the gas phase and in the liquid phase.For all the properties we see that the pure solvent effect, i.e., the difference between the gas phase and the solute properties, gives an increase in the properties.For the induced electric field, i.e., the difference between the solute and effective properties, a decrease in the properties was found.Therefore, the electric field induced in the solvent due to the macroscopic electric field produced a screening of the electric field whereas the field from the charge distribution of the solvent molecules produces an enhancement of the electric field.Furthermore, we have presented results for the refractive index, third-harmonic generation ͑THG͒, and electric field induced second harmonic generation ͑EFISH͒ for pure water and acetonitrile.We find in general good agreement between the calculated and experimental values for the refractive index and the THG susceptibility.For the EFISH susceptibility the differences between experiment and theory is larger due to neglect of the orientational effect arising from the static electric field.

FIG. 1 .
FIG. 1.The dipole moment of water in Debye.Top line is the averaged results and the bottom line is the gas phase value.
a Results are taken from ͒ b Results are taken from ͒ c TABLE II.The dipole moments of water and acetonitrile in gas and liquid phase in atomic units.

TABLE III .
The frequency-dependent polarizability of water and acetonitrile in the gas and liquid phase in atomic units.For water the frequency is ϭ0.0428 and for acetonitrile ϭ0.0885 a.u.

TABLE IV .
The first hyperpolarizability ␤(Ϫ2;) for water and acetonitrile in the gas and liquid phase.For water the frequency is ϭ0.0428 and for acetonitrile ϭ0.0885 a.u.All results are in atomic units.

TABLE V .
The second hyperpolarizability, ␥(Ϫ2;,0), for water and acetonitrile in the gas and liquid phase.For water the frequency is ϭ0.0428 and for acetonitrile ϭ0.0885 a.u.All results are in atomic units.