2504
The Journal of Physical Chemistry, Voi. 82, No. 23, 1978
R. A. Nemenoff, J. Snir, and H. A. Scheraga
A Revised Empirical Potential for Conformational, Intermolecular, and Solvation Studies. 2. Parameterization and Testing for Water and Saturated Organic Molecules' Raphael A. Nemenoff,28Joseph Snir, and Harold A. Scheraga*2b Department of Chemistty, Cornell University, Ithaca, New York 14853, and Biophysics Department, WeizmannInstitute of Science, Rehovoth, Israel (Received December 21, 1977; Revised Manuscript Received August 9, 1978)
The revised model of the EPEN potential (empiricalpotential based on the interactions of electrons and nuclei), EPENI2, has been parameterized for water, saturated hydrocarbons, alcohols, ethers, and amines. For water, the potential was parameterized by fitting to theoretical values of the energies of various dimer configurations. Using a Monte Carlo technique, with this potential, reasonable results were obtained for the thermodynamic properties of liquid water; however, the potential lacks the longer range tetrahedral order necessary for a proper description of the structural properties. For saturated organic molecules, the potential was parameterized by fitting to gas-phase dipole moment and conformation data, and to crystal data. The potential was tested by performing additional calculations on gas-phase conformations and crystal structures. Good agreement with experiment was obtained for the dihedral angles of stable conformers, relative energies between conformers, barriers to internal rotation, and dipole moments. Crystal structure calculations were performed with the maximum number of degrees of freedom consistent with maintaining the observed symmetry. The revised model (EPEN/2) gives better agreement with experiment than did the original version of this potential function.
I. Introduction It has been shown3 that the EPEN potential (empirical potential based on the interactions of electrons and nuclei), which had been parameterized for saturated organic molecule^,^^^ could not reproduce the structural and thermodynamic properties of liquid water. The reasons for these difficulties, which are critical for a potential developed for studies of solvation, have been discussed, and a revised model (EPEN/2) has been p r ~ p o s e d .In ~ this paper, we present the results of the reparameterization of EPEN (with this new model) for water, saturated hydrocarbons, alcohols, ethers, and amines. In subsequent the potential is extended to unsaturated molecules. In section 11, the development of the parameters for water will be presented. In section 111, the parameterization of the remaining saturated systems is discussed. Tests of the new potential are described in section IV, and a summary and conclusions are given in section V. Most of the calculations were performed on the Golem B computer at the Weizmann Institute; the remaining ones were carried out on the IBM 370/168 computer at Cornel1 University. 11. Development of Parameters for Water The previous version of EPEN provided a.good representation of the interaction between water molecules for low-energy configurations where an 0-H-0 hydrogen bond is present, as in ice4 and in small clusters of water molecule^.^ However, for orientations of water molecules where lone-pair-lone-pair interactions dominate (see Figure 1 of ref 3), the potential predicted an incorrect attractive behavior. Since a good description of the entire potential energy surface is required in order to model the properties of a liquid, errors in these regions would lead to incorrect structural and thermodynamic results in a Monte Carlo simulation of the liquid. When such calculations were performed using the original EPEN potential, the system was found to collapse with a large negative internal energy, and constantly decreasing volume.1° As discussed in paper 1 of this ~ e r i e sto , ~describe the complete potential energy surface of two water molecules adequately, a potential function should be parameterized on data representative of the complete surface. Such data 0022-3654/78/2082-2504$0 1.OO/O
can be supplied only by quantum mechanical calculations. Matsuoka et a1.l1 have performed ab initio quantum mechanical calculations with an extended basis set, and configuration interaction, for 66 water-dimer configurations. The dimer energies were fitted to an empirical potential, which will be designated here as the MCY potential. Monte Carlo calculations were performed in both the NVT12and NPT13ensembles, where N , V, P, and T are the number of molecules, volume, pressure, and temperature, respectively, and the potential was found to give a reasonably good description of both the structural and thermodynamic properties of the liquid. The model used in the MCY potential is specific for water, and not readily transferable to other molecules. Since it is of interest to study molecular conformation of a solute molecule in liquid water, it is necessary to have a potential which treats water and solute molecules in a consistent way. We have used the same data for the 66 dimer configurations, and fit the energies to the revised EPEN model, which is readily extendable to solute molecules. The EPEN model for a water molecule has seven centers of interaction. The oxygen nucleus has a charge of +6 au, and each proton a charge of +1au. We have used the same geometry for the water molecule as Matsuoka and coworkers.ll (The 0-H bond length is 0.9572 A, and the H-0-H bond angle is 104.52O.) There are four electron pairs in the molecule; two 0-H bonding electron pairs are constrained to lie along the 0-H bonds. The two lone pairs are constrained to lie in the plane perpendicular to the plane of the molecule, containing the H-0-H bisector, and are placed symmetrically with respect to the molecular plane (see Figure 1). There is a Coulombic interaction between all point charges of the form qiq,j Ecoul = 332.0719Z - kcal/mol i j rij where q i and qj are the charges in atomic units, and rij is the distance between qi and q , in angstrom units. There is a nonbonded Buckinghamli4 interaction only between the electrons of the form electrons
0 1978 American Chemical Society
The Journal of Physical Chemistry, Vol. 82, No. 23, 1978 2505
H20 and Saturated Organic Molecules
Figure 1. EPEN/2 model for the water molecule.
where A , B, and C are nonbonded parameters. In the original EPEN all electrons in the water molecule had been assigned the same set of parameters for the nonbonded potential. To provide the potential with greater flexibility, a distinction is made here between bonding and lone pair electrons, and each type of electron is assigned a separate set of A , B, and C parameters. For interaction between the bonding-pair and lone-pair electrons, the following combination rules were used:
All = (AitAJ'12 Bij = (Bit
+ Bjj)/2
CLj = (C,,Cjl)l"
(34 (3b) (3c)
where i and j refer to bonding or lone-pair electrons. There are thus nine parameters to be determined by fitting to the dimer energies: the distance of the 0-H bonding electron pairs from the 0 nucleus (along the 0-H bond); the distance of the lone pairs from the 0 nucleus, constrained to be symmetric and in a plane perpendicular to the molecular plane; the lone pair-0-lone pair angle; A , B, and C for the bonding electrons; and A , B, and C for the lone-pair electrons. There are ten adjustable parameters in the MCY potential. The nine parameters were derived by minimizing the square of the deviation of the EPENIB intermolecular energy from the ab initio energy, Ehk,, of Matsuoka et al." for the 66 dimer configurations, plus the square of the deviation of the calculated dipole moment from the experimental one, using the algorithm of P0we1l.l~ The dipole moment was calculated using the following expression p =
4.80325X qlRl (debye)
Figure 2. Repulsive configurations of water molecules which were examined with EPEN and the MCY potential.
(4)
1
where Ri is the vector position (in angstrom units) of point charge i, and p is the dipole moment vector. A leastsquares rms deviation of 0.35 was obtained; Matsuoka et al.ll obtained a value of 0.184 for the MCY potential, compared to their a b initio values. To test the new EPEN potential function, a Monte Carlo calculation was performed on 64 water molecules in the N P T ensemble, using periodic boundary conditions. Although the calculation was not carried to convergence, it was clear that this version of EPEN did not lead to a good description of the liquid. In particular, the molar volume was much too low, the internal energy was too negative, and the number of nearest neighbors was much too high. To analyze the problem, EPEN was compared with the MCY potential in regions of the dimer potential energy surface which had not been used in parameterization. Figure 2 shows the configurations examined, and Figure 3 is a comparison of potential energy vs. O.-O distance for both EPEN and the MCY potential for one of these configurations. Although the EPEN potential was repulsive in all regions where the MCY was repulsive, the EPEN curve was less steep than that of the MCY potential. Thus, even though EPEN and MCY agreed for the 66 quantum-mechanically-calculated dimer energies
2
Flgure 3. Dimer energy as a function of 0-0 distance for configuration b of Figure 2. Curve 1 is the EPEN energy; curve 2 is the energy calculated with the MCY potential.
(within the deviations cited above), the two functions differed in other regions of the potential energy surface, and these differences led to qualitatively different results when applied to liquid water. To overcome this difficulty, the parameterization was carried out by using more dimer configurations, selected particularly from the repulsive regions of the energy surface. Instead of generating the dimer energies for these new configurations by ab initio calculations, they were obtained directly from the empirical MCY potential. Sixteen additional dimer configurations of the type shown in Figure 2 were added to the data set, and the square of the deviations of the energies of the 82 dimer configurations and of the dipole moment was reminimized as a function of the nine adjustable EPEN parameters. A least-squares deviation of 0.536 was obtained primarily due to the inclusion of repulsive configurations; the final parameters are given in Table I.
2506
The Journal of Physical Chemistry, Vol. 82, No. 23, 1978
R. A. Nemenoff, J. Snir, and H. A. Scheraga
TABLE I: Revised (EPEN/SI Parameters for Water parameter
value
distance of 0-Helectron pair from 0 nucleus distance of lone pair from 0 nucleus lone pair-0-lone pair angle A (bonding electron) B (bonding electron) C (bonding electron) A (lone-pair electron) B (lone-pair electron) C (lone-pair electron)
0.56716 A 0.27207 W 102.0" 697.20 7 kcal/mol 3.24826 A d ' 51.609 kcal/mol A 6 18294.604 kcal/mol 5.82272 A - '
dipole moment for water monomer
+'H
0.000 kcal/mol A 6 expt
1.85 D',b
EPEN/2 2.249
D
( +I H
Figure 4. EPENIP model for saturated hydrocarbon fragment.
a Reference 16. The dipole moment was given low weighting during the fitting procedure.
P
Monte Carlo calculations in the NPT ensemble were carried out again, using 100 water molecules, with periodic boundary conditions. The results were greatly improved over those obtained with the original EPEN potential or with the modified EPEN potential (parameterized on only 66 dimer configurations). The new potential (EPEN/2) gives good agreement with experiment for several thermodynamic properties such as the internal energy and the molar volume (see Table 11). However, some of the structural properties of liquid water are not reproduced well. The second maximum of the 0-0 radial distribution function is located at too great an 0-0 distance, and the average number of nearest neighbors of a water molecule is too high. These two results indicate that the new EPEN model of water does not possess enough tetrahedral order (leading to too low a heat capacity). The O-.O radial distribution function is similar to the one obtained with an early version of the MCY potential, parameterized on Hartree-Fock energies without any configuration interaction or dispersion ~0rrection.l~ Further study is required to understand why the new EPEN does not reproduce these structural properties, whereas the MCY potential does, although the two models give close agreement for most regions of the water dimer potential surface. Although larger deviations of the EPEN/2 energy from the MCY energy occur in strongly repulsive regions (greater than 3 kcal/mol), such configurations do not make a significant contribution to the Monte Carlo results.1° It appears that the simple fitting of a function to the 82 dimer energies (and the experimental dipole moment) does not guarantee that the potential will reproduce the properties of liquid water well.
,, , \
-i
X
t1
Figure 5. EPEN/2 model for amines. P represents the sum of unit vectors N-H,, N-H,, and N-X. The lone pair lies along N-P.
Reparameterization with a larger number of configurations may be required. It is our opinion that the present EPEN model is capable of being reparameterized to reproduce these properties. Conformational properties of molecules containing a saturated oxygen fragment, and crystal structures of molecules containing oxygen, are not sensitive to small changes in the parameters for water. Therefore, changes in the parameters for water will not affect any of the other parameterizations substantially. This point is discussed in paper 4.'
IV. Tests of the Potential The parameterized revised potential (EPEN/2) was tested on a variety of saturated systems. These tests were made on the properties of both single molecules, and on systems containing more than one molecule. A. Single Molecule Tests. Conformational studies on a variety of saturated organic molecules were carried out. The conformational space was searched for zero-gradient
-
TABLE 11: Properties of Liquid Water' orooertv
exotlb value
(U, ) / N , dkcal/mol
( E )/N, kcal/ mol C p , p ~ ~cal/mol ,f deg C p , cal/mol deg ( W I N . cm3/mo1 ~ , 16-6/atm h no. of nearest neighbors position of 2nd i a x in go ...o,A
-8.14 18.0 18.07 41 -5.2 4.6
-
Monte Carlo valuesC EPEN/2 MCY -10.26 i 0.18 -8.3 i: 0.3e 6.5 r 1
10.1 i I f 16.44 i 0.05 26 i 4 6.97 i 0.03 -6.0
-
-9.07 t 0.13 - 7 . 1 i 0.3e 17 i 7 20.6 c 7g 23.81 i 0.16 47 i 16 4.05 i 0.03 -4.8
References to a All quantities refer to T = 298 K, P = 1 atm. See ref 1 3 for definitions of thermodynamic quantities. experimental results are given in ref 13. Both Monte Carlo calculations were carried out for N = 100 water molecules, in This is the value of the potential energy obtained the NPT ensemble. Results for the MCY calculation are from ref 13. in the MC calculation, e Calculated using estimates of the kinetic energy and quantum corrections discussed in Appendix to ref 13. f This is the contribution to the heat capacity from fluctuations in U,. g Calculated using estimate of vibrational heat capacity as discussed in Appendix to ref 13. This is the isothermal compressibility.
The Journal of Physical Chemistry, Vol. 82,No. 23, 1978 2507
H,O and Saturated Organic Molecules
TABLE VII: Comparison of Experimental and Calculated Results on Free Molecules Used to Test EPEN/2 molecule property experiment EPEN/2 isobutanea neopentaned n-butane ethylamine'
isopropylamine!
CH, barrier, kcal/mol dipole moment, D CH, barrier, kcal/mol trans-gauche barrier, kcal/mol gauche-gauche barrier, kcal/mol CH, barrier in g conform, kcal/mol CH, barrier in t conform, kcal/mol NH, barrier (g + g), kcal/mol NH, barrier (g + t), kcal/mol AE( g-t), kcal/mol P,
3.9 i 0.75b 0.13c 4.3-4.7e 3.6-4.2g 6.1,h 7.4* 3.743 1.6' 2.1' -0.297' 1.23k 4.4,'" 2.7'" 2.1'" 2.1m -0.12n 5.0,'" 4.1,m 2.7'" 2.5'" 1.29-1.32' 1.09q 0.98' 1.38' 1.693r 3.300' 2.702' 0.146' 1.165' 1.174' 1.5' 55.4 c 2t 28.3 c 2t 159.5 f I t 3.05t
D
CH, barrier in g conform, kcal/mol NH, barrier (g t), kcal/mol NH, barrier ( t -t t), kcal/mol A E ( g-t), kcal/mol CH, barrier, kcal/mol NH, barrier, kcal/mol -j.
tert-butylamine methanolp
P, D
CH, barrier, kcal/mol ( i i to C-0), D p Y (Ito C-0), D /.tx
methyl ethyl ether'
Pt,
D
Ma,
D
CH,-C barrier, kcal/mol CH,-0 barrier, kcal/mol
Pb, D Ptr
4.282 0.04 4.579 4.039 12.80 4.59 4.53 1.94 2.86 -0.306 1.470 5.12 2.57 3.49 -0.264 5.69 2.93 1.47 0.988 0.600 1.865 1.960 3.912 2.646 0.298 1.575 1.603 2.180 54.58 38.80 167.47 3.34
D
AE(g-t), kcal/mol 2-aminoe thanolt (0CC)-(CCN) dihedral angle, deg (CC0)-(COH) dihedral angle, deg (CCN)-(CNH) dihedral angle, deg P, D a Geometry taken from ref 34. Reference 35. Reference 36.. Geometry taken from ref 37. e Reference 38. f Geometry taken from ref 24. g Reference 39. Reference 40. Geometry and energy differences taken from ref 41. Reference 42. Reference 43. Geometry used was the same as for ethylamine, with one or both of the methylene hydrogens replaced by a methyl group. Reference 44. Reference 45. O Reference 46. Geometry taken from ref 29. Geometry and dipole Reference 48. ' Geometry and other experimental quantities taken from ref 49. Reference 47. moment taken from ref 50.
'
points using a modified version of a Newton-Raphson minimizer,33with all dihedral angles varied simultaneously. The convention used for defining the dihedral angle is the same as described in the original EPEN work.18 A positive value of the dihedral angle corresponds to rotating the distal functional group clockwise relative to the near group. For hydrocarbons, the zero value of dihedral angle corresponds to an eclipsed conformation (see Figure 6a). In alcohols or ethers the zero value refers to the 0-H or 0-C bond, respectively, eclipsing one of the C-H bonds in a methyl group, or the C-C bond in a methylene group (see Figure 6b). In amines, the lone pair of the nitrogen eclipsing the C-H bond for methyl groups, and the C-C bond for methylene groups, defines the zero of dihedral angle (see Figure 6c). Gauche refers to a value of the dihedral angle of approximately 60°. Experimental geometries were used wherever possible. A comparison of the EPEN/2 results with experimental values is presented in Table VII. The results are discussed below. (1)Alkanes. The potential was parameterized on the trans-gauche energy difference in butane. As a test, the trans-gauche barrier and the gauche-gauche barrier for this molecule were calculated. The trans-gauche barrier occurs at a value of 1 1 9 . 9 O for the central C-C dihedral angle, and the calculated barrier height is in good agreement with e ~ p e r i m e n t . For ~ ~ the gauche-gauche barrier, the calculated value is too high.3gv40The position of the gauche-gauche barrier has the two terminal methyl groups in close contact. It might be expected that the rigid-geometry approximation used in EPEN/2 would break down at such close contacts, and bond angle distortions would occur.
a
c
H H
Flgure 6. Newman proJectionsillustrating convention for definition of dihedral angles. Positive dihedral angle corresponds to clockwise rotation of the distal group relative to the near group. Dotted lines correspond to lone pairs: (a) rotation of two methylene groups, 4 = 0 when C eclipses C; (b) rotation of hydroxyl group, 4 = 0 when H eclipses C; (c) rotation of amine group, 4 = 0 when lone pair eclipses C.
The methyl rotational barriers calculated for isobutane and neopentane are in good agreement with experiment.36J8 (2) Amines. Conformational studies were carried out on ethylamine, isopropylamine, and tert-butylamine. For
2508
The Journal of Physical Chemistry, Vol. 82, No. 23, 1978
.e.
