3084
J. Phys. Chem. 1995, 99, 3084-3092
Molecular Solvation Potential. A New Tool for the Quantum Mechanical Description of Hydration in Organic and Bioorganic Molecules C. Alhambra and F. J. Luque* Departament de F a d c i a , Unitat Fisicoquimica, Facultat de Farmhcia, Avgda Diagonal sn, Barcelona 08028, Spain
Modesto Orozco* Departament de Bioquimica Facultat de Quimica, Universitat de Barcelona, Marti i FranquLs I , Barcelona 08028, Spain Received: June 17, 1994; In Final Form: November 3, 1994@
A new method for the representation of interactions between solutes and their hydration waters is presented. The algorithm is based on a generalization of the molecular interaction potential. The method, after a careful parametrization of the van der Waals interactions between the solute and the waters, provides a very accurate quasi-quantum mechanical picture of the hydration characteristics of the solute. Some applications of the method to the study of hydration patterns of relevant molecules are reported.
Introduction The molecular electrostaticpotential (MEP) provides the exact value for the electrostatic component of the interaction energy between an unperturbed charge distribution and a proton. Furthermore, the MEP also gives a value general picture of the reactive properties of the molecules with a reasonable computational cost. This explains its large number of applications in very different fields of physics, chemistry, biochemistry, and pharmacology.’ The MEP has been recently generalized as shown in eq 1. The new potential, named molecular interaction potential (MIP; 2 ) , includes an empirical van der Waals term that makes it possible to avoid the nuclear discontinuity of the MEP, enabling the study of interactions with both negative and positive fractional charges. The MIP is not a rigorously defied parameter like the MEP but can provide a more realistic picture of the nonbonded molecular interactions than the MEP.*v3 Therefore, it can be regarded not as a substitute of the MEP but as a useful complement for a priori estimates of molecular
interaction^.^ V,,(r,) =
q; A
Y1A
+
1A
c A
ZAQI Irl
-
-RAI
(3) where c stands for the atomic hardness and R , for the van der Waals radius. The MIP provides a representation of the noncovalent interactions between the molecule and a monopole. This is very informative in some processes, but a finer description of the reactive features would be desirable, especially when the interaction involves not a point monopole but a neutral molecule, @
Abstract published in Advance ACS Abstracts, February 1, 1995.
for which the electrostatic interaction is determined from the dipole or higher order multipoles. Accordingly, it seemed suitable to modify the MIP to properly account for these interactions between the “quantum” molecule and polyatomic classical groups. The interaction with a water molecule is especially interesting. The modified potential, whose development is reported here, is named molecular solvation potential (MSP). A detailed microscopic understanding of the solvation phenomena can be obtained from classical Monte Carlo (MC) and molecular dynamics (MD)4 or quantum mechanicslmolecular mechanics (QM/MM) method^.^ The former methods generally use an atom-centered Coulombic expression for the electrostatic interactions, which often cannot provide sufficient accuracy.6 On the other hand, the latter methods can be very expensive from a computational point of view if MD or MC representations of the solvent are used, even when the quantum mechanical subsystem is treated at the semiempirical level. In this context, the MSP was conceived as a compromise between accuracy and expensiveness. The calculation of the MSP requires the previous resolution of the Schrodinger equation in gas phase (whatever the method used in the computation, e.g., semiempirical or ab initio, onedeterminantal SCF, many-determinantal CI, etc.), or even incorporating environment effects (as in SCRF7 computations). From the molecular wave function the three-center integrals in eq 1 can be rigorously evaluated and the MSP be determined for any solute-water arrangement. Then, the quality of the MSP will mainly depend on the solute wave function and an the classical water modeL8 In this paper we present the development of the MSP and its refinement from SCF calculations, as well as its application at both ab initio and semiempirical levels to the study of some relevant systems.
Methods The TIP3P model8 of water was used to evaluate the MSP due to (i) its proven ability to reproduce properties of pure watefi and (ii) the similarity between TIP3P and ab initio 6-31G* electrostatic charges which guarantees a good coupling between TIP3P and 6-31G* charge densities. Accordingly, the MIP (eq 1) was slightly modified to account for the interaction of the
0022-365419512099-3084$09.0010 0 1995 American Chemical Society
J. Phys. Chem., Vol. 99, No. IO, 1995 3085
Molecular Solvation Potential quantum mechanical solute with this classical three-center water (eq 4). Since the probe molecule in the MSP is a three-center water (TP3P), the value of the MSP at point rl (the position of the oxygen) depends on the orientation of the water molecule, defined by three rotation angles. Accordingly, to define the MSP unequivocally, the orientation of the water molecule needs to be optimized. This is done by minimizing the MSP with respect to external rotations at rl using a gradient method. In order to accelerate this process, the relative position of the water was initially optimized by using a classical potential function, where a Coulombic expression replaces the third term on the right-hand side of eq 4. Electrostatic chargesg were used to represent the solute charge distribution in this preliminary calculation. After this preoptimization, the position of the water was further refined using the rigorous expression of the MSP in eq 4.
SCHEME 1
I/
H3C"2(n)
HQ"z(TI)
a
H,COH(TI)
TABLE 1: Polarization ( E 1) and al, Initio HF/SCF (ESCF) Interaction Energies (in kc$rnol) between Different Molecules and a Classical TIP3P Water"&
where the index l i refers to the position of the oxygen and hydrogens of the water molecule. The T P 3 P model assigns charges of -0.834 and +0.417 e to the oxygen and hydrogen, respectively. The combination rules in eqs 3 and 4 were used in the calculation of the van der Waals term. As in the original TIP3P model, all the van der Waals interactions were assigned to the oxygen (R, = 1.768 A and E = 0.1521 kcaymol). The solute van der Waals parameters were determined by fitting the MSP to SCF results derived from ab initio SCF 6-31G*1° computations, where both solute and water were treated quantum mechanically. In these calculations the basis set superposition error (BSSE), which was small in all the cases, was corrected using Boys's method." The MSP was determined at the ab initio 6-31G*1° and AM1 l2 semiempirical levels. Standard methodsla were used to compute the three-center integrals in eq 4 for ab initio calculations, but at the semiempirical level the NDDO strategy13 was used. The selection of the AM1 Hamiltonian was due to previous studies13bwhere we demonstrated the superiority of the AMl/NDDO strategy with respect to other semiempirical methods in order to obtain the MEP. The MSP was computed using our MOPETE program.I4 Ab initio calculations were performed using HONDO-77. l5 Semiempirical calculations were carried out using a modified version of MOPAC.16 Optimization of van der Waals parameters was performed by fitting MSP and ab initio 6-31G* data with the SIMPLEXFIT ~r0gram.l~ Results and Discussion Suitability of the MSP Formalism. The MSP includes two different energy contributions: electrostatic, which is rigorously computed from the MEP, and dispersion-repulsion, which is approximated by a r-12-r--6 relationship. Therefore, the main source of systematic errors is expected to arise from an incomplete representation of the polarization between the solute and the solvent molecules. This shortcoming may be partially solved by using a solvent-adapted wave function within the SCRF framework,' which includes the polarizing effect of the bulk solvation waters. However, even in this case the specific polarizing effect of the water molecule directly bound to the solute is still neglected in the MSP.
molecule Hz0-Hz0 3" (TI) 3"
('W
CH4 HzCO (Tl) HzCO (T2) HzCO (T3) CH3"2 (TI) C " f z (T2) CH30H (Tl) CH3OH (T2)
EFml -0.30 -0.28 -0.11 -0.05 -0.25
ESCP
-0.19 -0.31
-6.12 -5.13 -2.00 -0.12 -1.85 -2.40 -4.04 - 1.83
-0.02
-5.78
-0.35 -0.23
-4.29 -4.22
-0.44
r ('4 2.16 3.04 3.31 3.15 3.63 3.21 2.96 3.29 3.08 3.53 2.91
For structures, see Scheme 1. Geometries were determined at the 6-31G* level (see text). Distances are in angstroms. To elucidate whether or not the neglect of such a specific polarization term was justifiable, the interaction energy with a TIP3P water was computed at the SCF level using the 6-31G* basis set for a series of prototypical molecules, whose geometry was previously optimized at this computational level. After optimizing the relative position of the T P 3 P water with regard to the molecule, the polarization was estimated applying eq 5, where V stands for the classical potential generated by the charge distribution of the TIP3P water, Yois the unperturbed molecular wave function, which was used as a guess in the SCF procedure, and Y is the final wave function of the solute in the presence of the TIP3P water molecule.
Epl = (YIH
+ VlY) - (VPIH + Vlwo)
(5)
The values obtained for this series of molecules and for different orientations of the TIP3P water are displayed in Table 1. The magnitude of Epol is small in all the cases, with a maximum (in absolute values) of -0.42 kcaUmol for the HzCO and an average value of only -0.25 kcdmol for the 11 interactions. These values are only a small fraction of the total interaction energy (near 7% in average) at the 6-31G* level. It might be expected that the magnitude of this polarization effect is larger if more flexible wave functions are used. However, it seems clear that the specific polarization effect of a classical T P 3 P water directly bound to the solute is not very large, and the neglect of such an explicit term in the MSP potential function
J. Phys. Chem., Vol. 99, No. 10, 1995
Alhambra et al. TABLE 2: van der Waals Parameters (E, kcal/mol; R , h;; See Eq 4) for MSP Calculation& atom E R atom E R C 0.0013 3.263 0, 0.0053 2.360 0 0.0121 2.319 CA 0.0002 3.037 N 0.0181 2.488 CAI 0.0027 2.948 CC 0.0001 3.009 NA 0.0652 2.244 0, 0.0010 2.746 H 0.0000 0.000 N, 0.0219 2.199
4
?ORYIC Fa
I
t
For structures, see Scheme 2. The relative orientation of the dimers used to derive the parameters are shown.
9
2'
3-
4-
5-
-6'
7-
8-
9
2.0
2.5
3.0
35
R(W-0)
4.0
45
5.0
(A)
4
t
PYRINIDIHE Tl
2-
0'
2-
4-
6-
Figure 1. Comparison of the ab initio HF/SCF energy profiles corresponding to the approach of a water to the solute molecule and the profiles obtained by using the OPLS and AMBER parameters in the van der Waals tern. Solute molecules are formic acid, formamide, and pyrimidine. is not a source of large errors in MSP calculations. Furthermore, it is worth noting for our purposes that this polarization effect might be captured indirectly in the van der Waals term during the parametrization process (see below). Parametrization of the van der Waals Term. The MSP includes a van der Waals term that should account for all the nonelectrostatic interactions between the solute and the water. When combination rules are used, the calculation of the van der Waals contribution only requires the assignment of the van der Waals Parameters for the atoms of the solute. This may appear a priori a trivial task, since they could be taken from existing force fields. This simple strategy was adopted in the development of the MIP2 as well as in other standard QM/MM
computation^.^ However, it must be stressed that the transferability of van der Waals parameters from empirical force fields to mixed classical-quantum mechanical calculations is not demonstrated but needs to be assessed. To test the transferability of the standard van der Waals parameters to the MSP strategy, the energy profiles corresponding to the approach of a water to different molecules were determined at the ab initio 6-31G* SCF level. Those profiles were then compared with energy values derived from the MSP, in which either AMBER" or OPLS18,19standard van der Waals parameters were considered. Figure 1 shows selected results for formic acid, formamide, and pyrimidine. The results clearly state that the use of standard parameters underestimates the optimum solute-water distance, which becomes shorter than the SCF distance by more than 0.3 A. This leads to errors of several kcaUmol in the interaction energy in the inner regions. Although the errors are smaller at larger distances, the use of van der Waals parameters derived from classical force fields should be avoided in mixed quantum mechanics/molecular mechanics calculations. The reason for this failure lies in the differences in strength of the electrostatics at short distances between classical-classical or quantum mechanical-classical molecular interactions. Classical van der Waals parameters avoid a close contact between positive and negative atoms, whose charges are typically in the range of f l e. Nevertheless, they are not adequate to modulate the approach of a classical negative atom to the large positive charge generated by the nuclei of a molecule treated quantum mechanically. There are several possible ways to solve the problem of the transferability of classical van der Waals parameters to the MSP. We tried to scale the standard van der Waals parameters, but we were unable to obtain reliable results. Another option was the addition of an extra overlap term in the MSP to represent more accurately the repulsion between the solute and the classical water. However, after several pilot calculations we ruled out the inclusion of additional terms and decided to reoptimize fully the van der Waals parameters. For this purpose the MSP was fitted to the 6-31G* SCF energy profiles for the approach of a quantum mechanical water to the molecule. Each profile was defined by means of 7-10 points, and several orientationsz0of the water were considered for the approach to the molecules used in the parametrization: H20 (4),HCOOH (3), H2CO (2),HCONH2 (4),CHq (l), "3 ( 2 ) , and pyrimidine (2). At every point the electrostatic part of the MSP was subtracted from the SCF energy, and the residual energy was fitted to the van der Waals term. Table 2 shows the fitted van der Waals parameters. The optimized parameters clearly differ from those reported in classical force fields. The van der Waals radii (&) are larger, and the hardness parameters (c) are smaller. This means that bigger, but more deformable, atoms should be defined in mixed simulations. The only concern is the hardness for sp2 carbons (Cc, CA, CAI), which is especially small. This could lead to errors in the representation of nonbonded interactions
J. Phys. Chem., Vol. 99, No. 10, 1995 3087
Molecular Solvation Potential
TABLE 3: Optical Interaction Energies (kcaYmo1) and Equilibrium Distances (A) for All the Profiles Used in the Fitting and for Five Test Dimers! molecule ESCF EMSP RSCF RMSP HzO
“3
CHq HzCO HCOOH
x
h Reg. coef.= 0.96 Corr. coef.= 0 . 9 8
-12
-10
-8
-6
-4
-2
0
2
4
6
HCONHz 8
IO
WSP (Xcal/mol) Figure 2. Ab initio HF/SCF and MSP energy profiles for selected molecule orientations. Fitted van der Waals parameters shown in Table 2 were used in evaluating the dispersion-repulsion term of the MSP.
pyrimidine
SCHEME 2
CH3OCH3 CH3CONH2
II
I 0
‘
CH3OH CH3NHZ
NHiCONHz
II,
pyridine
T1 T2 T3 T4 T1 T2 T2 T3 T1 T2 T3 T1 T2 T3 T4 T1 T2 T1 T2 T1 T2 T1 T2 T1 T2
-4.65 -4.70 -4.22 -3.23 -5.56 -2.03 -0.17 -2.28 -3.76 -8.80 -3.91 -1.32 -3.95 -4.81 -2.04 -5.27 -4.17 -0.85 -4.29 -4.22 -1.83 -5.78 -3.80 -6.61 -7.76 -7.91 -5.29 -2.82
-4.35 -4.91 -4.80 -3.39
-5.50 -2.15 -0.14 -2.23 -3.52 -7.16 -3.57 -1.53 -5.11 -5.04 -2.30 -4.95 -4.32 -1.10 -4.19 -4.67 -1.99 -5.53 -2.75 -6.29 -9.05 -9.87 -6.51 -2.70
2.93 2.94 3.05 3.12 3.09 3.38 4.16 3.25 3.09 2.89 3.11 3.23 3.09 3.11 3.22 3.01 3.12 3.57 3.53 2.97 3.29 3.08 2.94 2.95 2.89 2.87 3.10 4.07
2.99 2.94 2.95 3.13 3.11 3.30 4.22 3.21 3.03 3.05 3.10 3.05 3.12 2.93 3.08 2.99 3.06 3.44 3.54 2.96 3.32 3.08 3.03 2.96 2.95 2.87 3.09 4.12
a For CH30CH3the water is on the oxygen atom; for the two dimers (Tl, T2) of CHsCONHz water is placed between 0 and N atoms; for NHzCONHzwater is between N and 0 atoms (Tl) and on the N atom (T2); for pyridine the water is on the N atom.
HCOOH
HCONHz
CH4
11 I{ ‘0’
PYRYMIDINE
“3
between the water and sp2 centers, especially when the water is placed perpendicular to the sp2 plane. However, calculations of the interaction between a water with the n-charge distribution of benzene (see below), and with H2CO (results not shown) allowed us to verify the satisfactory behavior of these parameters to reflect both “in-plane’’ and “out-of-plane’’ interactions with water. Figure 2 stresses the quality of the optimized van der Waals parameters for some molecules. The agreement between SCF and MSP results is very good in all the cases, especially in the minima region, and is clearly better with respect to the use of empirical parameters (see Figure 1). Since the Parameters in Table 2 were derived from a limited set of molecules, their transferability to other molecules is questionable. To test this point, the SCF and MSP equilibrium
geometries and interaction energies for different dimers were determined. MSP calculations (see Table 3) reproduce equilibrium water-molecule distances with an average error less than 0.06 8, and the optimum interaction energy with an average error around 0.6 kcal/mol. The set of parameters in Table 2 is only a preliminary force field, which may be improved by fitting to a larger set of molecules. However, results in Table 3 are quite satisfactory and suggest a reasonable degree of transferability for the optimized van der Waals parameters. For instance, the transfer of van der Waals parameters of water and methane to methanol, which is a priori a rather drastic transfer for an united atom model, leads to quite reasonable results. The suitability of van der Waals parameters for transfer to larger molecules is discussed below. The assessment of the MSP to describe water-molecule interactions is provided by Figure 3, where the 6-31G* SCF and MSP energies for more than 150 interactions (for molecules in Tables 2 and 3) are compared. A mean unassigned error (rms) of 0.57 kcdmol is found for the whole set of interaction energies, which lie in the range -10 to +6 kcdmol. The regression coefficient, which is close to 0.96, shows no systematic deviation between SCF and MSP results. Finally, the correlation coefficient of 0.98 reflects the goodness of fit. Use of the MSP to Represent Hydration Propensities of Molecules. The potential use of the MSP to study hydration processes can be illustrated by means of few examples: (i) the determination of the isopotential hydration map in the molecular plane of a water, (ii) the estimation of the interaction energy between a water and a benzene, and (iii) the calculation of the MSP maps of the DNA nucleic bases. Calculation of the MSP of a Water Molecule. The MSP of a water molecule was computed at 7500 points in the molecular plane of the water using a 6-31G* wave function for
Alhambra et al.
3088 J. Phys. Chem., Vol. 99, No. 10, 1995
I N
.
.
T
T
9
L1
I
N
0
N
Molecular Solvation Potential
11
J. Phys. Chem., Vol. 99,No. IO, 1995 3089 H20
-
9 7 -
5 -
3 1 -
u
-3 -5 -
-1
-7
:
4 ti
1,,.=3.59
1,,,=5.26
E,,.=-0.76
Em,,=-1.22
6-
-
1-
-13 -15
-15 -13
-11
-9
-7
-5
-3
-1
I
3
5
7
9
11
13
15
Figure 4. MSP isopotential map of a water molecule. Energy difference (in kcal/mol) between isopotential contour lines is 1 kcal/mol. Distance is given in atomic units.
2-
0-
2-
the solute water. [This calculation took about 12 CPU hours in an IBM-3090, while the same calculation at the SCF level (with the same basis) would take up to 15 days of CPU time.] The isopotential maps in Figure 4 provide a complete picture of the hydration structure around a water. This picture exhibits many of the features reported by other authors from spatial distribution functions.21 Two identical regions of optimum interaction appear along the 0-H line, for which the TIP3P water acts as hydrogen-bond acceptor. The energy for the optimum interaction in this region is -4.9 kcdmol, which compares with the SCF value of -4.6 kcdmol. A broader area of favorable interaction appears in the oxygen region. Here the TIP3P water acts as hydrogen-bond donor, and MSP values around -4 kcdmol are reached. This region is not really a minimum, which is located outside the molecular plane (above and below) following approximately the expected direction of the lone pairs. The MSP value in the true minima is -4.9 kcal/mol, in agreement with the symmetry in HOH-OH2 and H20--HOH interactions. Calculation of the Interactions between a Water and Benzene. The ability of the MSP to reproduce the interactions with n-charge distributions was explored for benzene. The MSP was computed along the c6 axis at distances from the molecular plane to the oxygen of the TIP3P water ranging from 2.6 to 5.0 8, (T1 in Figure 5). The orientation of the water at every point was fully optimized (see Methods). The energy profile was flat, but a clear minimum of -2.10 kcal/mol was detected at a distance of 3.37 A. This value compares well with accurate experimental data,22which suggest an equilibrium distance of 3.347 A and an optimum interaction between -1.63 and -2.78 kcal/mol. It is interesting to note that the small E assigned to the aromatic carbon (see above) is enough to avoid the collapse of the water into the benzene. Finally, it is worth noting that the MSP profile demonstrates the importance of electrostatic effects in the water-benzene interactions, as shown by the low decay of the energy with the distance. (Even at 5 A the interaction energy is almost - 1 kcdmol.) Two benzene-water dimers were also examined to determine the suitability of the MSP to describe aryl interactions, which correspond to the approach of a water to benzene (i) along the line bisecting the C-C bond in the molecular plane (T2) and (ii) along a C-H bond (T3). Figure 5 shows that the minimum for the T2 profile (a bifurcated interaction) is located at 3.59 A
R
(A)
Figure 5. Three water orientations around a benzene ring studied with the MSP. Equilibrium distances (r-) are in 8, and energies (Emin)in
kcdmol. MSP profiles for these three benzene-water orientations are displayed below. from the center of the ring, the stabilization energy being -0.76 kcal/mol, whereas for the T3 minimum the distance is 5.26 A and the stabilization energy is -1.22 kcdmol. There is not direct experimental data on aryl interaction energy, but the values found are in agreement with the experimental evidence and also with previous theoretical calculation^,^^^^^ which reported energy values around - 1 k c d mol. Indirect support comes from the fact that both kinds of aryl interactions are less stable than the interaction, which agrees with experimental data.22 In order to gain further evidence of the accuracy of the MSP, we optimized the waterbenzene dimers at the 6-31G* SCF level. The optimum SCF distances (water-center of ring) for T2 and T3 were 3.6 and 5.3 A, which match the MSP values. The agreement is slightly poorer for the interaction energies (after BSSE correction), since the preference for linear interactions is not seen in the SCF results. However, the SCF values (-0.8 (T2) and -0.7 k c d mol (T3) are quite similar to the MSP values, a finding which supports the suitability of the MSP. The ability of the MSP to reproduce interactions involving benzene could appear as surprising since polarization effects are not directly included in MSP formalism. However, these results are in good agreement with previous studies,24which suggest that even a simple and non-polarizable r-1-r-6-r-12 potential can reproduce qualitatively several benzene interactions. Calculation of the MSP of the Bases of the DNA. To test the potential use of the MSP in drug design or biochemical fields, we computed the MSP of the DNA bases in the molecular plane. To reduce the cost of the simulations, all these calculations were performed at the AM1 l 2 semiempirical level. The corresponding MSP maps (Figure 6) contain more than 7500 points for pyrimidines and 10 OOO points for purines, where each point corresponds to a MSP-optimized water molecule.
n
Alhambra et al.
3090 J. Phys. Chem., Vol. 99, No. 10, 1995 GUANINE
ADENiKE
I
28
t-
i4
z3 16 12
I
t
c -
I -
'r
-&I r -'
17
!-
"i
I
TEYM IN E
75
:5
!
11
-
-
7 -
c
I -l
I
-15
-11
-7
-3
!
1
~
l
5
I
1
8
,
13
,
,
r
,
17
Figure 6. MSP isopotential maps of the DNA bases (adenine, guanine, cytosine, thymine). Energy difference (in kcal/mol) between isopotential contour lines is 1 kcal/mol. Distance is given in atomic units.
van der Waals parameters for the bases were transferred from the values in Table 2. The maps are in general smooth and regular, and a detailed picture of the finest features is provided. This demonstrates the quality of the optimization algorithm used to find the optimum orientation of the water at each point. As expected from their free energies of hydrati~n,~ the bases interact strong1 with the water in the molecular plane. Thus, even at 5-6 from the molecule the interaction energies are close to -1 kcal/ mol (Figure 6). The minima are typically located in the lone pair region or in the vicinity of the acidic protons. It is interesting to note the stability of bifunctional waters, which interact simultaneously as donor and acceptor of hydrogen bonds with the bases. Examples of this kind of interaction are the waters acting as hydrogen bond donor for carbonyl oxygens and acceptor for nitrogens of guanine, cytosine, and thymine.
f
Another example is the water which acts as donor for N7 and as acceptor for the N6 amino group of adenine. In contrast, bifurcated hydrogen bonds, such as HzO--HzN, are less stable than the linear interactions H20--H-N-H (see maps for guanine or adenine). All the MSP minima for the different bases were located using the ab initio 6-31G* basis set. The MSP-6-31G* optimized geometries were used to compute AM1-MSP and ab initio 6-31G* SCF interaction energies (BSSE'O was corrected in all the SCF calculations). Results in Table 4 clearly show the suitability of the MSP-6-3 lG* to represent SCF interaction energies. Thus, the rms deviation between 6-31G* SCF and MSP values is only 0.37 kcal/mol. The ratio SCFMSP is 0.97, which reveals a systematic overestimation of the SCF values around 3%. Finally, the correlation coefficient of 0.99 between MSP and SCF values indicates an almost perfect correlation. It
J. Phys. Chem., Vol. 99, No. IO, 1995 3091
Molecular Solvation Potential TABLE 4: Ab Initio HF/SCF(6-31G*), MSP(6-31G*), and MSP(AM1) Interaction Energy (in kcal/rnol) for Minima Found in Isopotential Maps of Figure 6 CYt N4 N3-02 N1-02 thy 04 N3 N1-02 gua N9 N3 N1-06 N7-06 ade N9 N3 Nl-N6 N6-N7
SCF(6-31G*)
MSP(6-31G*)
MSP(AM1)
-4.01 -7.78 -9.96
-3.72 -8.16 -9.80
-2.75 -6.97 -6.75
-4.20 -6.15 -8.07
-3.93 -6.36 -8.45
-3.49 -4.36 -6.03
-5.08 -5.03 -8.34 -6.44
-5.56 -5.31 -9.03 -6.57
-4.93 -4.67 -6.52 -5.65
-7.08 -4.67 -5.18 -5.95
-7.22 -4.85 -5.77 -6.45
-6.06 -4.94 -5.02 -5.43
must be emphasized that the cost of the MSP calculation is only a fraction of that required for the SCF calculation. For instance, 1000 single point water-adenine interactions at the SCF level implies about 16 CPU days in an IBM 3090 600VF computer, while the same calculation at the MSP level takes only 900 s. The excellent agreement between MSP and SCF results for the four DNA bases confirms (i) the suitability of the MSP to represent hydration interactions in biologically-relevant molecules and (ii) the general goodness of the “first-generation” van der Waals parameters reported in Table 2. The agreement between MSP-AM1 and ab initio 6-31G* SCF values is satisfactory, but not as good as that found at the MSP6-31G* level, as shown by the rms deviation of 1.4 kcal/mol. This can be realized from three different sources of error: (i) the intrinsic shortcomings of the AM1 wave function, (ii) the problems of the semiempirical NDDO MEP in inner regions,12c (iii) the van der Waals parameters obtained from ab initio calculations and the core-water interactions, which can be quite different from those in semiempirical calculations. However, in spite of all these shortcomings, the MSP-AM1 performs reasonably well and offers a good alternative for very large molecules, whose wave functions cannot be computed at the ab initio 6-31G* level. In this paper we have presented the development and optimization of a new parameter to study solvation phenomena: the molecular solvation potential. The MSP provides at a very reduced computational cost an accurate picture of the structural and energetic characteristics of hydration for medium and large molecules. The method can be of general use in the calculation of other interactions, where a particle (or particles) needs to be represented at a high level of quality, while the rest of the system can be roughly represented at the classical level. The method could evidently be improved. For instance, a more precise van der Waals re-parametrization would increase the quality of the result. The introduction of a classical waterwater potential function would allow us to perform MC, MD, or even free energy perturbation calculations using the MSP as an effective Hamiltonian. In addition, the use of solvent-adapted wave functions obtained from SCRF calculation would allow us to include the statistically averaged polarization effects. Specific polarization interactions were not found to be extremely relevant in the systems studied. However, if their inclusion seems necessary, they might be easily introduced by means of induced dipoles, avoiding the costly recalculation of the SCF wave function for each solvent movement. The method is very
efficient computationally and has a large number of possibilities, which make it a very useful parameter for the study of molecular interactions in large systems. Acknowledgment. We thank the Centre de Supercomputaci6 de Catalunya (CESCA) and the Direcci6n General de Investigaci6n Cientifica y TBcnica (DGICYT PM92-0055-C02-01) for financial support. We thank Dr. J. Tirado-Rives for making available a copy of his SIMPLEXFIT program. Finally, we also thank R. Rycroft for his technical assistance in preparing the manuscript. References and Notes (1) (a) Scrocco, E.; Tomasi, J. Top. Curr. Chem. 1973, 42, 95. (b) Politzer, P., Truhlar, D. G., Eds. Chemical Applications of Atomic and Molecular Electrostatic Potentials; Plenum: New York, 1981. (c) Politzer, P.; Daiker, K. C. In The Force Concept in Chemistry; Deb, B. M., Ed.; Van Nostrand Reinhold New York, 1981; p 294. (d) Nhay-Szab6, G.; S u r j h , P. R. In Theoretical Chemistry of Biological Systems; Nhay-Szab6, G., Ed.; Elsevier: Amsterdam, 1986, p 1. (e) Tomasi, J.; Alagona, G.; Bonaccorsi, R.; Ghio, C.; Cammi, R. In Theoretical Models of Chemical Bonding, Maksic, Z., Ed.; Springer: Berlin, 1991; Vol. 4. (f) Politzer, P.; Murray, J. S. In Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; VCH Publishers: New York, 1991; Vol. 2, p 273. (2) Orozco, M.; Luque, F. J. J. Comput. Chem. 1993, 14, 587. (3) Orozco, M.; Luque, F. J. Biopolymers 1993, 33, 1851. (4) (a) McCammon, J. A.; Harvey, S. C. Dynamics of Proteins and Nucleic Acids; Cambridge University: Cambridge, 1987. (b) Bash, P. A,; Singh, U. C.; Langridge, R.; Kollman, P. A. Science 1987, 236, 564. (c) van Gunsteren, W. F.; Weiner, P. K., Eds. Computer Simulation of Biornolecular Systems; ESCOM Science Publishers: Leiden, 1989. (d) Karplus, M.; Petsko, G. A. Nature 1990, 347, 631. (e) van Gunsteren, W. F.; Berendsen, H. J. C. Angew. Chem. Int. Ed. Engl. 1990, 29, 992. (f) Jorgensen, W. L. Chemtracts: Org. Chem. 1991,4, 91. (5) (a) Weiner, S. J.; Seibel, G. L.; Kollman, P. A. Proc. Natl. Acad. Sci. U.S.A. 1986, 83, 649. (b) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (c) Warshel, A. Computer Modeling of Chemical Reactions in Enzymes and Solutions; Wiley: New York, 1991. (d) Luzhkov, V.; Warshel, A. J. Comput. Chem. 1992, 13, 199. (e) Gao, J. J. Phys. Chem. 1992,96,537. (f) Gao, J.; Pavelites, J. J. J. Am. Chem. SOC. 1992, 114, 1912. (g) Gao, J.; Xia, X. Science 1992, 258, 631. (6) (a) Williams, D. E. J. Comput. Chem. 1988, 9, 745. (b) Essex, J. W.; Reynolds, C. A.; Richards, W. G. J. Am. Chem. SOC. 1992,114,3634. (c) Carlson, H. A.; Nguyen, T. B.; Orozco, M.; Jorgensen, W. L. J. Comput. Chem. 1993, 14, 1240. (d) Orozco, M.; Jorgensen, W. L.; Luque, F. J. J. Comput. Chem. 1993, 14, 1498. (e) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollman, P. A. J. Am. Chem. Soc. 1993,115,9620. (f) Rauhut, G.; Clark, T. J. Comput. Chem. 1993, 14, 503. (9) Alemfin, C.; Orozco, M.; Luque, F. J. Chem. Phys., in press. (7) (a) Miertus, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981,55, 117. (b) Karelson, M. M.; Tamm, T.; Katritzsky, A. R.; Cato, S. J.; Zerner, M. Tetrahedron Comput. Method. 1989, 12, 295. (c) Negre, M.; Orozco, M.; Luque, F. J. Chem. Phys. Lett. 1992,196.27. (d) Chudinov, G. E.; Napolov, D. V.; Basilevsky, M. U. Chem. Phys. 1992, 160, 41. (e) Wang, B.; Ford, G. P. J. Comput. Chem. 1992, 97, 4162. (0 Cramer, C. J.; Truhlar, D. G. Science 1992,256,213. (g) Rauhut, G.; Clark, T.; Steinke, T. J. Am. Chem. SOC. 1993, 115, 9174. (h) Bachs, M.; Luque, F. J.; Orozco, M. J. Comput. Chem. 1994, 15, 446. (i) Rivail, J. L.; Rinaldi, D. Chem. Phys. 1976, 18, 223. (j) Rivail, J. L.; Cartier, A. Mol. Phys. 1978.36, 1085. (k) Floris, F.; Persico, M.; Tani, A,; Tomasi, J. Chem. Phys. Lett. 1992, 199, 518. (8) (a) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981; p 331. (b) Warshel, A.; Russel, S. Q. Rev. Biophys. 1984, 17, 283. (c) Berendsen, H.J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987,91,6269. (d) Jorgensen, W. L.; Chandrasekhar,J.; Madura, J.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (e) van Gunsteren, W. F.; Luque, F. J.; Timms, D.; Torda, A. E. Annu. Rev. Biophys. Biomol. Struct. 1994, 23, 847. (9) (a) Momany, F. A. J. Phys. Chem. 1978,82, 592. (b) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984,5, 129. (c) Orozco, M.; Luque, F. J. J. Comput. Chem. 1990, 11, 909. (d) Williams, D. E. In Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; VCH Publishers: New York, 1991; Vol 2, p 219. (10) Hariharan,P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213. (11) Boys, S . F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (12) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. P. J. Am. Chem. Soc. 1985, 105, 3902. (13) (a) Ferenczy, G. G.; Reynolds, C. A; Richards, W. G. J. Comput. Chem. 199O,II,431. (b) Alhambra, C.; Luque, F. J.; Orozco, M. J. Comput. Chem. 1994, 15, 12.
3092 J. Phys. Chem., Vol. 99, No. IO, 1995 (14)Luque, F.J.; Orozco, M. MOPETE Computer program, University of Barcelona, 1993. (15) Dupuis, M.; Rys, J.; King, H. F. QCPE. Bull. 1977,58, 17. (16)Stewart, J. J. P. QCPE Bull. 1984,4 , 109.Version modified by J. Bofill and S. Olivella, 1990,and F. J. Luque and M. Orozco, 1993. (17) Tirado-Rives, J. SIMPLEXFIT computer program, Yale University, 1993. (18) (a) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, P. A. J. Compur. Chem. 1977,8, 581. (b) Weiner, S. J.; Kollman, P. A.; Case, P. A.; Singh, U. C.; Ghio, C.; Alagona, G. J. Am. Chem. Soc. 1984,106,765. (19) (a) Jorgensen, W.L.; Swenson, C. J. J. Am. Chem. Soc. 1985,107,
Alhambra et al. 1489.(b) Jorgensen, W. L. Bioorganic Simulation System (BOSS), Version 3.4,Yale University, 1993. (20)The number of different orientations is indicated in parentheses. (21) Svishchev, I. M.; Kusalik, P. G. J. Chem. Phys. 1993,99, 3049. (22) Suzuki, S.;Green, P. G.; Burmgarner, R. E.; Dasgupta, S.; Goddard III, W. A.; Blake, G. A. Science 1992,257,943. (23)Thomas, K. A.; Smith,G. M.; Thomas, T. B.; Feldmann, R. J. Proc. Natl. Acad. Sci. U.S.A. 1982, 79, 4843. (24)Jorgensen,W. L.;Severance,D.L. J. Am.Chem.Soc. 1990,112,4768, JP941514L
