J. Phys. Chem. 1996, 100, 17791-17796
17791
Electronic Polarization of a Water Molecule in Water. A Combined Quantum Chemical and Statistical Mechanical Treatment Nigel W. Moriarty and Gunnar Karlstro1 m* Department of Theoretical Chemistry, Chemical Center, P.O.B. 124, S-221 00 Lund, Sweden ReceiVed: May 21, 1996; In Final Form: August 2, 1996X
The polarization of an ab initio quantum chemically described water molecule due to the surrounding statistical mechanical Monte Carlo water bath is studied. The quantum chemically described water molecule is surrounded by 89 statistically described water molecules situated in a spherical cavity in a dielectric continuum. The total interaction between the quantum chemically described water molecule and the total surroundings is on average 0.5 kJ/mol larger than for the interaction for a classically described water molecule. This increase is despite the fact that the potential function used between the quantum chemically described water molecule and a classical water molecule is almost identical to the potential between the classical molecules. The difference is due mainly to the polarization of the quantum chemically and the classically described water molecules.
1. Introduction In recent years both quantum chemical and statistical mechanical techniques have been developing with great speed. Today it is possible to describe relatively complex reaction processes in the gas phase and the behavior of large proteins in water by quantum chemical methods. These quantum chemical methods provide for the interaction of the molecules by the variation of the wave function of the entire system. This is in contrast to a basic assumption in the statistical mechanical simulations that the molecules remain largely unchanged regardless of their interaction with other molecules in the system. As a consequence, it is not possible to effectively study processes where the interacting molecules change their character due to this interaction using statistical methods. A first step to circumvent this problem has been taken in the so-called polarizable models, like the NEMO approach,1-8 where each molecule or molecular segment is polarized due to the electric field experienced. One must keep in mind, however, that the models used to describe the polarization are normally very simple; spherically isotropic polarizability tensors omitting all hyperpolarizabilities. An alternative procedure would be to treat the central part of the system (the molecule or molecules of interest) quantum mechanically and treat the surrounding by some statistical mechanical means. Methods of this type, where the quantum chemically described system is treated with a semiempirical quantum chemical method, have been developed9-14 (including reviews15-20). At least two closely related major problems may be identified in this context. Normally the semiempirical quantum chemical models are parameterized to describe molecules at equilibrium gas phase conditions where experimental data are frequent, but in this type of application one is often interested in nonequilibrium properties and conditions. The other problem is that for each new part of the configurational space that is investigated a new quantum chemical calculation, where the coupling between the surrounding and the quantum chemical system is included, must be performed. This requires large computer resources and has prevented ab initio quantum chemical methods from being used as the quantum chemical * Author for correspondence. X Abstract published in AdVance ACS Abstracts, October 15, 1996.
S0022-3654(96)01476-1 CCC: $12.00
method. Nevertheless, in order to make a theoretical description of chemical reactions occurring in condensed phases possible, it is necessary to embark on such a research direction. In this work we will describe such an effort and apply it to study a quantum chemically described water molecule surrounded by 89 classically described water molecules. The entire system is surrounded by a spherical cavity in a dielectric continuum. 2. The Model Most often when one wishes to describe a chemical process occurring in a condensed phase, it is easy to identify a central part of the system where the action takes place: where the bonds are broken and formed or where the electron is released and absorbed. It is obvious from a theoretical point of view that these processes must be described by quantum chemical means, preferably with an ab initio method since semiempirical methods are not parametrized for the purpose of describing chemical reactions in a solvent. The central quantum chemical part, which we will call the kernel of the system, couples to the surrounding molecules both by electrostatic couplings and via the Pauli principle. Consequently, the kernel of the system must experience the electrostatic potential and the field from the surrounding molecules as well as the exchange repulsion due to the electrons in the surrounding molecules. In most approaches of this type formulated thus far only the electrostatic part of the coupling has been considered. At the same time the surrounding molecules must experience the electric field from the kernel of the system. Here we can identify a problem. Our intention is that the surrounding molecules are supposed to be classical, and thus no wave functions are needed to describe them; but in order to describe their exchange interaction with the kernel, a wave function for the classical part is needed. This addition modifies the quantum Hamiltonian with the overlap contributions of the surrounding water molecules. The total molecular system (the quantum mechanically and classically described molecules) is embedded in a cavity in a dielectric continuum. The total system is visualized in Figure 1. We have chosen to embed the molecular system in a continuum rather than use periodic boundary conditions, as is normally done in statistical mechanical simulations. This strategy has been chosen since there is a natural origin in the © 1996 American Chemical Society
17792 J. Phys. Chem., Vol. 100, No. 45, 1996
Figure 1. Diagram of the total system including kernel molecule, classical molecules, and dielectric continuum.
system where large charge displacements may occur, and the boundary conditions have been been chosen in order not to introduce any artificial couplings between different system regions. Naturally, implicit in this strategy is the wish to keep the kernel in the center of the cavity to minimize the effects of the change from a molecular fluid to a continuum. However, since the number of solvent molecules is limited, artifical packing effects may be introduced on these solvent molecules by the boundary, and these packing effects may be enhanced by the presence of a fixed molecular system in the center of the cavity. This problem can be solved by fixing the kernel to a spring anchored at the center of the cavity. Analysis of the physics of the system presented in Figure 1 shows there are three different types of theoretical descriptions that must be linked together. These are a quantum mechanical region, a classical discrete molecular region, and a classical continuum region. An outline of the formalism used to describe each of these regions as well as how these regions are linked to each other is the basis of the next section. 3. The Formalism This fairly lengthy section is divided into four subsections. The first two of these subsections will deal with the description of the kernel and of the discrete molecular surrounding, respectively. The third subsection will describe how the quantum mechanical system is linked to the discrete molecular system. Finally in the fourth subsection we will present how these two parts are linked to the continuum. As will be seen below, the main problem will be to make reasonable compromises between the accuracy of the calculations and the fact that only a limited amount of computer resources are available. 3.1. The Quantum Chemical Kernel. The method that is presented below is based on the simplest possible ab initio level of approximation: The Hartree-Fock level. The actual description given below refers to the particular problem studied here, the polarization of a water molecule due to other water molecules, but it can be generalized to any molecular system. The first step in any quantum chemical calculation is to choose a basis set and a geometry. In this work we intended to link the kernel of the system to a statistical mechanical region where water molecules are described with what we have call a NEMO potential.1-8 These potentials are obtained with an ANO basis set21 constructed from 10s, 6p, and 3d functions which
Moriarty and Karlstro¨m are contracted to 4s, 3p, and 1d on oxygen and a 6s4p basis set contracted to 3s and 2p on hydrogen. The O-H bond length used in the NEMO approach has the value 0.9572 Å, and the ∠HOH angle a value of 104.52°. For the sake of consistency we would like our quantum chemical system to resemble the classical water molecules as much as possible. Thus we choose the same geometry and basis set for this part as was used in the construction of the NEMO potential. Since the molecule will interact with the surrounding molecules, it cannot be expected to possess any particular symmetry. The total basis set consists of 37 basis functions, and the corresponding Fock matrix will contain 703 elements. This means that if a supermatrix formalism is used to construct the Fock matrix, 247 456 different supermatrix elements must be processed in order to construct a Fock matrix. (A supermatrix element (PQRS) is defined by the ordinary two-electron integrals (pqrs) according to (PQRS) ) (pqrs) - 0.25*((prqs + (psqr)).) On the average five iterations are needed to solve the Hartree-Fock (HF) equations for one statistical mechanical configuration, meaning that approximately 106 operations must be performed on the computer in the Fock matrix construction for one particular configuration. Thus a large part of the computational time would be used to construct the Fock matrix, and if the kernel was twice as big, it would be impossible to perform the calculation. There is consequently an urgent need to reduce the size of the basis set and the Fock matrix. In some other approaches it has been tried to perform a new quantum chemical calculation only after a specific number of new conformations for the classical system have been generated. This approach has the disadvantage of not responding to the changes in the wave function due to the changing environment. A different solution involving the HF calculation at every configuration but reducing the need for computer resources is required. To describe the water molecule reasonable well, it is desirable to retain a basis set of the above quality. In this particular case we are interested in the polarization of the water molecule, and one way may be to design a basis set suitable particularly for the purpose of describing the polarization of the water molecule. In order to do so, a set of calculations were performed where the water molecule was subjected to different perturbations such as a field in different directions or a potential acting on the electrons of the type
Uel ) A(1 - e-(r/r0) )12 2
(1)
In eq 1, A and r0 are constants and r is the distance from the electron to the center of the cavity. Such potentials have previously been applied in ab initio quantum chemical dielectric cavity calculations to prevent the electrons from penetrating out into the dielectric medium.22 Seven different perturbations with a field in the positive and negative axis directions (0.01 au) and one where the water molecule was confined to a cavity with r0 ) 7.7 au and an A value of 100 were used. The density matrices from these calculations were added, in a fashion analogous to the ANO basis sets construction,23 and the resulting density matrix was diagonalized. The 13 orbitals with the highest occupation numbers were selected (the highest not included occupation number was on the order of 10-7). The one-electron integrals and the supermatrix elements were then transformed to this molecular orbital basis, leading to a reduction of the number of Fock matrix elements to 91 (13% of the original size) and a reduction of the number of supermatrix elements to 4596 (1.9% of the original value). By this technique the loss of accuracy can be kept very small and substantially
Electronic Polarization of a Water Molecule in Water reduce the time used for handling the intramolecular part of the quantum chemical subsystem. Thus in summary, in the absence of the surrounding, an ordinary Hartree-Fock calculation with the transformed approximate molecular orbital basis is performed. 3.2. The Statistical Mechanical Region. The purpose of this region is to link the kernel to the continuum and give a good description of the short and intermediate range part of the solvation of the kernel. This is achieved by modeling the medium by an ordinary Monte Carlo procedure, where the new configurations are chosen according to the procedure of Metropolis et al.24 In all, 89 classical water molecules are located in the dielectric cavity, together with the quantum chemically described water molecule. In a set of reference calculations with only classical molecules in the cavity 90 water molecules were used. The NEMO potential in the original form was used2 to describe the interaction between the classical molecules. This means that each water molecule has three centers that carry a polarizability, and the energy calculation needs to be solved self-consistently (see ref 2) even for classical systems. This also means there is very little to gain by storing the additive part of the molecular interaction energies, and an alternative way of moving the molecules has been tested. In this procedure all molecules are moved in every new configuration (obviously the displacement parameters must then be given smaller values than when moving only one molecule). In this way, the configuration space is sampled more efficiently. It is, however, not the purpose of this work to investigate this problem. Apart from the multimolecule displacement and the boundary condition, the treatment of this region is similar to ordinary Monte Carlo simulations. A similar approach but with only classical molecules has been tested by Wallqvist.25 3.3. Coupling between the Classical and the Quantum Mechanical Region. In the Introduction it was mentioned that there are two types of interactions between the “classical molecules” and the “quantum mechanical molecule”: electrostatic and exchange repulsion interactions. The electrostatic properties of each of the classical molecules is described with a set of point charges and a set of polarizabilities. The ideal approach would obviously be to calculate the interaction between these point charges and the dipoles generated by the polarizabilities and the electronic and nuclear density of the quantum molecule. The density matrix for the kernel molecule has the same size as the Fock matrix, 703 in the contracted basis set and 3655 in the uncontracted basis set. Each classical molecule consists of 4 point charges and 3 (with 3 directions) polarizabilities, and there are 89 classical molecules. In all, a total of over 4 000 000 interactions need to be calculated for each configuration, constituting a severe bottleneck in the calculations if not approximated. When the NEMO potential between the classical molecules was constructed,1 the charge density of the classical molecules was expanded into a multicenter multipole expansion.26,27 Formally there exists one center for each atom and for each formal bond (there is a formal bond between each pair of atoms). The density corresponding to a pair of basis functions is assigned to a center connected with the basis functions in the following way. If both basis functions are centered on the same atom, then this atom is the expansion center. If the basis functions belongs to different atoms, then the expansion center is chosen as the formal bond between the atoms. For more details see refs 26 and 27. In each of the expansion points a point charge, three unit dipoles, and the five quadrupole tensor components are located. The interaction between these quantities and the classical molecules are used in the calculations as an ap-
J. Phys. Chem., Vol. 100, No. 45, 1996 17793 proximation of the complete number of interactions mentioned earlier. This reduces the number of elements representing the quantum mechanical molecule from 3655 to 60. This greatly reduces the number of interactions (from 4 000 000 to 65 000) and makes the calculation far more tractable than using the densities of the quantum molecule. The density corresponding to each of the primitive density matrix elements is then expanded in the multipole moments of the multicenter multipole expansion. Finally the representation of the densities generated from the molecular orbitals, which are used in the actual calculations, in terms of the multicenter multipole expansion is calculated. It should be noted that the multicenter multipole expansion is only calculated once and not for each new Monte Carlo configuration. The polarizabilities of both the classical and quantum molecules will mutually polarize each other, inducing dipoles. These must be determined in a self-consistent fashion, and this is done concurrently with the self-consistent field of the density in the quantum molecule. In the previous section, we accounted for the electrostatic coupling between the classical and the quantum molecules. However there is also a force originating from the Pauli principle, normally termed exchange repulsion, due to the overlap between the wave function of different molecules. To calculate this quantity, each classical molecule is embued with a wave function used to determine the exchange repulsion when the distance between a classical molecule and the quantum kernel is less than 7.7 au. The wave function of the classical molecules is represented with a set of orbitals; in this case the lowest five orbitals of the quantum system. When the exchange is turned on, the overlap intergrals between the “classical molecule” and the kernel are calculated and an extra term is added to the Hamiltonian describing the kernel according to
hij ) hij(1 + rij)
(2)
where rij is a correction term that depends on the overlap between the wave functions of the quantum system and the wave function corresponding to the classical molecules. To describe the calculation of the rij, we begin by noting that the ideas behind the definition of this quantity are similar to the ideas behind the description of the exchange repulsion in the NEMO model. We start by introducing the quantities Sik and sik(m,n). Sik is the overlap between orbital i on the quantum mechanical system and occupied orbital k in the classical system. sik(m,n) is the contribution to this overlap from basis functions centered on atom of type m in the quantum system and atom of type n in the classical system. Using these quantities, we may define the contribution to rij from occupied orbital j in the classical part of the system according to
∆rij ) ζ(k)∑{Siksjk(m,n) + Sjksik(m,n)}d(m,n)
(3)
n,m
where ζ(k) is the orbital energy for orbital k and d(m,n) are parameters (given in Table 1) defining the importance of a certain type of overlap for the exchange repulsion. The total rij is then obtained by summing up the ∆rij for all occupied orbitals on molecules in the statistical mechanical part of the system. In practice, only molecules located closer to the kernel than a certain cutoff distance (7.7 au) are considered. This approximation neglects interactions that give rise to approximately 0.01 kcal/mol. In addition to the exchange repulsion terms involving Sij2, a small additional term is added to the total energy of the system according to
17794 J. Phys. Chem., Vol. 100, No. 45, 1996
Moriarty and Karlstro¨m
TABLE 1: Parameters of the QMSTAT Program name
value
exchange repulsion-Fock matrix contribution exchange repulsion-d(H,H)a exchange repulsion-d(O,H)a exchange repulsion-d(H,O)a exchange repulsion-d(O,O)a exchange repulsion-Rb hydrogen bond exponent-ac hydrogen bond factor-bc O-O separation-cd O-H separation-cd
-0.014 027 4.771 842 1.3 1.0 1.0 1079.966 454 -0.098 069 1.635 431 3.49 1.84
a
From eq 3. b From eq 4. c From eq 5. d From eq 6. extra ∆Eexrep ) R∑Sij6
(4)
i,j
This term does not effect the Fock matrix, but is introduced in order to make the potential curves steep enough for short distances. Apart from the terms above, all of which correspond to interactions described at the Hartree-Fock level of approximation, there is a dispersive coupling between the classical molecules and the quantum mechanical kernel. This term is added in a phenomenological way and is similar to the dispersive coupling between the classical molecules.1 Apart from the terms described above an extra correction term was needed in order to make the hydrogen bonds donated and accepted by the quantum mechanical system equally strong. This term has the form
ae(-bROH)
(5)
added to the potential energy when the quantum mechanical system is acting as hydrogen bond donor and subtracted when it is acting as acceptor. The term is small: 0.09 kcal/mol slightly inside the potential energy minimum. One final inclusion is a non-Lennard-Jones repulsive term of the form
(cij/Rij)21
(6)
designed purely to improve the OQM-OSM separation during the simulation. Typically, the increase in repulsion even at relatively small O-O bond separations, say 2.5 Å, is less than 1.1 kcal/mol, but at the water dimer optimal bond length it is less than 0.06 kcal/mol. The O-H term is included in a similar fashion as in the NEMO model, while the H-H term is set to zero. The O-H term has very little effect on the simulation of bulk water. With these parameters the potential functions between a classical and a quantum mechanical water molecule are almost identical with respect to variations of the O-O distances for optimal values of angular parameters. 3.4. The Coupling with the Dielectric Medium. Finally we will now describe the coupling between the classical and quantum systems with the dielectric medium. Friedman28 suggested the possibility of including the effect of a dielectric surrounding in statistical mechanical simulations by using the image charge approximation. This technique has been used in this work. In the image charge approximation the response of the dielectric medium to a charge in the cavity is represented by a charge in dielectric medium, a dipole in the cavity leads to a charge and a dipole in the medium, and so on. To facilitate the calculations only terms originating from charges and dipole moments are considered in this work. This is no approximation for the classical molecules since they are represented by charges
Figure 2. Energy potential in hartrees and angstroms between one quantum water molecule and one classical water molecule using the QMSTAT model. The dimer is in the minimum configuration, and the legend refers to the quantum molecule.
and polarizabilities (dipoles), but this is an approximation for the kernel. However, since the kernel is far from the dielectric boundary, one may expect this effect to be very small. In short one can say that each molecule in the cavity gives rise to an image molecule in the dielectric and that the size of the system consequently is doubled by the use of this technique. The model presented has been implemented in a program called QMSTAT. 4. Discussion and Results The main reason for parameterizing the interaction between a classical and a quantum mechanical water molecule in order to resemble the interaction between two classical or two quantum mechanically described molecules was to facilitate the interpretation of the obtained results. In principle, three types of effects on the structural data obtained in this study could be expected. Firstly, one may expect differences caused by differences in the pair potentials. (It was the purpose of the work described in the previous section to minimize these effects.) Secondly, one may expect an influence from the dielectric discontinuity. This effect can be studied by comparing the results from a purely classical system in the cavity with the results obtained from a periodic boundary condition study.2 One may expect this effect to be small if the results from a purely classical simulation are compared with the results from our combined system. Thirdly, one may expect many-body effects present in the quantum chemical system and not present in the classical system to be of some importance. Here it should be noted that the most important part of the many-body effect (the induction) is described in both the classical and the quantum mechanical system. The parametrization of the model was performed on the water dimer. One water molecule was treated quantum mechanically, while the other was the classical model from the NEMO model. The quantum molecule was tested in both the hydrogen donor and acceptor role. The agreement between the two energy potential curves displayed in Figure 2 is sufficient for the model. The minimum for each of the curves is greater than or equal to
Electronic Polarization of a Water Molecule in Water
Figure 3. Pair correlation function in arbitrary units versus angstroms for the quantum oxygen and the classical oxygen atoms. The NEMO results3 are included for comparison.
Figure 4. Pair correlation function in arbitrary units versus angstroms for the oxygen and hydrogen atoms. The NEMO results3 are included for comparison.
2.9 Å. The steepness of the curve at shorter distances than the minimum is critical to the correct simulation of the bulk water, and several simulations were performed to determine the suitability of this area of the surface. A Monte Carlo simulation of 2000 macrosteps (after equilibration) using the parameters determined via the dimer calculations was performed using one quantum water molecule and 89 classical water molecules. The total interaction of the combined quantum and statistical system was 0.5 kJ/mol larger than the purely classical system. There are large variations between configurations, giving a standard deviation of nearly 6 kJ/mol. These variations are between consecutive steps and separate regions of the simulation. The pair correlation functions, g(r), of the QMSTAT simulation displayed in Figures 3, 4, and 5 give the average distance between the quantum oxygen and the classical oxygen, gOQMOSM(r), as 2.7 Å. This is in contrast to the dimer results indicating a slightly too large bond length for the O-O distance (the subscripts of the atoms are omitted from subsequent discussion by using the convention of placing the quantum atom first in the subscript). This shortening is due to the many-body effects.
J. Phys. Chem., Vol. 100, No. 45, 1996 17795
Figure 5. Pair correlation function in arbitrary units versus angstroms for the quantum hydrogen and the classical hydrogen atoms. The NEMO results3 are included for comparison.
There is some additional structure on the QMSTAT curves which may be due to the reduced sampling available to the system compared to a purely classical system. The g(r) does not approach a constant value at large separation due to the geometry of the simulated system. This makes the normalization of the curves impossible. The agreement between the results of Wallqvist et al.3 is good for distances less than 4.5 Å. For large distances one may expect an influence from the dielectric medium. To investigate the effect of the dielectric medium, a Monte Carlo simulation was performed with a cavity radius 5% larger. The original radius choice was based on Monte Carlo simulations in which the radius was allowed to fluctuate as part of the random changes. It was noted that the radius changed very slowly, and a suitable value was chosen for all subsequent simulations. The increased radius had no effect on the g(r) results at distances smaller than 4.5 Å. For larger distances, the results were not improved. The pair correlation functions for the O-H distances, gOH(r) and gHO(r), are in reasonable agreement (see Figure 4). This agreement is lost further than 3.5 Å from the quantum water. Again it is likely that the extra structure found in the gOH(r) and gHO(r) curves is due to the dielectric surrounding. The leftmost peak is very well related and the second peak has reasonable correlation with the experimental results despite the statistically inferior results further out. The statistics may also give rise to the large differences in height between the two QMSTAT curves displayed. The H-H pair correlation function (see Figure 5) also agrees well with the calculated results of Wallqvist et al.3 The agreement extends to 4.5 Å and includes the first and second peaks and the intervening trough. The relative heights are also reasonable. The dipole moment of the quantum water molecule is shown in Figure 6. The peak of a Gaussian fit is at 2.88 D, and the exponent is 1/(0.36)2. This result is 0.02 D higher than NEMO calculations. This is also the difference between the QMSTAT and NEMO calculations for the monomer: 2.06 and 2.04 D, respectively. The consistently higher dipole is expected and is due to the more diffuse basis set. The electronic polarization of a water in water is 0.82 D using the QMSTAT method. The angle between the dipole of the molecule along the bisector of the ∠HOH angle (fixed dipole) and the total dipole
17796 J. Phys. Chem., Vol. 100, No. 45, 1996
Moriarty and Karlstro¨m The success of the combined quantum chemical and statistical mechanical treatment to reproduce pair correlations, g(r), dipole moment, and electronic polarization results of the NEMO calculations gives credence to the method. The close agreement of NEMO with experiment1 adds confidence and suggests that useful results can be obtained from the method suggested in this paper. The main problem in using the suggested method probably originates in the parameterization of the interaction between the classical solvent molecule and the quantum system. References and Notes
Figure 6. Dipole moment of the quantum water molecule in debyes including a Gaussian curve fit.
Figure 7. Angle between the total dipole moment of the quantum water molecule and the fixed dipole, µz, in degrees including a Gaussian curve fit.
is shown in Figure 7. The peak of the Gaussian is at 6.5°, and the half-width is 2.4°. This shows that the total dipole is almost parallel to the bisector. 5. Conclusions In this work we have shown that using some straightforward approximations, it is possible to model a molecular system using ab initio quantum chemical methods in a statistical mechanical surrounding. The quantum chemical molecule produces a 3-fold increase in computer time over the purely classical system, a much smaller increase than typical in ab initio calculations.
(1) Wallqvist, A.; Karlstro¨m, G. Chem. Scr. 1989, 29A, 131. (2) Wallqvist, A.; Ahlstro¨m, P.; Karlstro¨m, G. J. Phys. Chem. 1990, 94, 1649. (3) Wallqvist, A.; Ahlstro¨m, P.; Karlstro¨m, G. J. Phys. Chem. 1991, 95, 4922. (4) Wallqvist, A.; A° strand, P.-O. J. Chem. Phys. 1995, 102, 6559. (5) A° strand, P.-O.; Karlstro¨m, G. J. Phys. Chem. 1991, 95, 6395. (6) A° strand, P.-O.; Wallqvist, A.; Karlstro¨m, G.; Linse, P. J. Chem. Phys. 1991, 95, 8419. (7) A° strand, P.-O.; Wallqvist, A.; Karlstro¨m, G. J. Chim. Phys. 1991, 88, 2457. (8) A° strand, P.-O.; Wallqvist, A.; Karlstro¨m, G. J. Chem. Phys. 1994, 100, 1262. (9) Fois, E. S.; Sprik, M.; Parinello, M. Chem. Phys. Lett. 1994, 223, 411. (10) Wei, D.; Salahub, D. R. Chem. Phys. Lett. 1994, 224, 291. (11) Muino, P.; Callis, P. R. J. Chem. Phys. 1994, 100, 4093. (12) Gao, J. J. Am. Chem. Soc. 1994, 116, 9324. (13) Car, R.; Parrinello, M. J. Chem. Phys. 1985, 55, 2471. (14) Gao, J.; Xia, X. Science 1992, 258, 631-635. (15) Car, R.; Parrinello, M. In Simple Molecular Systems at Very High Density; Polian, A., Ed.; Plenum: New York, 1989; p 455. (16) Payne, M. C.; Bristowe, P. D.; Joannopolous, J. D. Mol. Simul. 1989, 4, 79. (17) Payne, M. C. J. Phys. Condens. Matter 1989, 1, 2199. (18) Remler, D. A.; Madden, P. A. Mol. Phys. 1990, 70, 921. (19) Gao, J. In ReViews in Computational Chemistry; VCH: New York, 1995; pp 119-185. (20) Warshel, A. Computer Modeling of Chemical Reactions in Enzymes and in Solutions; John Wiley & Sons: New York, 1991. (21) Pierloot, K.; Dumez, B.; Widmark, P.; Roos, B. O. Theor. Chim. Acta 1995, 90, 87. (22) Karlstro¨m, G. J. Phys. Chem. 1989, 93, 4952. (23) Widmark, P.-O.; Malmqvist, P.-A° .; Roos, B. O. Theor. Chim. Acta 1990, 77, 291. (24) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, J. J. Chem. Phys. 1953, 21, 4952. (25) Wallqvist, A. Mol. Simul. 1993, 13, 101. (26) Karlstro¨m, G. On the evaluation of intermolecular potentials. In Proceedings of 5th Seminar on Computational Methods in Quantum Chemistry, van Duijnen, W. C., Nieuwpoort, P. T., Eds.; Lab. of Chem. Phys., Univ. of Gronningen, Nijenborgh 16, 9747 AG Gronningen, The Netherlands, 1981; p 353. Max-Plank-Institut fu¨r Physik und Astrophysik; Institut fu¨r Astrophysik, Munich. (27) Karlstro¨m, G.; Linse, P.; Wallquist, A.; Jo¨nsson, B. J. Am. Chem. Soc. 1983, 105, 3777. (28) Friedman, H. L. Mol. Phys. 1975, 29, 1533. (29) Coulson, C. A.; Eisenberg, D. Proc. R. Soc. London 1966, A291, 445. (30) Soper, A. K.; Phillips, M. G. Chem. Phys. 1986, 107, 47.
JP9614761
