J . Phys. Chem. 1990, 94, 1656-1660
1656
velocity autocorrelation function, defined as A ( t ) = (v'(O)*v'(t))/(v'')
(20)
This function and its Fourier cosine transform are shown in Figure 13. The initial rapid decay of the correlation function is followed by oscillatory behavior on the time scale of tenths of picoseconds. Analyzing the density of states, we see an initial sharp peak at around 60 cm-I with a broad shoulder centered at 175 cm-'. These values are in excellent agreement with experimental and correspond to the bending and stretching motions of the hydrogen bond (0-H-0) unit.
Conclusions The nonempirical NEMO potential construct is convenient and fast compared to direct S C F calculations. Inclusion of specific potential terms representing physically meaningful energy terms minimizes the actual number of parameters in a given potential fit as well as increases in the interpretation of each contribution. Many-body polarization interactions are important for successfully using the pair potential in a liquid simulation. This property will also be of great interest when using NEMO water as solvent for
*
biomolecules and other solutes since the model will be able to respond nonlinearly to changes in its environment. The NEMO liquid water properties are shown to agree to a high degree with experimental data and other common water potentials. There is approximately 0.5 kcal/mol lacking in internal energy, but this could be accounted for by the rigidity of the model. The dynamics is slower than that measured experimentally, but the type of hindered translational motion characteristic for hydrogen-bonded water is well reproduced. The total induced dipole moment is high, about 0.8 D, and can be traced back to the Hartree-Fock S C F approximation. It is clear that NEMO water is capable of reproducing liquid water properties to the level of other water model potentials. This implies that a nonempirical force field of reasonable accuracy and consistency could be constructed from the NEMO procedure. The added reliability of the entire force field warrants the extra computational effort. Acknowledgment. We thank Drs. B. Jonsson and 0.Teleman for a careful reading of the manuscript. Registry No. H20, 7732-18-5.
Analytical Intermolecular Potential Functions from ab Initio Self-Consistent-Field Calculations for Hydration of Methylamine and Methylammonlum Ion Sung So01 Wee, Seungmoak Kim, Mu Shik Jhon, Department of Chemistry, Korea Advanced Institute of Science and Technology, P.O. Box 150. Chong-yang Ni, Seoul, 131 -650 Korea
and Harold A. Scheraga* Baker Laboratory of Chemistry, Cornell University, Ithaca. New York 14853-1301 (Received: May I , 1989; In Final Form: August 4, 1989)
The interaction energies between water (H20)and methylamine (CH,NH,) and methylammonium ion (CH3NH3'), respectively, are of interest because these are prototypes of side-chain groups of proteins. Such intermolecular potential functions are necessary for Monte Carlo calculations to obtain free energies of hydration. The solute-water interaction energies were derived from ab initio SCF calculations, using a 6-31G basis set. The monomer geometries of CH3NH2,CH3NH3+,and H 2 0 were kept rigid, with the bond lengths and bond angles taken from experiment. The computed ab initio SCF energies of the solute-water complexes were used to obtain simple analytic potential functions in the form of inverse powers of interatomic distances. The structures of the solute-water complexes were optimized with the computed potential functions and compare very well with those calculated with the 6-31G basis set; both contain nearly linear hydrogen bonds, with H 2 0 being a proton donor to CH3NH2but a proton acceptor from CH3NH3+.
1. Introduction
TABLE I:
This is a continuation of our previous paper' to obtain intermolecular potential functions from a b initio S C F calculations for use in Monte Carlo calculations to compute bulk properties. In this work, our initial efforts are focused on model compounds that represent the side-chain groups of proteins. Our previous paper' was concerned with the hydration of CH4, CH30H, CH3COOH, and CH,COO-, and we now treat CH3NH2and CH3NH3+in this paper. For this purpose, we carried out a b initio S C F calculations of the potential energy surface of the CH3NH2-.H20 and CH3NH3+-.H20complexes with the 6-31G basis set and the restricted Hartree-Fock method. It should be noted that the potentials depend on the quality of the basis set and on the method for correcting for the correlation energy in the ab initio S C F calculations. A minimal basis set such as STO-3G or 3-21G exaggerates the anisotropy of the charge distribution (hence, overestimating the interaction energy) and leads to a poor rep( I ) Kim. s.;Jhon, M. s.;Scheraga, H. A. J. Phys. Chem. 1988, 92, 7216.
0022-3654/90/2094- 1656$02.50/0
Geometries of CH,NHh CHqNHIt, and H 2 0
bond lengths, A
R(O-H) 0.9572
angles, deg
H20
LHOH 104.52
CH3NH2and CH3NH3+
R(C-N) 1.47 R(N-H) 1.03 R(C-H) 1.11
LHNH 106.0 LCNH 11 1.5 LHCH 108.4 LHCN 110.5
resentation of the hydrogen-bonding structure. Clementi and co-workers2 have developed potentials for the water dimer with a large polarization basis set and also used configuration interaction (CI), to compute the correlation energy. However, to reduce the required amount of computer time, the monomer (2) Popkie, H.; Kistenmacher, H.; Clementi, E. J. Chem. Phys. 1973, 59, 1325. (3) Matsuoka, 0.;Clementi, E.; Yoshimine, M. J. Chem. Phys. 1976,64, 1351.
0 1990 American Chemical Society
S C F Calculations for Hydration of CH3NH2and CH3NH3+
The Journal of Physical Chemistry, Vol. 94, No. 4, 1990 1657
TABLE 11: Cartesian Coordinates' and Chargesbof Solute Molecules atom X Y Z 9
solute water interaction energies were calculated with the 6-31G basis set with the Gaussian 82 program." The solute molecules were placed at fixed points in space with the nitrogen atoms of CH3NH2 and CH3NH3+taken as the origin of the Cartesiancoordinate system. For each specified direction and orientation, the relative position of the water molecule was varied with respect to the fixed solute molecules to calculate the interaction energies as a function of the Ne-0 distance (RN+). Also, the symmetries of the CH3NH2 and CH3NH3+ molecules were used. The CH3NH2(staggered form) has u8 symmetry, and thus only half of the CH3NH2 molecule was considered; for the CH3NH3+ molecule (staggered form) with DW symmetry, only one N-H and C-H direction of the CH3NH3+molecule was considered since the NH3+ and CH3 groups have C3, symmetry and the corresponding energy surfaces have the same symmetry. In this manner, 335 configurations for CH3NH2 H 2 0 and 337 configurations for CH3NH3+ H 2 0 were considered. All computations were carried out on the IBM 3083/JX3 system a t the Korea Advanced Institute of Science and Technology. Second, the ab initio S C F energies were fit to a potential function of the form
CHSNHZ
N
0.000
C HI(C) H2(C) H3(C) H4(N) H5(N)
0.000 -1.041 0.521 0.521 -0.492 -0.492
0.000 0.000 0.000 0.902 -0.902 -0.823 0.823
0.000 -1.472 -1.862 -1.862 -1.862 0.378 0.378
-0.808 -0.259 0.126 0.165 0.165 0.305 0.305
N C Hl(C) HVC) H3(C) H4(N) H5(N) H6(N)
0.000 0.000 -1.041 0.521 0.521 0.960 -0.480 -0.480
CHoNH3' 0.000 0.000 0.000 0.902 -0.902 0.000 -0.832 0.832
0.000 -1.470 -1.860 -1.860 -1.860 0.378 0.378 0.378
-0.871 -0.343 0.267 0.267 0.267 0.471 0.471 0.471
'In angstroms.
In electronic charge units (ecu).
geometries of the solute and water were kept rigid.2 The main purpose of the present work is to obtain an accurate description of the potential functions for these solute water systems. It is clearly desirable to have simple and effective intermolecular potential functions for the complexes. We present analytic expressions that reproduce the ab initio potential energy hypersurfaces for CH3NH2+ H 2 0 and CH3NH3++ H 2 0 .
+
+
+
+
2. Method of Calculation 2.1. Model Compounds. In Table I, the monomer geometries of the CH3NH2,CH3NH3+,and H 2 0 molecules are The geometry of CH3NH2was taken from the joint analysis of microwave spectra and electron diffraction data of Iijima et aL4 The conformations of CH3NH2and CH3NH3+were taken only as the staggered form, since the staggered form of both molecules is the most stable at room temperature.1° The Cartesian coordinates and charges of the CH3NH2 and CH3NH3+ molecules are represented in Table 11, where the Cartesian coordinates were computed from the experimental geometries and the charges were taken from a Mulliken population analysis with the 6-31G basis set." The coordinates of the H 2 0 molecule were based on the model of Clementi and c o - ~ o r k e r s , ~ and the charges were taken from a Mulliken population analysis with the 6-31G basis set; these coordinates and charges are given in Table I of our previous paper.' The ab initio S C F energies of CH3NH2, CH3NH3+, and H 2 0 with the 6-31G basis set are -95.164 922, -95.533 844, and -75.983 998 au, respectively. 2.2. Calculation Procedure. The details of the calculation procedure can be summarized in two consecutive steps; the first is an ab initio S C F calculation for the solute water complex, and the second is an optimization of the parameters of the potential functions to reproduce the ab initio S C F potential hypersurface of the complexes. First, we computed the interaction energies of the solute molecules with H 2 0 placed a t a sufficiently large number of positions and orientations relative to the solute molecules so as to provide a reasonable sampling of the potential energy surface representing the interactions of the solute water complex. The
+
+
(4) lijima, T.; Jimbo, H.; Taguchi, M. J . Mol. Strucr. 1986, 144, 381. (5) The NH, group of the CH3NH3+ion is tetrahedrally coordinated. (6) Benedict, W. S.; Gailar, N.; Plyler, E. K. J . Chem. Phys. 1956, 24, 1139. (7) Takagi, K.; Kojima, T. J . Phys. SOC.Jpn. 1971, 30, 1145. (8) Higginbotham, H. K.; Bartell, L. S. J. Chem. Phys. 1965, 42, 1131. (9) Kuchitsu, K.; Cyvin, S. J. Molecular Srrucrure and Vibration; Elsevier: Amsterdam, 1972; Chapter 12. (IO) The staggered forms of CH3NH2and CH3NH3+ions are 2.7 and 2.8 kcal/mol more stable than the eclipsed forms, respectively, with the 6-31G basis set. (,I 1). Binkley, J. S . ; Frisch, M. J.; DeForees, D. J.; Raghavachari, K.; Whiteside, R. A.; Schelgel, H. B.; Fluder, E. M.; Pople, J. A. New Version of Gaussian Program for IBMIMVS System; Department of Chemistry, Carnegie-Mellon University: Pittsburgh, PA, 1985.
where the i's are the atomic centers on H 2 0 and the j's are the atomic centers on the solutes. A, and Blj are the fitted parameters; qi and qj are the computed net charges for the i a n d j atoms; and r, is the distance between the pair of atoms i and j . The interaction energy V(SCF) is EscF(r) - ESCF(SO), expressed in kilocalories/mole, where &cF(r) is the interaction energy of the complex and EscF(50)is the interaction energy at RN+ = 50 A. The SCF energies of the complexes at R N a = 50 A are -1 7 1.148'912 au for CH3NH2+ H 2 0 and -171.517957 au for CH3NH3++ HzO, while the sums of each of the monomer energies are -171.148 919 au for CH3NH2+ H 2 0 and -171.517841 au for CH3NH3++ H2O. The parameters were optimized with the aid of a standard nonlinear least-squares method.l2,l3 The data were not weighted equally in order to favor the lower energy ones, which were more accurate than the high-energy ones. The weighting was carried out with the following equations N
6 = Cwi[V(SCF)j - V(12-4-1)J2
(2)
i=l
wi = 1
+ CY exp[-(V(SCF),
- V(SCF)&/kT]
(3) where V(SCF), denotes the a b initio SCF energy; V(12-4-l)i is the energy from the potential function in eq 1; i is the number of configurations between the solute and a water molecule; V(SCF), is the lowest ab initio SCF energy, which is -8.42 kcal/mol for CH3NH2 H 2 0 and -24.5 kcal/mol for CH3NH3+ H 2 0 , respectively; k is the Boltzmann constant; and T i s the absolute temperature (298 K). The value of a was 150, which was found by adjustment to give a good low-energy fit and still provide reasonable results for the repulsive values. The fit was made to 6-31G energies below 10 kcal/mol for CH3NH2 H 2 0 and CH3NH3+ H20. The parameters of the potential were obtained by minimizing the sum of the squares of the deviations (6) between the ab initio energy and the potential energy.
+
+
+
+
3. Results and Discussion T h e optimized parameters for the intermolecular potential functions are given in Table 111. No physical significance should be attributed to them since the ab initio SCF energies have simply been expressed in terms of inverse powers of r in the potential functions. The resulting potential functions, however, provide a very simple way to determine (with good accuracy) the positions (12) Powell, M. J. D. In Numerical Methods f o r Nonlinear Alqebraic Equations; Rabinowitz, P., Ed.; Gordon and Breach: London, 1970. (13) Dennis, J. E., Jr.; Gay, D. M.; Welsch, R. E. ACM Trans. Math. Software 1981, 7, 348, 369.
Wee et al.
1658 The Journal of Physical Chemistry, Vol. 94, No. 4, 1990 TABLE 111: Parameters for the Potential Function atom Dair
N-0
c-0
H(C)-0 H(N)-0 N-H(0) C-H(0) H(C)-H(O) H(N)-H(O)
+ H20
1.033 X 1.165 X 5.477 5.011 X 2.050 X 8.890 X 1.878 X 2.082 X
CH3NH3' N-0
c-0
H(C)-0 H(N)-0 N-H(0) C-H(0) H(C)-H(O) H(N)-H(O)
a
I n kcal Ai2/mol.
8b
A"
CH'NH2
IO6 IO6
IO3 IO4 10'
4.107 1.987 5.596 6.471 6.100 6.156 1.097 4.859
lo2
X X 10' X X lo-' X X lo-' X 10'
+ H20 IO' IOs
3.657 X 6.129 X 0.000 1.599 X 1.712 x 1.576 x 5.849 x 3.461 X
IO-! 105 105
103 IO4
1.270 X 1.295 X 1.213 x 1.015 X 1.085 X 1.756 X 4.088 X 1.290 X
IO2 IO2 10-4
IO-! IO2 IO2 IO1 lo2
-I
kcal A4/mol.
and orientations corresponding to the minimum-energy configurations of the H 2 0 molecule around the CH3NH2and CH3NH3+ molecules. The potential functions should provide reasonable energy hypersurfaces for the complexes. The standard deviations for the fit are u < ~ 1.04; , U
