- Email: [email protected]
SCRF models: influence of the number of discrete solvent molecules
Chemical Physics 244 Ž1999. 151–162 www.elsevier.nlrlocaterchemphys
Hydration of cytosine using combined discreterSCRF models: influence of the number of discrete solvent molecules Carlos Aleman ´
)
Departament d’Enginyeria Quımica, E.T.S. d’Enginyers Industrials, UniÕersitat Politecnica de Catalunya, Diagonal 647, Barcelona ´ ` E-08028, Spain Received 11 January 1999
Abstract The solvation of cytosine has been investigated using a combined discreterSCRF model. Molecular structures and gas phase binding energies were obtained at the HFr6-31GŽd. level for ten cytosine–nH 2 O complexes, where n ranges from 1 to 13. Furthermore, gas-phase binding were also computed at the MP2r6-31GŽd. level for complexes with n F 7. The Polarizable Continuum Model self-consistent reaction-field ŽSCRF. were applied to such complexes in order to investigate the effect of the explicit water on the free energy of solvation and on the electronic structure of the solute. The results were compared with those provided by SCRF methods revealing that the discreterSCRF models constitute a very reliable strategy. q 1999 Published by Elsevier Science B.V. All rights reserved.
1. Introduction Most of the current investigations in theoretical chemistry are based on the study of molecules immersed in a solvent phase. The theoretical treatment of the solvent can be broadly classified in two categories: Ži. discrete models; and Žii. continuum models. In the discrete models, solvent molecules are treated explicitly allowing the description of both the microscopic structure of the solvent and the specific solute–solvent interactions. In practice, if the system is described using quantum mechanics, the applicability of the discrete models is restricted to a selected number of configurations of a solute surrounded by a few solvent molecules w1,2x. Conversely, force-field simulations using empirical interaction potentials al)
Tel.: q34-93401-6680r6681; fax: q34-93401-6600r7150; e-mail: [email protected]
low calculations with a large number of both solvent molecules and configurations w3x. However, they cannot account for the modifications induced by the solvent in the electronic structure of the solute. In this connection, hybrid methods including a classical simulation for the solvent terms and solute–solvent interactions coupled to a quantum mechanical description of the solute provide excellent results w4,5x. On the other hand, the continuum models allow the incorporation of solvent effects but neglect the microscopic solvation structure in the vicinity of the solute. Self-consistent reaction-field ŽSCRF. methods based on continuum dielectrics have received much attention within the quantum mechanical framework. A large number of SCRF methods have been developed w6–23x, all of them providing special attention to the determination of the electrostatic contribution to the free energy of solvation Ž DGsol ., which is evaluated assuming that the infinite dielectric
0301-0104r99r$ - see front matter q 1999 Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 Ž 9 9 . 0 0 1 2 2 - 6
152
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
medium, the solvent, reacts against the solute charge distribution generating a reaction field, which in turn interacts with the solute. The main differences between the SCRF methods lie w20x: Ži. in the shape of the solutersolvent interface; Žii. in the definition of the reaction field; and Žiii. in the evaluation of the non-electrostatic contributions to DGsol . For the last few years, SCRF methods have been applied to a wide number of chemical systems w18x. However, electronic effects associated with specific solute–solvent interactions are not considered by these methods. In the middle of the 1970s, Claverie et al. w24x proposed the use of combined discreterSCRF models to overcome this deficiency. M ore recently, we have used com bined discreterSCRF methods to predict both the DGsol and electronic properties in aqueous solution of some natural products w25x and nucleic acid bases w26x. Comparison with the results provided by SCRF methods proved the reliability of discreterSCRF models, which for the study of electronic properties provide more information than either discrete or SCRF methods. However, it should be emphasized that such studies were performed considering a reduced number of explicit solvent molecules, i.e., the maximum number of explicit water molecules was 4 w25x or 5 w26x, depending on the hydrophilicity of the solute. Therefore, it would be interesting to investigate the changes on the electronic properties of the solute as well as in the predicted DGsol on the number of explicit solvent molecules used in the combined discreterSCRF models. In this paper, a systematic study about the influence of the number of explicit solvent molecules in the reliability of combined discreterSCRF models is presented. For this purpose, we investigated the hydration of cytosine nucleic acid base Žsee below., considering complexes with a number of water molecules ranging from 1 to 13. The work is outlined in the following sections. The next section summarises the quantum mechanical methods used in the calculations. Next, we present the results, which are divided in three parts. First, we show the gas-phase binding energies for cytosine–nH 2 O complexes, where n denotes the number of explicit water molecules, computed by the supermolecule approach. Second, we present the solvation of the complexes. In this section, the influence of the num-
ber of explicit water molecules on the DGsol of cytosine is investigated. After this, we analyze the electronic structure of the solute under the influence of the solvent. For this purpose, the atomic charges obtained for different complexes are compared. Finally, the conclusions of the work are presented.
2. Methods 2.1. Gas-phase calculations Molecular structures of four different cytosine– 1H 2 O complexes were optimized at the HFr631GŽd. level w27x. A frequency analysis was performed to verify the nature of the minimum state of the stationary points reached after geometry optimizations. In order to generate the cytosine–nH 2 O complexes with n s 3, 5, 7 and 9, the water molecules were distributed around the solute. First, we distributed the water molecules around the hydrophillic region of cytosine. Once this was surrounded by the first solvation shell, we put explicit solvent molecules in the hydrophobic region of the solute. The complexes were submitted to a complete geometry optimization at the HFr6-31GŽd. level. Single point calculations at the MP2r6-31GŽd. level w28x were performed for all the cytosine–nH 2 O complexes with n s 1, 3, 5 and 7. Cytosine–nH 2 O complexes with n s 11 and 13 were generated using as starting points the optimized complexes with n s 9 and 11, respectively, and adding two new water molecules in the closest positions to the solute. Geometry optimizations of the added water molecules were performed at the HFr6-31GŽd. level but fixing the cytosine and the remaining water molecules.
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
The basis set superposition error ŽBSSE. was corrected at both HF and MP2 levels using the standard counterpoise ŽCP. method w29x. The gasphase binding energy Ž D Egp . was computed as the difference of the total energy of the complex Ž Ecyt – nH2O . and the sum of the solute and nH 2 O energies Ž Ecyt , EnH2O .. D Egp s EcytynH 2 O y Ecyt y EnH 2 O
Ž 1.
153
The reaction field operator was described in terms of an apparent set of point charges ŽEq. Ž4.. spread over the solute–solvent interface, i.e., the solute cavity. The apparent charge density was determined by solving the Laplace equation at the solute–solvent interface. The electrostatic contribution to the DGsol was determined following linear free energy response theory ŽLFER. according to Eq. Ž5. w32x. VR s S i s Ž si . Sir< r y r 0 < s S i qir< r y r 0 <
Ž 4.
DGele s-C < H
y 1r2 Ž - C < V
In a recent work, we compared the reliability of some of the most common SCRF methods in the framework of combined discreterSCRF models w23x. We found that the polarizable continuum model ŽPCM. adapted to the AM1 semi-empirical method ŽPCMrAM1. w20,30x provided excellent results at a very reasonable computational cost. Moreover, the differences between the results provided by the semi-empirical PCMrAM1 and the PCM at the ab initio HFr6-31GŽd. level ŽPCMr6-31GŽd.. were smaller than 1 kcalrmol, which is the error associated to these methods w25x. Therefore, we selected the PCMrAM1 method to describe the bulk water. According to this method, the DGsol was determined as the addition of electrostatic and steric contributions ŽEq. Ž2... The steric contribution was computed as the sum of the cavitation and a van der Waals.
Ž 5.
2.3. Discreter SCRF model The free energy of solvation for the solute with n EW . explicit water molecules Ž DGsol was computed as the addition of the gas-phase binding energy between the solute and the explicit water molecules and the DGsol for the complex ŽEq. Ž6... EW DGsol s D Egp q DGsol Ž cytosine–nH 2 O .
Ž 6.
The cavitation term was determined using Pierotti’s scaled particle theory w31x, while the van der Waals term was evaluated by means of a linear relation with the molecular surface area w20,30x. The PCMrAM1 considers the solvent as a continuum dielectric with the cavity accurately modeled on the solute. The solvent reacts against the solute charge distribution, generating a reaction field Ž VR .. The electrostatic interaction between the solute and the solvent is introduced as a perturbation operator in the solute Hamiltonian ŽEq. Ž3...
It should be noted that in Eq. Ž6., the entropic contributions corresponding to the association of the cytosine and the cluster of n water molecules are neglected in both the gas phase and solution terms. Owing to the chemical chemical nature of the systems under study, i.e., a rigid solute surrounded by the first hydration shell, the entropic contributions are mainly due to the six external degrees of freedom which become internal degrees in the complex. Nevertheless, the entropic contribution in the gas phase and in solution can reasonably be assumed to be very similar, and therefore a considerable compensation is expected. EW Ideally, DGsol should be equivalent to the free energy of solvation of the molecular system comAC . puted in an all-continuum approach Ž DGsol , which is determined as the sum of the free energies of solvation for the solute and the water molecules involved in the molecular system:
Ž H 0 q VR . C s EC
AC s DGsol Ž cytosine. q n )DGsol Ž H 2 O . DGsol
DGsol s DGele q DGcav q DGvdW
Ž 2.
Ž 3.
Ž 7.
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
154
The free energy of solvation for cytosine provided Dr SCRF . by the combined discreterSCRF model Ž DGsol was determined according to Eq. Ž8.: D r SCRF DGsol Ž cytosine. EW s DGsol Ž cytosine–nH 2 O . y n )DGsol Ž H 2 O .
s D Egp q DGsol Ž cytosine–nH 2 O . y n )DGsol Ž H 2 O .
Ž 8. Dr SCRF DGsol
It is worth noting that has an incorrect limit for n ™ ` since D Egp remains constant when a large number of explicit solvent molecules is considered, whereas nUDGsolŽH 2 O. increases faster than DGsolŽcytosine–nH 2 O.. Gas-phase calculations were performed with the Gaussian-94 w33x computer package. PCMrAM1 calculations were performed with an adapted version of the MOPAC 93 w34x computer program 1.
3. Results and discussion 3.1. Gas-phase calculations. The optimized geometries of the four cytosine– 1H 2 O complexes are shown in Fig. 1, where the intermolecular geometric parameters are also displayed. All the complexes are true ground-state minima, which is verified from second derivatives calculations. The molecular geometry of one of these complexes, i.e., that in which the water molecule interacts with the H1 and O2 atoms of cytosine ŽFig. 1Žd.., was recently optimized at the MP2r6-31GŽd. level w35x. The small differences between the intermolecular geometric parameters obtained at the HFr6-31GŽd. and MP2r6-31GŽd. levels indicate that the former is a suitable method of performing geometry optimizations of hydrogen-bonded complexes. More specifically, the donor-acceptor distances predicted at the HF level are overestimated by ˚ with respect to the MP2 values, about 0.05–0.10 A whereas the donor–hydrogen-acceptor angles are underestimated by about 1–28.
1
Adapted to perform PCM calculations by F.J. Luque and M. Orozco
Fig. 2 shows the cytosine–nH 2 O complexes with n s 3, 5, 7, 9, 11 and 13. It is worth noting that water molecules form a network of hydrogen bonds when a large enough number of water molecules is considered. These stable interactions allow the surrounding of the hydrophobic region of cytosine with water molecules in the complexes with n s 11 and 13. Two different types of water molecules can be detected in the cytosine–nH 2 O complexes: Ži. those that form at least one hydrogen bond with the solute; and Žii. those that form hydrogen bonds with other water molecules but not with the solute. Table 1 lists the number of water molecules of each type for the different complexes investigated together with the gas-phase binding energies obtained from Eq. Ž1.. It is worth noting that the maximum number of solvent molecules that are hydrogen bonded to the cytosine, i.e., those that play a crucial role in specific solute– solvent interactions, is six. An inspection of the energies displayed in Table 1 and the molecular geometries depicted in Figs. 1 and 2 indicates that the contributions to the gas-phase binding energies of the water molecules involving one and two hydrogen bonds with cytosine are about 5 and 7–10 kcalrmol, respectively. In contrast, water molecules surrounding the hydrophobic region of cytosine present a weak interaction with the solute, their contribution to the gas-phase binding energy being only 0.5–1.5 kcalrmol. Comparison between the gas-phase binding energies estimated at the HF and MP2 levels ŽTable 1. indicates that the two methods provide very similar values from a qualitative point of view. Accordingly, the binding energy order predicted by the two methods for the monohydrated complexes is the same, the HF binding energies being systematically less attractive than those for MP2. Such underestimation corresponds, on average, to an error of about 8%. Similar results were obtained in previous works for related hydrogen bonded systems w25,26,36x. On the other hand, it is worth noting that the relative energy order predicted for the four monohydrated complexes can easily be rationalized considering the hydrogen bonding parameters Žsee Fig. 1.. A detailed analysis of both cytosine–water and water–water hydrogen bonding interactions for cytosine–nH 2 O complexes is presented in Table 2, where the hydrogen-acceptor distances and the donor–hy-
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
155
Fig. 1. Molecular geometries of the four cytosine–1H 2 O complexes optimized at the HFr6-31GŽd. level. The hydrogen bonding parameters are also displayed.
drogen-acceptor angles are listed for the complexes with n s 3, 5, 7 and 9. The hydrogen bonding ˚ larger for the cytosine– distances are about 0.1 A water interactions than for the water–water interactions. On the other hand, an inspection of the complexes with n s 7 and 9 reveals that the donor–hydrogen-acceptor angles are less favorable for the water–water interactions than for the cytosine–water interactions. Similar conclusions are reached from the analysis of the complexes with n s 11 and 13. 3.2. Aqueous solution calculations The DGsol values predicted by the PCMrAM1 method for all the complexes investigated are listed in Table 3, where both the electrostatic and steric
contributions are also displayed. It is worth nothing that the electrostatic contribution is the most important contribution to DGsol . This attractive contribution mainly depends on the accessibility to the bulk solvent of the polar atoms contained in both the cytosine and explicit water molecules. Accordingly, the electrostatic contribution increases with the surface accessible to the bulk solvent. This feature is clearly evidenced in the monohydrated complexes. Thus, in the complex displayed in Fig. 1Žc., the two hydrogen atoms contained in the water molecule are accessible to the bulk solvent, whereas in the other three monohydrated complexes only one atom of the water molecule points to the bulk solvent. Therefore, the attractive electrostatic interactions are larger for the complex displayed in Fig. 1Žc. than for the
156
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
Fig. 2. Molecular geometries of the cytosine–nH 2 O complexes optimized at the HFr6-31GŽd. level: Ža. n s 3; Žb. n s 5; Žc. n s 7; Žd. n s 9; Že. n s 11 and Žf. n s 13.
others. On the other hand, the steric term provides a repulsive contribution since the polarizable nature of
water makes this solvent to be highly structured, and therefore, generation of cavity is difficult from an
C. Alemanr ´ Chemical Physics 244 (1999) 151–162 Table 1 Number of water molecules that are hydrogen bonded to cytosinea Ž NHB S . and to other water moleculesb Ž NHBW . in cytosine– nH 2 O complexes. Gas-phase binding energies c ŽEq. Ž1.. computed at the HFr6-31GŽd. and MP2r6-31GŽd. levels of theory Complex d
NHB S NHBW D EHFr6-31GŽd. D E MP2r6-31GŽd.
Cytosine–1H 2 O Ža. Cytosine–1H 2 O Žb. Cytosine–1H 2 O Žc. Cytosine–1H 2 O Žd. Cytosine–3H 2 O Cytosine–5H 2 O Cytosine–7H 2 O Cytosine–9H 2 O Cytosine–11H 2 O Cytosine–13H 2 O
1 1 1 1 3 4 5 6 6 6
– – – – – 1 2 3 5 7
y7.1 y9.3 y5.4 y9.8 y20.7 y31.0 y36.1 y46.0 y47.0 y50.1
y7.3 y10.3 y6.2 y10.4 y22.8 y34.1 y39.6 – – –
157
Table 3 Free energies of solvationa predicted for the cytosine– nH 2 O complexes by the PCMrAM1 method. Electrostatic and steric Ž DGcav q DGvdW . contributionsa are included. Complex
DGsol
DGele
DGcav q DGvdW
Cytosine–1H 2 O Ža. Cytosine–1H 2 O Žb. Cytosine–1H 2 O Žc. Cytosine–1H 2 O Žd. Cytosine–3H 2 O Cytosine–5H 2 O Cytosine–7H 2 O Cytosine–9H 2 O Cytosine–11H 2 O Cytosine–13H 2 O
y16.3 y16.1 y21.0 y17.3 y19.8 y18.3 y21.4 y23.8 y30.5 y35.8
y20.5 y20.3 y25.2 y21.6 y25.3 y25.2 y29.3 y32.8 y40.8 y47.3
4.2 4.2 4.2 4.3 5.5 6.9 7.9 9.0 10.3 11.5
a
In kcalrmol.
a
Number of water molecules with one or more hydrogen bonding interactions with cytosine. b Number of water molecules that are hydrogen bonded to other water molecules but not cytosine. c In kcalrmol. d The labels of the cytosine–1H 2 O complexes correspond to those displayed in Fig. 1.
energy point of view. The repulsive cavitation contribution is not compensated by the van der Waals interactions, which, due to the small size of the cavity, are not attractive enough.
Table 2 Cytosine . . . water a and water . . . water hydrogen bonding parametersb measured in the optimized molecular geometries of the cytosine–nH 2 O complexes with n s 3, 5, 7 and 9 Complex
Cytosine–3H 2 O
Cytosine–5H 2 O
Cytosine–7H 2 O
Cytosine–9H 2 O
a
Cytosine . . . water
Water . . . water
Distances
Angles
Distances
Angles
O8 . . . Hw s 2.13 H12 . . . 6Ow s 2.04 H13 . . . Ow s 2.11 O8 . . . Hw s 2.00 H11 . . . Ow s 2.08 N3 . . . Hw s 1.99 H12 . . . Ow s 2.17 H13 . . . Ow s 2.03 O8 . . . Hw s 1.90 H11 . . . Ow s 1.96 N3 . . . Hw s 2.01 H12 . . . Ow s 2.10 H13 . . . Ow s 2.11 O8 . . . Hw s 1.91 O8 . . . Hw s 1.94 H11 . . . Ow s 2.11 N3 . . . Hw s 2.07 H12 . . . Ow s 2.13 H13 . . . Ow s 2.07
O8 . . . Hw–Ow s 166.6 N7–H12 . . . Ow s 161.8 N7–H12 . . . Ow s 172.5 O8–Hw–Ow s 147.9 N1–H11 . . . Ow s 146.5 N3 . . . Hw–Ow s 165.8 N7–H12 . . . Ow s 167.6 N7–H13 . . . Ow s 153.3 O8 . . . Hw–Ow s 166.3 N1–H11 . . . Ow s 174.6 N3 . . . Hw–Ow s 165.8 N7–H12 . . . Ow s 173.1 N7–H13 . . . Ow s 146.0 O8 . . . Hw–Ow s 167.8 O8 . . . Hw–Ow s 178.8 N1–H11 . . . Ow s 177.9 N3 . . . Hw–Ow s 169.6 N7–H12 . . . Ow s 167.2 N7–H13 . . . Ow s 148.8
2.00 – – 1.84 1.94 1.93 – – 1.88 1.85 2.00 1.86 2.12 1.85 1.97 1.83 1.99 1.85 2.13
164.2 – – 159.3 173.5 167.9 – – 159.1 156.9 165.0 159.9 147.5 157.6 160.7 154.7 165.8 160.3 147.6
The numeration for the atoms of cytosine is displayed in the Scheme I while H W and O W refer to the hydrogen and oxygen atoms of the water molecules. b ˚ and degrees, respectively. Hydrogen-acceptor distances and donor–hydrogen-acceptor angles are in A
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
158
AC EW Table 4 shows the values of DGsol and DGsol computed according to Eqs. Ž6. and Ž7., respectively, as well as the difference between them, which can be considered as a measure of the reliability of the combined discreterSCRF model. As expected from EW estimated the results displayed in Table 1, the DGsol from the HFr6-31GŽd. and MP2r6-31GŽd. gasphase binding energies are very similar. FurtherAC more, the agreement between DGsol and both EW EW Ž DGsol HF. and DGsol ŽMP2. is encouraging, the maximum difference between these parameters being 12.2 and 5.5 kcalrmol, respectively. Such differences correspond to an error of 12.4% and 14.8% for the cytosine–13H 2 O and cytosine–3H 2 O complexes, respectively. Fig. 3 represents the variation of the error Žin %. corresponding to the difference AC EW Ž between DGsol and both DGsol HF . and EW Ž DGsol MP2. of the number of explicit solvent molecules. Inspection to the data derived from HFr6-31GŽd. gas-phase binding energies indicates that good results are achieved for cytosine–nH 2 O complexes with n s 5, 7 and 9, i.e., the error is lower than 6.5%, whereas the error increases when the number of explicit water molecules is either too
Table 4 Free energies of solvationa,b predicted using an all continuum AC . EW . approachc Ž DGsol and explicit solvent molecules Ž DGsol for AC cytosine– nH 2 O complexes. The differences b w DGsol y AC EW Ž EW Ž DGsol HF.x and w DGsol y DGsol MP2.x are listed in parenthesis in the 2nd and 3rd columns, respectively Complex
AC DGsol
EW Ž DGsol HF. d
EW Ž DGsol MP2. d
Cytosine–1H 2 O Ža. Cytosine–1H 2 O Žb. Cytosine–1H 2 O Žc. Cytosine–1H 2 O Žd. Cytosine–3H 2 O Cytosine–5H 2 O Cytosine–7H 2 O Cytosine–9H 2 O Cytosine–11H 2 O Cytosine–13H 2 O
y24.9 y24.9 y24.9 y24.9 y37.1 y49.3 y61.5 y73.7 y85.9 y98.1
y23.4 ( y 1.5) y25.4 (0.5) y26.4 (1.5) y27.2 (2.2) y40.5 (3.4) y49.3 (0.0) y57.5 ( y 4.0) y69.8 ( y 3.9) y77.5 ( y 8.4) y85.9 ( y 12.2)
y23.6 ( y 1.3) y26.4 (1.5) y27.2 (2.3) y27.7 (2.8) y42.6 (5.5) y52.4 (3.1) y61.0 ( y 0.5) – – –
a
AC EW DGsol and DGsol computed according to Eqs. Ž6. and Ž7., respectively. b In kcalrmol. c The free energy of solvation of the water molecule used to AC compute DGsol is y6.1 kcalrmol. d EW Ž EW Ž DGsol HF. and DGsol MP2. refer to the values computed determined from HFr6-31GŽd. and MP2r6-31GŽd. gas-phase binding energies, respectively.
Fig. 3. Variation of the error corresponding to the difference AC EW Ž EW Ž between DGsol and both DGsol HF. ŽI. and DGsol MP2. ŽB. upon the number of explicit water molecules.
large or too small. The errors corresponding to the AC EW Ž differences between DGsol and DGsol MP2. suggest similar trends. These results indicate that the strategy of including explicit water molecules is very reliable when the number of solvent molecules included in the calculation is sufficient to surround the hydrophillic region of the solute. In contrast, when the number of water molecules is too large, such as in cytosine– 11H 2 O and cytosine–13H 2 O complexes, the error can be greater than 12%. Table 1 shows that these complexes contain an important number of explicit water molecules without cytosine–water hydrogen bonding interactions Ž NHB W s 5 and 7, respectively.. Inspection of Fig. 2Že. and 2Žf. reveals that these water molecules are mainly located in the hydrophobic region of the solute or behind the water molecules hydrogen-bonded to the cytosine. On the other hand, when the number of explicit solvent molecules is too small, such as in cytosine–1H 2 O and cytosine–3H 2 O complexes, the error can be greater than 9%. In these cases, the number of water molecules is not sufficient to surround the hydrophillic region of cytosine completely ŽFigs. 1 and 2Ža... The possible errors in the combined discreter SCRF model can mainly be attributed to four different sources. The first source is the complexity of the molecular surface in complexes. The building-up of the surface, as determined in PCMrAM1 w37x, defines a large number of small spheres not assigned specifically to any atom, but needed to describe
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
properly a smooth solute–solvent interface and, consequently, to estimate correctly the solvent reaction field. The number of such spheres increases with the complexity of the surface, providing a less accurate prediction of DGsol . The second source of error is related to the use of solvent-accessible surfaces which are not the most appropriate to describe aggregation processes in solution. Thus, previous studies of intermolecular interactions in solution indicated that the solvent-accessible surface is not able to describe properly the area of the cavity not accessible to the bulk solvent w38x. Furthermore, the error associated with this surface model increases with the contact distance until the rupture of the complex. Inspection of Fig. 2 and Table 4 indicates that the errors associated to these deficiencies are minimized when the number of water molecules is sufficient to fill up the hydrophillic region of cytosine, since the number of void spaces in all the possible cytosine–water and water–water dimers is very small. However, the complexity of the surface increases when the number of water molecules is either too small or too large. Furthermore, if the number of water molecules is too large, then the contact distance in some of the cytosine–water and water–water dimers involved increases and the worst description of the volume inaccessible to the bulk solvent is obtained. The third
159
source of error is the neglection of the vibrational and rotational nuclear movements. As was stated in Section 2, the vibrational contribution is solution is compensated by that of the gas phase. However, the rotational nuclear motions provide an unfavorable contribution to DGsol due to the frictional solvent effects on the rotation of the solute system. It should be emphasized that the rotational nuclear motions do not have a direct effect on the solute charge distribution and have been included in the parametrization of the nonelectrostatic terms of the PCMrAM1 method, However, the parameters were developed for single solutes and not for complexes. Finally, it should be noted that we are considering frozen geometries. This approximation is not very important for rigid solutes like cytosine but could lead to dramatic errors in the study of flexible compounds. Dr SCRF The values of DGsol predicted for cytosine ŽEq. Ž8.. using the complexes with 5, 7 and 9 explicit water molecules are displayed in Table 5. The DGsol values estimated for cytosine w39–42x and its methylated derivative w39,43–45x using other computational approaches are also included in the Dr SCRF values derived table for comparison. The DGsol from HFr6-31GŽd. gas-phase binding energies are lower than those obtained from MP2r6-31GŽd. energies. They range from y14.8 to y18.8 kcalrmol
Table 5 Free energy of solvationa Ž DGsol . predicted for cytosine using the combined discreterSCRF methodb . DGsol values predicted by other computational approaches for cytosine and its methylated derivative are also included for comparison DGsol
Method
Computational approach
DiscreterSCRF DiscreterSCRF DiscreterSCRF DiscreterSCRF DiscreterSCRF SCRF SCRF LD LD MC QMrMM SCRF SCRF SCRF OPLSrFEP OPLSrAMBER
Cytosine–5H 2 O complex using HFr6-31GŽd. gas-phase binding energy and PCMrAM1. y18.8 Cytosine–5H 2 O complex using MP2r6-31GŽd. gas-phase binding energy and PCMrAM1. y21.9 Cytosine–7H 2 O complex using HFr6-31GŽd. gas-phase binding energy and PCMrAM1. y14.8 Cytosine–7H 2 O complex and MP2r6-31GŽd. gas-phase binding energy and PCMrAM1 y18.3 Cytosine–9H 2 O complex and HFr6-31GŽd. gas-phase binding energy and PCMrAM1 y14.9 Cytosine using PCMrAM1 y18.8 Cytosine using SM2rAM1 y20.6 Cytosine using the noniterative LD model y14.2 Cytosine using the iterative LD model y18.0 Cytosine using mixed quantum mechanics and classical Monte Carlo simulations y16.3 Methylated derivative of cytosine using PCMrAM1 y16.1 Methylated derivative of cytosine using PCMr6-31GŽd. y15.1 Methylated derivative of cytosine using SM2rAM1 y18.7 Methylated derivative of cytosine using free energy preturbations with the OPLS force-field y20.1 Methylated derivative of cytosine using free energy preturbations with the AMBER force-field y12.7
a
In kcalrmol. DrSCRF DGsol values predicted using Eq. Ž8..
b
Ref. Present work Present work Present work Present work Present work Present work w39x w41x w42x w40x w43x w43x w39x w44x w45x
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
160
and from y18.3 to y21.9 kcalrmol, respectively. Dr SCRF On the other hand, the DGsol values predicted for cytosine using 11 and 13 explicit solvent molecules are y10.4 and y6.6 kcalrmol, respectively. These low values are mainly due to the underestimation of the gas-phase binding energies estimated at the HF level. Dr SCRF The DGsol values displayed in Table 5, i.e., those obtained using the 5, 7 and 9 water molecules, are within the range of previous estimations predicted by other computational approaches. Thus, the DGsol values predicted by SM2, PCM and iterative LD models Žy20.6, y18.8 and y18.0 kcalrmol, respectively. are very close to those obtained using the MP2r6-31GŽd. gas-phase energies. On the other hand, the value obtained from Monte Carlo quantum mechanical and molecular mechanical ŽMC QMrMM. simulations Žy16.3 kcalrmol. is close to Dr SCRF the DGsol values estimated from HFr6-31GŽd. gas-phase energies. These results suggest that the values predicted by both combined discreterSCRF model using HFr6-31GŽd. gas-phase binding energies and molecular mechanics methods are slightly underestimated. Unfortunately, there is not a quantitative experimental DGsol value for cytosine. It should be emphasized that the solvation studies performed for its methylated derivative w39,43–45x also provide a wide range of DGsol values. Table 6 shows the Mulliken charges of isolated cytosine in a gas phase provided by the AM1 hamil-
tonian, isolated cytosine in an aqueous solution obtained using the SCRF method and cytosine in an aqueous solution predicted by the combined discreterSCRF method for the complexes with 3, 7 and 13 explicit water molecules. Comparison between the charges obtained in the gas phase and aqueous solution points out the complexity of the polarizing effect of water on the solute charge distribution. The more important features are the enhancement of the negative charges of the O2 and N3 atoms, as well as the increase of the positive charge of the hydrogen atoms bonded to N1 and N4. Furthermore, there is also an appreciable increase in the charge of the hydrogen atoms located in the hydrophobic region of the cytosine. Comparison of the atomic charges derived for cytosine using the combined discreterSCRF method reveals very similar values for all the complexes investigated. This trend suggests that the number of explicit water molecules is not crucial to obtain a suitable description of the electronic effects induced by specific solute–solvent interactions. For instance, the largest difference between the charges obtained for the cytosine–3H 2 O and cytosine–13H 2 O complexes is about 0.02 e Žwhere e is the charge of an electron. for the O2 atom. Of particular interest is the comparison of the atomic charges provided by the combined discreter SCRF model and those obtained using the SCRF method without including any explicit solvent
Table 6 Mulliken atomic chargesa for cytosine computed in the gas phase and in aqueous solution. Charges in aqueous solution were obtained using the PCMrAM1 method and the combined discreterSCRF method considering 3, 7 and 13 explicit water molecules. Charges obtained for the cytosine–7H 2 O complex but without the SCRF are also listed Žin italics. Atom
Cytosine Žgas phase.
Cytosine ŽPCMrAM1.
Cytosine–3H 2 O
Cytosine–7H 2 O
Cytosine–13H 2 O
N1 C2 N3 C4 C5 C6 H– ŽN1. O5ŽC2. N4 H– ŽC5. H– ŽC6. H–N3 H–N5
y0.324 0.349 y0.299 0.237 y0.348 0.065 0.252 y0.344 y0.372 0.153 0.152 0.246 0.232
y0.308 0.364 y0.423 0.255 y0.348 0.086 0.301 y0.503 y0.349 0.179 0.186 0.278 0.281
y0.324 0.395 y0.443 0.266 y0.344 0.087 0.287 y0.479 y0.341 0.176 0.174 0.263 0.277
y0.327 (y0.334) 0.400 (0.396) y0.451 (y0.410) 0.273 (0.267) y0.345 (y0.351) 0.090 (0.081) 0.275 (0.273) y0.461 (y0.396) y0.342 (y0.362) 0.173 (0.159) 0.173 (0.154) 0.262 (0.260) 0.273 (0.261)
y0.320 0.401 y0.455 0.272 y0.347 0.092 0.275 y0.458 y0.345 0.175 0.184 0.259 0.274
a
In units of electron.
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
molecule. It is worth noting that charges provided the SCRF method are somewhat enhanced for some atoms with respect to those obtained for the cytosine–nH 2 O complexes. For instance, atomic charges predicted for C2, N3, H1 and O2 differ by 0.036, 0.028, 0.026 and 0.042 e, respectively, with respect to those obtained for the cytosine–7H 2 O complex. On the other hand, Table 6 also shows the atomic charges for cytosine–7H 2 O complex but without applying the SCRF, i.e., considering a discrete method. It is worth noting that the atomic charges obtained from the discrete and discreterSCRF methods are quite similar for all the atoms except for N3 and O2, which differ by 0.041 and 0.065 e, respectively. However, an inspection to the results listed in Table 6 indicates that polarizing effect in cytosine seems to be overestimated by SCRF methods.
4. Conclusions In this study, we investigated the influence of the number of explicit water molecules in the solvation of cytosine using a combined discreterSCRF model. For this purpose, ten cytosine–nH 2 O complexes, where n ranges from 1 to 13, were optimized at the HFr6-31GŽd. level. The gas-phase binding energies were computed at both HFr6-31GŽd. and MP2r631GŽd. levels of theory. A combined discreterSCRF model using the PCMrAM1 method was applied to compute the DGsol of the ten complexes. From the AC EW comparison of DGsol and DGsol , we found that the best results, i.e., the smallest error between the two magnitudes, are obtained when the explicit water molecules are distributed around the hydrophillic region of cytosine filling up all the solvation sites. On the other hand, when the number of water molecules is too small or too large both the complexity of the molecular surface and the use of solventaccessible surfaces increase the differences between AC EW DGsol and DGsol . The analysis of the change in the atomic charges of cytosine provided by the combined discreterSCRF model emphasizes the important role of the specific interactions between the solute and the solvent. Furthermore, the results indicate that SCRF methods tend to overestimate the polarizing effect of the solvent on the solute suggesting that more investigation is needed in this area.
161
Acknowledgements Author is indebted to Dr. M. Orozco and F.J. Luque for making available to him their version of MOPAC93 adapted to perform PCM calculations. Author thanks to CESCA for computational facilities. References w1x O. Tapia, F. Colonna, J.G. Angyan, J. Chem. Phys. 87 Ž1990. 875. w2x G. Alagona, A. Pullman, E. Scrocco, J. Tomasi, Int. J. Pept. Protein Res. 5 Ž1973. 251. w3x M.P. Allen, D.J. Tildesley, in: Computer Simulation of Liquids, Oxford University Press, Oxford, 1987. w4x M.J. Field, P.A. Basch, M. Karplus, J. Comput. Chem. 11 Ž1990. 700. w5x J. Gao, J. Am. Chem. Soc. 116 Ž1994. 1563. w6x L. Rivail, D. Rinaldi, Theor. Chim. Acta 32 Ž1973. 57. w7x M.W. Wong, M.J. Frisch, K.B. Wiberg, J. Am. Chem. Soc. 113 Ž1991. 4776. w8x S. Miertus, E. Scrocco, Tomasi, J. Chem. Phys. 55 Ž1981. 117. w9x M.M. Karelson, T. Tam, A.R. Katritzky, S.J. Cato, M. Zerner, Tetrahedron Comput. Methodol. 2 Ž1989. 295. w10x B. Wang, G.P. Ford, J. Comput. Chem. 7 Ž1992. 4162. w11x G.E. Chudinov, D.V. Napolov, M.U. Basilevsky, Chem. Phys. 160 Ž1992. 41. w12x C.J. Cramer, D.G. Truhlar, Science 256 Ž1992. 213. w13x V. Luzhkov, A.J. Warshel, Comput. Chem. 13 Ž1993. 199. w14x G. Rauhut, T. Clark, T. Steinke, J. Am. Chem. Soc. 115 Ž1993. 9174. w15x C. Kolle, K.J. Jug, Comput. Chem. 18 Ž1997. 1. w16x A. Klamt, G. Schuurmann, J. Chem. Soc., Perkin Trans. 2 Ž1993. 799. w17x T.N. Truong, E.V. Stefanovich, Chem. Phys. Lett. 240 Ž1995. 253. w18x J. Tomasi, M. Persico, Chem. Rev. 94 Ž1994. 2027. w19x M. Cossi, V. Barone, R. Cammi, J. Tomasi, Chem. Phys. Lett. 255 Ž1996. 327. w20x F.J. Luque, Y. Zhang, C. Aleman, J. Gao, M. Orozco, J. Phys. Chem. 100 Ž1996. 4269. w21x J.B. Foresman, T. Keith, K.B. Wiberg, J. Soonian, M.J. Frisch, J. Phys. Chem. 100 Ž1996. 16098. w22x E. Canes, B. Menucci, J. Tomasi, J. Chem. Phys. 107 Ž1997. 3032. w23x C.M. Cortis, J.-M. Langlois, M.D. Breachy, R.A. Friesner, J. Chem. Phys. 105 Ž1996. 5472. w24x P. Claverie, J.P. Daudey, J. Langlet, B. Pullman, D. Piazzola, J. Phys. Chem. 82 Ž1978. 405. w25x C. Aleman, S.E. Galembeck, Chem. Phys. 232 Ž1998. 151. w26x C. Aleman, Chem. Phys. Lett. 302 Ž1999. 461. w27x P.C. Hariharam, J.A. Pople, Theor. Chim. Acta 28 Ž1978. 213.
162
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
w28x C. Moller, M.S. Plesset, Phys. Rev. 46 Ž1934. 618. w29x S.F. Boys, F. Bernardi, Mol. Phys. 19 Ž1970. 533. w30x F.J. Luque, M.J. Negre, M. Orozco, J. Phys. Chem. 97 Ž1993. 4386. w31x R.A. Pierotti, Chem. Rev. 76 Ž1976. 717. w32x C.J.F. Boetcher, Theory of Electric Polarization, Elsevier, Amsterdam, 1973. w33x Gaussian-94, Revision B3., M.J. Frisch, G.W. Trucks, H.B. Schlegel, P.M.W. Gill, B.G. Johnson, M.A. Robb, J.R. Cheeseman, T.A. Keith, G.A. Petersson, J.A. Montgomery, K. Raghavachari, M.A. Al-Laham, V.G. Zakrzewski, J.V. Ortiz, J.B. Foresman, C.Y. Peng, P.Y. Ayala, W. Chen, M.W. Wong, J.L. Andres, E.S. Replogle, R. Gomperts, R.L. Martin, D.J. Fox, J.S. Binkley, D.J. Defrees, J. Baker, J.J.P. Stewart, M. Head-Gordon, C. Gonzalez, J.A. Pople, Gaussian, Pittsburgh PA, 1995. w34x J.J.P. Stewart, MOPAC 93 Rev. 2, Fujitsu, 1993. w35x L. Gorb, J. Leszczynski, Int. J. Quantum. Chem. 70 Ž1998. 855.
w36x C. Aleman, J.J. Navas, S. Munoz-Guerra, J. Phys. Chem. 99 Ž1995. 17653. w37x J.L. Pascual-Ahuir, E. Silla, E. Tomasi, J. Bonaccorsi, J. Comput. Chem. 8 Ž1987. 778. w38x J. Pitarch, V. Moliner, J.L. Pascual-Ahuir, E. Silla, I. Tunon, ˜ J. Phys. Chem. 100 Ž1996. 9955. w39x C.J. Cramer, D.G. Truhlar, Chem. Phys. Lett. 198 Ž1992. 74. w40x Gao, J. Biophys. Chem. 51 Ž1994. 253. w41x J. Bentzien, J. Florian, T.M. Glennon, A. Warshell, ACS Symposium Series 712, in: J. Gao, M.A. Thompson, ŽEds.., American Chemical Society, New York, 1998, pp. 16–34. w42x J. Florian, J. Sponer, A. Warshell, J. Phys. Chem. B 103 Ž1999. 884. w43x M. Orozco, C. Colominas, F.J. Luque, Chem. Phys. 209 Ž1996. 19. w44x A.H. Elcock, W.G. Richards, J. Am. Chem. Soc. 115 Ž1993. 7930. w45x P.A. Bash, U.C. Singh, R. Langridge, P.A. Kollman, Science 236 Ž1987. 564.
Hydration of cytosine using combined discreterSCRF models: influence of the number of discrete solvent molecules Carlos Aleman ´
)
Departament d’Enginyeria Quımica, E.T.S. d’Enginyers Industrials, UniÕersitat Politecnica de Catalunya, Diagonal 647, Barcelona ´ ` E-08028, Spain Received 11 January 1999
Abstract The solvation of cytosine has been investigated using a combined discreterSCRF model. Molecular structures and gas phase binding energies were obtained at the HFr6-31GŽd. level for ten cytosine–nH 2 O complexes, where n ranges from 1 to 13. Furthermore, gas-phase binding were also computed at the MP2r6-31GŽd. level for complexes with n F 7. The Polarizable Continuum Model self-consistent reaction-field ŽSCRF. were applied to such complexes in order to investigate the effect of the explicit water on the free energy of solvation and on the electronic structure of the solute. The results were compared with those provided by SCRF methods revealing that the discreterSCRF models constitute a very reliable strategy. q 1999 Published by Elsevier Science B.V. All rights reserved.
1. Introduction Most of the current investigations in theoretical chemistry are based on the study of molecules immersed in a solvent phase. The theoretical treatment of the solvent can be broadly classified in two categories: Ži. discrete models; and Žii. continuum models. In the discrete models, solvent molecules are treated explicitly allowing the description of both the microscopic structure of the solvent and the specific solute–solvent interactions. In practice, if the system is described using quantum mechanics, the applicability of the discrete models is restricted to a selected number of configurations of a solute surrounded by a few solvent molecules w1,2x. Conversely, force-field simulations using empirical interaction potentials al)
Tel.: q34-93401-6680r6681; fax: q34-93401-6600r7150; e-mail: [email protected]
low calculations with a large number of both solvent molecules and configurations w3x. However, they cannot account for the modifications induced by the solvent in the electronic structure of the solute. In this connection, hybrid methods including a classical simulation for the solvent terms and solute–solvent interactions coupled to a quantum mechanical description of the solute provide excellent results w4,5x. On the other hand, the continuum models allow the incorporation of solvent effects but neglect the microscopic solvation structure in the vicinity of the solute. Self-consistent reaction-field ŽSCRF. methods based on continuum dielectrics have received much attention within the quantum mechanical framework. A large number of SCRF methods have been developed w6–23x, all of them providing special attention to the determination of the electrostatic contribution to the free energy of solvation Ž DGsol ., which is evaluated assuming that the infinite dielectric
0301-0104r99r$ - see front matter q 1999 Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 Ž 9 9 . 0 0 1 2 2 - 6
152
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
medium, the solvent, reacts against the solute charge distribution generating a reaction field, which in turn interacts with the solute. The main differences between the SCRF methods lie w20x: Ži. in the shape of the solutersolvent interface; Žii. in the definition of the reaction field; and Žiii. in the evaluation of the non-electrostatic contributions to DGsol . For the last few years, SCRF methods have been applied to a wide number of chemical systems w18x. However, electronic effects associated with specific solute–solvent interactions are not considered by these methods. In the middle of the 1970s, Claverie et al. w24x proposed the use of combined discreterSCRF models to overcome this deficiency. M ore recently, we have used com bined discreterSCRF methods to predict both the DGsol and electronic properties in aqueous solution of some natural products w25x and nucleic acid bases w26x. Comparison with the results provided by SCRF methods proved the reliability of discreterSCRF models, which for the study of electronic properties provide more information than either discrete or SCRF methods. However, it should be emphasized that such studies were performed considering a reduced number of explicit solvent molecules, i.e., the maximum number of explicit water molecules was 4 w25x or 5 w26x, depending on the hydrophilicity of the solute. Therefore, it would be interesting to investigate the changes on the electronic properties of the solute as well as in the predicted DGsol on the number of explicit solvent molecules used in the combined discreterSCRF models. In this paper, a systematic study about the influence of the number of explicit solvent molecules in the reliability of combined discreterSCRF models is presented. For this purpose, we investigated the hydration of cytosine nucleic acid base Žsee below., considering complexes with a number of water molecules ranging from 1 to 13. The work is outlined in the following sections. The next section summarises the quantum mechanical methods used in the calculations. Next, we present the results, which are divided in three parts. First, we show the gas-phase binding energies for cytosine–nH 2 O complexes, where n denotes the number of explicit water molecules, computed by the supermolecule approach. Second, we present the solvation of the complexes. In this section, the influence of the num-
ber of explicit water molecules on the DGsol of cytosine is investigated. After this, we analyze the electronic structure of the solute under the influence of the solvent. For this purpose, the atomic charges obtained for different complexes are compared. Finally, the conclusions of the work are presented.
2. Methods 2.1. Gas-phase calculations Molecular structures of four different cytosine– 1H 2 O complexes were optimized at the HFr631GŽd. level w27x. A frequency analysis was performed to verify the nature of the minimum state of the stationary points reached after geometry optimizations. In order to generate the cytosine–nH 2 O complexes with n s 3, 5, 7 and 9, the water molecules were distributed around the solute. First, we distributed the water molecules around the hydrophillic region of cytosine. Once this was surrounded by the first solvation shell, we put explicit solvent molecules in the hydrophobic region of the solute. The complexes were submitted to a complete geometry optimization at the HFr6-31GŽd. level. Single point calculations at the MP2r6-31GŽd. level w28x were performed for all the cytosine–nH 2 O complexes with n s 1, 3, 5 and 7. Cytosine–nH 2 O complexes with n s 11 and 13 were generated using as starting points the optimized complexes with n s 9 and 11, respectively, and adding two new water molecules in the closest positions to the solute. Geometry optimizations of the added water molecules were performed at the HFr6-31GŽd. level but fixing the cytosine and the remaining water molecules.
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
The basis set superposition error ŽBSSE. was corrected at both HF and MP2 levels using the standard counterpoise ŽCP. method w29x. The gasphase binding energy Ž D Egp . was computed as the difference of the total energy of the complex Ž Ecyt – nH2O . and the sum of the solute and nH 2 O energies Ž Ecyt , EnH2O .. D Egp s EcytynH 2 O y Ecyt y EnH 2 O
Ž 1.
153
The reaction field operator was described in terms of an apparent set of point charges ŽEq. Ž4.. spread over the solute–solvent interface, i.e., the solute cavity. The apparent charge density was determined by solving the Laplace equation at the solute–solvent interface. The electrostatic contribution to the DGsol was determined following linear free energy response theory ŽLFER. according to Eq. Ž5. w32x. VR s S i s Ž si . Sir< r y r 0 < s S i qir< r y r 0 <
Ž 4.
DGele s-C < H
y 1r2 Ž - C < V
In a recent work, we compared the reliability of some of the most common SCRF methods in the framework of combined discreterSCRF models w23x. We found that the polarizable continuum model ŽPCM. adapted to the AM1 semi-empirical method ŽPCMrAM1. w20,30x provided excellent results at a very reasonable computational cost. Moreover, the differences between the results provided by the semi-empirical PCMrAM1 and the PCM at the ab initio HFr6-31GŽd. level ŽPCMr6-31GŽd.. were smaller than 1 kcalrmol, which is the error associated to these methods w25x. Therefore, we selected the PCMrAM1 method to describe the bulk water. According to this method, the DGsol was determined as the addition of electrostatic and steric contributions ŽEq. Ž2... The steric contribution was computed as the sum of the cavitation and a van der Waals.
Ž 5.
2.3. Discreter SCRF model The free energy of solvation for the solute with n EW . explicit water molecules Ž DGsol was computed as the addition of the gas-phase binding energy between the solute and the explicit water molecules and the DGsol for the complex ŽEq. Ž6... EW DGsol s D Egp q DGsol Ž cytosine–nH 2 O .
Ž 6.
The cavitation term was determined using Pierotti’s scaled particle theory w31x, while the van der Waals term was evaluated by means of a linear relation with the molecular surface area w20,30x. The PCMrAM1 considers the solvent as a continuum dielectric with the cavity accurately modeled on the solute. The solvent reacts against the solute charge distribution, generating a reaction field Ž VR .. The electrostatic interaction between the solute and the solvent is introduced as a perturbation operator in the solute Hamiltonian ŽEq. Ž3...
It should be noted that in Eq. Ž6., the entropic contributions corresponding to the association of the cytosine and the cluster of n water molecules are neglected in both the gas phase and solution terms. Owing to the chemical chemical nature of the systems under study, i.e., a rigid solute surrounded by the first hydration shell, the entropic contributions are mainly due to the six external degrees of freedom which become internal degrees in the complex. Nevertheless, the entropic contribution in the gas phase and in solution can reasonably be assumed to be very similar, and therefore a considerable compensation is expected. EW Ideally, DGsol should be equivalent to the free energy of solvation of the molecular system comAC . puted in an all-continuum approach Ž DGsol , which is determined as the sum of the free energies of solvation for the solute and the water molecules involved in the molecular system:
Ž H 0 q VR . C s EC
AC s DGsol Ž cytosine. q n )DGsol Ž H 2 O . DGsol
DGsol s DGele q DGcav q DGvdW
Ž 2.
Ž 3.
Ž 7.
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
154
The free energy of solvation for cytosine provided Dr SCRF . by the combined discreterSCRF model Ž DGsol was determined according to Eq. Ž8.: D r SCRF DGsol Ž cytosine. EW s DGsol Ž cytosine–nH 2 O . y n )DGsol Ž H 2 O .
s D Egp q DGsol Ž cytosine–nH 2 O . y n )DGsol Ž H 2 O .
Ž 8. Dr SCRF DGsol
It is worth noting that has an incorrect limit for n ™ ` since D Egp remains constant when a large number of explicit solvent molecules is considered, whereas nUDGsolŽH 2 O. increases faster than DGsolŽcytosine–nH 2 O.. Gas-phase calculations were performed with the Gaussian-94 w33x computer package. PCMrAM1 calculations were performed with an adapted version of the MOPAC 93 w34x computer program 1.
3. Results and discussion 3.1. Gas-phase calculations. The optimized geometries of the four cytosine– 1H 2 O complexes are shown in Fig. 1, where the intermolecular geometric parameters are also displayed. All the complexes are true ground-state minima, which is verified from second derivatives calculations. The molecular geometry of one of these complexes, i.e., that in which the water molecule interacts with the H1 and O2 atoms of cytosine ŽFig. 1Žd.., was recently optimized at the MP2r6-31GŽd. level w35x. The small differences between the intermolecular geometric parameters obtained at the HFr6-31GŽd. and MP2r6-31GŽd. levels indicate that the former is a suitable method of performing geometry optimizations of hydrogen-bonded complexes. More specifically, the donor-acceptor distances predicted at the HF level are overestimated by ˚ with respect to the MP2 values, about 0.05–0.10 A whereas the donor–hydrogen-acceptor angles are underestimated by about 1–28.
1
Adapted to perform PCM calculations by F.J. Luque and M. Orozco
Fig. 2 shows the cytosine–nH 2 O complexes with n s 3, 5, 7, 9, 11 and 13. It is worth noting that water molecules form a network of hydrogen bonds when a large enough number of water molecules is considered. These stable interactions allow the surrounding of the hydrophobic region of cytosine with water molecules in the complexes with n s 11 and 13. Two different types of water molecules can be detected in the cytosine–nH 2 O complexes: Ži. those that form at least one hydrogen bond with the solute; and Žii. those that form hydrogen bonds with other water molecules but not with the solute. Table 1 lists the number of water molecules of each type for the different complexes investigated together with the gas-phase binding energies obtained from Eq. Ž1.. It is worth noting that the maximum number of solvent molecules that are hydrogen bonded to the cytosine, i.e., those that play a crucial role in specific solute– solvent interactions, is six. An inspection of the energies displayed in Table 1 and the molecular geometries depicted in Figs. 1 and 2 indicates that the contributions to the gas-phase binding energies of the water molecules involving one and two hydrogen bonds with cytosine are about 5 and 7–10 kcalrmol, respectively. In contrast, water molecules surrounding the hydrophobic region of cytosine present a weak interaction with the solute, their contribution to the gas-phase binding energy being only 0.5–1.5 kcalrmol. Comparison between the gas-phase binding energies estimated at the HF and MP2 levels ŽTable 1. indicates that the two methods provide very similar values from a qualitative point of view. Accordingly, the binding energy order predicted by the two methods for the monohydrated complexes is the same, the HF binding energies being systematically less attractive than those for MP2. Such underestimation corresponds, on average, to an error of about 8%. Similar results were obtained in previous works for related hydrogen bonded systems w25,26,36x. On the other hand, it is worth noting that the relative energy order predicted for the four monohydrated complexes can easily be rationalized considering the hydrogen bonding parameters Žsee Fig. 1.. A detailed analysis of both cytosine–water and water–water hydrogen bonding interactions for cytosine–nH 2 O complexes is presented in Table 2, where the hydrogen-acceptor distances and the donor–hy-
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
155
Fig. 1. Molecular geometries of the four cytosine–1H 2 O complexes optimized at the HFr6-31GŽd. level. The hydrogen bonding parameters are also displayed.
drogen-acceptor angles are listed for the complexes with n s 3, 5, 7 and 9. The hydrogen bonding ˚ larger for the cytosine– distances are about 0.1 A water interactions than for the water–water interactions. On the other hand, an inspection of the complexes with n s 7 and 9 reveals that the donor–hydrogen-acceptor angles are less favorable for the water–water interactions than for the cytosine–water interactions. Similar conclusions are reached from the analysis of the complexes with n s 11 and 13. 3.2. Aqueous solution calculations The DGsol values predicted by the PCMrAM1 method for all the complexes investigated are listed in Table 3, where both the electrostatic and steric
contributions are also displayed. It is worth nothing that the electrostatic contribution is the most important contribution to DGsol . This attractive contribution mainly depends on the accessibility to the bulk solvent of the polar atoms contained in both the cytosine and explicit water molecules. Accordingly, the electrostatic contribution increases with the surface accessible to the bulk solvent. This feature is clearly evidenced in the monohydrated complexes. Thus, in the complex displayed in Fig. 1Žc., the two hydrogen atoms contained in the water molecule are accessible to the bulk solvent, whereas in the other three monohydrated complexes only one atom of the water molecule points to the bulk solvent. Therefore, the attractive electrostatic interactions are larger for the complex displayed in Fig. 1Žc. than for the
156
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
Fig. 2. Molecular geometries of the cytosine–nH 2 O complexes optimized at the HFr6-31GŽd. level: Ža. n s 3; Žb. n s 5; Žc. n s 7; Žd. n s 9; Že. n s 11 and Žf. n s 13.
others. On the other hand, the steric term provides a repulsive contribution since the polarizable nature of
water makes this solvent to be highly structured, and therefore, generation of cavity is difficult from an
C. Alemanr ´ Chemical Physics 244 (1999) 151–162 Table 1 Number of water molecules that are hydrogen bonded to cytosinea Ž NHB S . and to other water moleculesb Ž NHBW . in cytosine– nH 2 O complexes. Gas-phase binding energies c ŽEq. Ž1.. computed at the HFr6-31GŽd. and MP2r6-31GŽd. levels of theory Complex d
NHB S NHBW D EHFr6-31GŽd. D E MP2r6-31GŽd.
Cytosine–1H 2 O Ža. Cytosine–1H 2 O Žb. Cytosine–1H 2 O Žc. Cytosine–1H 2 O Žd. Cytosine–3H 2 O Cytosine–5H 2 O Cytosine–7H 2 O Cytosine–9H 2 O Cytosine–11H 2 O Cytosine–13H 2 O
1 1 1 1 3 4 5 6 6 6
– – – – – 1 2 3 5 7
y7.1 y9.3 y5.4 y9.8 y20.7 y31.0 y36.1 y46.0 y47.0 y50.1
y7.3 y10.3 y6.2 y10.4 y22.8 y34.1 y39.6 – – –
157
Table 3 Free energies of solvationa predicted for the cytosine– nH 2 O complexes by the PCMrAM1 method. Electrostatic and steric Ž DGcav q DGvdW . contributionsa are included. Complex
DGsol
DGele
DGcav q DGvdW
Cytosine–1H 2 O Ža. Cytosine–1H 2 O Žb. Cytosine–1H 2 O Žc. Cytosine–1H 2 O Žd. Cytosine–3H 2 O Cytosine–5H 2 O Cytosine–7H 2 O Cytosine–9H 2 O Cytosine–11H 2 O Cytosine–13H 2 O
y16.3 y16.1 y21.0 y17.3 y19.8 y18.3 y21.4 y23.8 y30.5 y35.8
y20.5 y20.3 y25.2 y21.6 y25.3 y25.2 y29.3 y32.8 y40.8 y47.3
4.2 4.2 4.2 4.3 5.5 6.9 7.9 9.0 10.3 11.5
a
In kcalrmol.
a
Number of water molecules with one or more hydrogen bonding interactions with cytosine. b Number of water molecules that are hydrogen bonded to other water molecules but not cytosine. c In kcalrmol. d The labels of the cytosine–1H 2 O complexes correspond to those displayed in Fig. 1.
energy point of view. The repulsive cavitation contribution is not compensated by the van der Waals interactions, which, due to the small size of the cavity, are not attractive enough.
Table 2 Cytosine . . . water a and water . . . water hydrogen bonding parametersb measured in the optimized molecular geometries of the cytosine–nH 2 O complexes with n s 3, 5, 7 and 9 Complex
Cytosine–3H 2 O
Cytosine–5H 2 O
Cytosine–7H 2 O
Cytosine–9H 2 O
a
Cytosine . . . water
Water . . . water
Distances
Angles
Distances
Angles
O8 . . . Hw s 2.13 H12 . . . 6Ow s 2.04 H13 . . . Ow s 2.11 O8 . . . Hw s 2.00 H11 . . . Ow s 2.08 N3 . . . Hw s 1.99 H12 . . . Ow s 2.17 H13 . . . Ow s 2.03 O8 . . . Hw s 1.90 H11 . . . Ow s 1.96 N3 . . . Hw s 2.01 H12 . . . Ow s 2.10 H13 . . . Ow s 2.11 O8 . . . Hw s 1.91 O8 . . . Hw s 1.94 H11 . . . Ow s 2.11 N3 . . . Hw s 2.07 H12 . . . Ow s 2.13 H13 . . . Ow s 2.07
O8 . . . Hw–Ow s 166.6 N7–H12 . . . Ow s 161.8 N7–H12 . . . Ow s 172.5 O8–Hw–Ow s 147.9 N1–H11 . . . Ow s 146.5 N3 . . . Hw–Ow s 165.8 N7–H12 . . . Ow s 167.6 N7–H13 . . . Ow s 153.3 O8 . . . Hw–Ow s 166.3 N1–H11 . . . Ow s 174.6 N3 . . . Hw–Ow s 165.8 N7–H12 . . . Ow s 173.1 N7–H13 . . . Ow s 146.0 O8 . . . Hw–Ow s 167.8 O8 . . . Hw–Ow s 178.8 N1–H11 . . . Ow s 177.9 N3 . . . Hw–Ow s 169.6 N7–H12 . . . Ow s 167.2 N7–H13 . . . Ow s 148.8
2.00 – – 1.84 1.94 1.93 – – 1.88 1.85 2.00 1.86 2.12 1.85 1.97 1.83 1.99 1.85 2.13
164.2 – – 159.3 173.5 167.9 – – 159.1 156.9 165.0 159.9 147.5 157.6 160.7 154.7 165.8 160.3 147.6
The numeration for the atoms of cytosine is displayed in the Scheme I while H W and O W refer to the hydrogen and oxygen atoms of the water molecules. b ˚ and degrees, respectively. Hydrogen-acceptor distances and donor–hydrogen-acceptor angles are in A
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
158
AC EW Table 4 shows the values of DGsol and DGsol computed according to Eqs. Ž6. and Ž7., respectively, as well as the difference between them, which can be considered as a measure of the reliability of the combined discreterSCRF model. As expected from EW estimated the results displayed in Table 1, the DGsol from the HFr6-31GŽd. and MP2r6-31GŽd. gasphase binding energies are very similar. FurtherAC more, the agreement between DGsol and both EW EW Ž DGsol HF. and DGsol ŽMP2. is encouraging, the maximum difference between these parameters being 12.2 and 5.5 kcalrmol, respectively. Such differences correspond to an error of 12.4% and 14.8% for the cytosine–13H 2 O and cytosine–3H 2 O complexes, respectively. Fig. 3 represents the variation of the error Žin %. corresponding to the difference AC EW Ž between DGsol and both DGsol HF . and EW Ž DGsol MP2. of the number of explicit solvent molecules. Inspection to the data derived from HFr6-31GŽd. gas-phase binding energies indicates that good results are achieved for cytosine–nH 2 O complexes with n s 5, 7 and 9, i.e., the error is lower than 6.5%, whereas the error increases when the number of explicit water molecules is either too
Table 4 Free energies of solvationa,b predicted using an all continuum AC . EW . approachc Ž DGsol and explicit solvent molecules Ž DGsol for AC cytosine– nH 2 O complexes. The differences b w DGsol y AC EW Ž EW Ž DGsol HF.x and w DGsol y DGsol MP2.x are listed in parenthesis in the 2nd and 3rd columns, respectively Complex
AC DGsol
EW Ž DGsol HF. d
EW Ž DGsol MP2. d
Cytosine–1H 2 O Ža. Cytosine–1H 2 O Žb. Cytosine–1H 2 O Žc. Cytosine–1H 2 O Žd. Cytosine–3H 2 O Cytosine–5H 2 O Cytosine–7H 2 O Cytosine–9H 2 O Cytosine–11H 2 O Cytosine–13H 2 O
y24.9 y24.9 y24.9 y24.9 y37.1 y49.3 y61.5 y73.7 y85.9 y98.1
y23.4 ( y 1.5) y25.4 (0.5) y26.4 (1.5) y27.2 (2.2) y40.5 (3.4) y49.3 (0.0) y57.5 ( y 4.0) y69.8 ( y 3.9) y77.5 ( y 8.4) y85.9 ( y 12.2)
y23.6 ( y 1.3) y26.4 (1.5) y27.2 (2.3) y27.7 (2.8) y42.6 (5.5) y52.4 (3.1) y61.0 ( y 0.5) – – –
a
AC EW DGsol and DGsol computed according to Eqs. Ž6. and Ž7., respectively. b In kcalrmol. c The free energy of solvation of the water molecule used to AC compute DGsol is y6.1 kcalrmol. d EW Ž EW Ž DGsol HF. and DGsol MP2. refer to the values computed determined from HFr6-31GŽd. and MP2r6-31GŽd. gas-phase binding energies, respectively.
Fig. 3. Variation of the error corresponding to the difference AC EW Ž EW Ž between DGsol and both DGsol HF. ŽI. and DGsol MP2. ŽB. upon the number of explicit water molecules.
large or too small. The errors corresponding to the AC EW Ž differences between DGsol and DGsol MP2. suggest similar trends. These results indicate that the strategy of including explicit water molecules is very reliable when the number of solvent molecules included in the calculation is sufficient to surround the hydrophillic region of the solute. In contrast, when the number of water molecules is too large, such as in cytosine– 11H 2 O and cytosine–13H 2 O complexes, the error can be greater than 12%. Table 1 shows that these complexes contain an important number of explicit water molecules without cytosine–water hydrogen bonding interactions Ž NHB W s 5 and 7, respectively.. Inspection of Fig. 2Že. and 2Žf. reveals that these water molecules are mainly located in the hydrophobic region of the solute or behind the water molecules hydrogen-bonded to the cytosine. On the other hand, when the number of explicit solvent molecules is too small, such as in cytosine–1H 2 O and cytosine–3H 2 O complexes, the error can be greater than 9%. In these cases, the number of water molecules is not sufficient to surround the hydrophillic region of cytosine completely ŽFigs. 1 and 2Ža... The possible errors in the combined discreter SCRF model can mainly be attributed to four different sources. The first source is the complexity of the molecular surface in complexes. The building-up of the surface, as determined in PCMrAM1 w37x, defines a large number of small spheres not assigned specifically to any atom, but needed to describe
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
properly a smooth solute–solvent interface and, consequently, to estimate correctly the solvent reaction field. The number of such spheres increases with the complexity of the surface, providing a less accurate prediction of DGsol . The second source of error is related to the use of solvent-accessible surfaces which are not the most appropriate to describe aggregation processes in solution. Thus, previous studies of intermolecular interactions in solution indicated that the solvent-accessible surface is not able to describe properly the area of the cavity not accessible to the bulk solvent w38x. Furthermore, the error associated with this surface model increases with the contact distance until the rupture of the complex. Inspection of Fig. 2 and Table 4 indicates that the errors associated to these deficiencies are minimized when the number of water molecules is sufficient to fill up the hydrophillic region of cytosine, since the number of void spaces in all the possible cytosine–water and water–water dimers is very small. However, the complexity of the surface increases when the number of water molecules is either too small or too large. Furthermore, if the number of water molecules is too large, then the contact distance in some of the cytosine–water and water–water dimers involved increases and the worst description of the volume inaccessible to the bulk solvent is obtained. The third
159
source of error is the neglection of the vibrational and rotational nuclear movements. As was stated in Section 2, the vibrational contribution is solution is compensated by that of the gas phase. However, the rotational nuclear motions provide an unfavorable contribution to DGsol due to the frictional solvent effects on the rotation of the solute system. It should be emphasized that the rotational nuclear motions do not have a direct effect on the solute charge distribution and have been included in the parametrization of the nonelectrostatic terms of the PCMrAM1 method, However, the parameters were developed for single solutes and not for complexes. Finally, it should be noted that we are considering frozen geometries. This approximation is not very important for rigid solutes like cytosine but could lead to dramatic errors in the study of flexible compounds. Dr SCRF The values of DGsol predicted for cytosine ŽEq. Ž8.. using the complexes with 5, 7 and 9 explicit water molecules are displayed in Table 5. The DGsol values estimated for cytosine w39–42x and its methylated derivative w39,43–45x using other computational approaches are also included in the Dr SCRF values derived table for comparison. The DGsol from HFr6-31GŽd. gas-phase binding energies are lower than those obtained from MP2r6-31GŽd. energies. They range from y14.8 to y18.8 kcalrmol
Table 5 Free energy of solvationa Ž DGsol . predicted for cytosine using the combined discreterSCRF methodb . DGsol values predicted by other computational approaches for cytosine and its methylated derivative are also included for comparison DGsol
Method
Computational approach
DiscreterSCRF DiscreterSCRF DiscreterSCRF DiscreterSCRF DiscreterSCRF SCRF SCRF LD LD MC QMrMM SCRF SCRF SCRF OPLSrFEP OPLSrAMBER
Cytosine–5H 2 O complex using HFr6-31GŽd. gas-phase binding energy and PCMrAM1. y18.8 Cytosine–5H 2 O complex using MP2r6-31GŽd. gas-phase binding energy and PCMrAM1. y21.9 Cytosine–7H 2 O complex using HFr6-31GŽd. gas-phase binding energy and PCMrAM1. y14.8 Cytosine–7H 2 O complex and MP2r6-31GŽd. gas-phase binding energy and PCMrAM1 y18.3 Cytosine–9H 2 O complex and HFr6-31GŽd. gas-phase binding energy and PCMrAM1 y14.9 Cytosine using PCMrAM1 y18.8 Cytosine using SM2rAM1 y20.6 Cytosine using the noniterative LD model y14.2 Cytosine using the iterative LD model y18.0 Cytosine using mixed quantum mechanics and classical Monte Carlo simulations y16.3 Methylated derivative of cytosine using PCMrAM1 y16.1 Methylated derivative of cytosine using PCMr6-31GŽd. y15.1 Methylated derivative of cytosine using SM2rAM1 y18.7 Methylated derivative of cytosine using free energy preturbations with the OPLS force-field y20.1 Methylated derivative of cytosine using free energy preturbations with the AMBER force-field y12.7
a
In kcalrmol. DrSCRF DGsol values predicted using Eq. Ž8..
b
Ref. Present work Present work Present work Present work Present work Present work w39x w41x w42x w40x w43x w43x w39x w44x w45x
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
160
and from y18.3 to y21.9 kcalrmol, respectively. Dr SCRF On the other hand, the DGsol values predicted for cytosine using 11 and 13 explicit solvent molecules are y10.4 and y6.6 kcalrmol, respectively. These low values are mainly due to the underestimation of the gas-phase binding energies estimated at the HF level. Dr SCRF The DGsol values displayed in Table 5, i.e., those obtained using the 5, 7 and 9 water molecules, are within the range of previous estimations predicted by other computational approaches. Thus, the DGsol values predicted by SM2, PCM and iterative LD models Žy20.6, y18.8 and y18.0 kcalrmol, respectively. are very close to those obtained using the MP2r6-31GŽd. gas-phase energies. On the other hand, the value obtained from Monte Carlo quantum mechanical and molecular mechanical ŽMC QMrMM. simulations Žy16.3 kcalrmol. is close to Dr SCRF the DGsol values estimated from HFr6-31GŽd. gas-phase energies. These results suggest that the values predicted by both combined discreterSCRF model using HFr6-31GŽd. gas-phase binding energies and molecular mechanics methods are slightly underestimated. Unfortunately, there is not a quantitative experimental DGsol value for cytosine. It should be emphasized that the solvation studies performed for its methylated derivative w39,43–45x also provide a wide range of DGsol values. Table 6 shows the Mulliken charges of isolated cytosine in a gas phase provided by the AM1 hamil-
tonian, isolated cytosine in an aqueous solution obtained using the SCRF method and cytosine in an aqueous solution predicted by the combined discreterSCRF method for the complexes with 3, 7 and 13 explicit water molecules. Comparison between the charges obtained in the gas phase and aqueous solution points out the complexity of the polarizing effect of water on the solute charge distribution. The more important features are the enhancement of the negative charges of the O2 and N3 atoms, as well as the increase of the positive charge of the hydrogen atoms bonded to N1 and N4. Furthermore, there is also an appreciable increase in the charge of the hydrogen atoms located in the hydrophobic region of the cytosine. Comparison of the atomic charges derived for cytosine using the combined discreterSCRF method reveals very similar values for all the complexes investigated. This trend suggests that the number of explicit water molecules is not crucial to obtain a suitable description of the electronic effects induced by specific solute–solvent interactions. For instance, the largest difference between the charges obtained for the cytosine–3H 2 O and cytosine–13H 2 O complexes is about 0.02 e Žwhere e is the charge of an electron. for the O2 atom. Of particular interest is the comparison of the atomic charges provided by the combined discreter SCRF model and those obtained using the SCRF method without including any explicit solvent
Table 6 Mulliken atomic chargesa for cytosine computed in the gas phase and in aqueous solution. Charges in aqueous solution were obtained using the PCMrAM1 method and the combined discreterSCRF method considering 3, 7 and 13 explicit water molecules. Charges obtained for the cytosine–7H 2 O complex but without the SCRF are also listed Žin italics. Atom
Cytosine Žgas phase.
Cytosine ŽPCMrAM1.
Cytosine–3H 2 O
Cytosine–7H 2 O
Cytosine–13H 2 O
N1 C2 N3 C4 C5 C6 H– ŽN1. O5ŽC2. N4 H– ŽC5. H– ŽC6. H–N3 H–N5
y0.324 0.349 y0.299 0.237 y0.348 0.065 0.252 y0.344 y0.372 0.153 0.152 0.246 0.232
y0.308 0.364 y0.423 0.255 y0.348 0.086 0.301 y0.503 y0.349 0.179 0.186 0.278 0.281
y0.324 0.395 y0.443 0.266 y0.344 0.087 0.287 y0.479 y0.341 0.176 0.174 0.263 0.277
y0.327 (y0.334) 0.400 (0.396) y0.451 (y0.410) 0.273 (0.267) y0.345 (y0.351) 0.090 (0.081) 0.275 (0.273) y0.461 (y0.396) y0.342 (y0.362) 0.173 (0.159) 0.173 (0.154) 0.262 (0.260) 0.273 (0.261)
y0.320 0.401 y0.455 0.272 y0.347 0.092 0.275 y0.458 y0.345 0.175 0.184 0.259 0.274
a
In units of electron.
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
molecule. It is worth noting that charges provided the SCRF method are somewhat enhanced for some atoms with respect to those obtained for the cytosine–nH 2 O complexes. For instance, atomic charges predicted for C2, N3, H1 and O2 differ by 0.036, 0.028, 0.026 and 0.042 e, respectively, with respect to those obtained for the cytosine–7H 2 O complex. On the other hand, Table 6 also shows the atomic charges for cytosine–7H 2 O complex but without applying the SCRF, i.e., considering a discrete method. It is worth noting that the atomic charges obtained from the discrete and discreterSCRF methods are quite similar for all the atoms except for N3 and O2, which differ by 0.041 and 0.065 e, respectively. However, an inspection to the results listed in Table 6 indicates that polarizing effect in cytosine seems to be overestimated by SCRF methods.
4. Conclusions In this study, we investigated the influence of the number of explicit water molecules in the solvation of cytosine using a combined discreterSCRF model. For this purpose, ten cytosine–nH 2 O complexes, where n ranges from 1 to 13, were optimized at the HFr6-31GŽd. level. The gas-phase binding energies were computed at both HFr6-31GŽd. and MP2r631GŽd. levels of theory. A combined discreterSCRF model using the PCMrAM1 method was applied to compute the DGsol of the ten complexes. From the AC EW comparison of DGsol and DGsol , we found that the best results, i.e., the smallest error between the two magnitudes, are obtained when the explicit water molecules are distributed around the hydrophillic region of cytosine filling up all the solvation sites. On the other hand, when the number of water molecules is too small or too large both the complexity of the molecular surface and the use of solventaccessible surfaces increase the differences between AC EW DGsol and DGsol . The analysis of the change in the atomic charges of cytosine provided by the combined discreterSCRF model emphasizes the important role of the specific interactions between the solute and the solvent. Furthermore, the results indicate that SCRF methods tend to overestimate the polarizing effect of the solvent on the solute suggesting that more investigation is needed in this area.
161
Acknowledgements Author is indebted to Dr. M. Orozco and F.J. Luque for making available to him their version of MOPAC93 adapted to perform PCM calculations. Author thanks to CESCA for computational facilities. References w1x O. Tapia, F. Colonna, J.G. Angyan, J. Chem. Phys. 87 Ž1990. 875. w2x G. Alagona, A. Pullman, E. Scrocco, J. Tomasi, Int. J. Pept. Protein Res. 5 Ž1973. 251. w3x M.P. Allen, D.J. Tildesley, in: Computer Simulation of Liquids, Oxford University Press, Oxford, 1987. w4x M.J. Field, P.A. Basch, M. Karplus, J. Comput. Chem. 11 Ž1990. 700. w5x J. Gao, J. Am. Chem. Soc. 116 Ž1994. 1563. w6x L. Rivail, D. Rinaldi, Theor. Chim. Acta 32 Ž1973. 57. w7x M.W. Wong, M.J. Frisch, K.B. Wiberg, J. Am. Chem. Soc. 113 Ž1991. 4776. w8x S. Miertus, E. Scrocco, Tomasi, J. Chem. Phys. 55 Ž1981. 117. w9x M.M. Karelson, T. Tam, A.R. Katritzky, S.J. Cato, M. Zerner, Tetrahedron Comput. Methodol. 2 Ž1989. 295. w10x B. Wang, G.P. Ford, J. Comput. Chem. 7 Ž1992. 4162. w11x G.E. Chudinov, D.V. Napolov, M.U. Basilevsky, Chem. Phys. 160 Ž1992. 41. w12x C.J. Cramer, D.G. Truhlar, Science 256 Ž1992. 213. w13x V. Luzhkov, A.J. Warshel, Comput. Chem. 13 Ž1993. 199. w14x G. Rauhut, T. Clark, T. Steinke, J. Am. Chem. Soc. 115 Ž1993. 9174. w15x C. Kolle, K.J. Jug, Comput. Chem. 18 Ž1997. 1. w16x A. Klamt, G. Schuurmann, J. Chem. Soc., Perkin Trans. 2 Ž1993. 799. w17x T.N. Truong, E.V. Stefanovich, Chem. Phys. Lett. 240 Ž1995. 253. w18x J. Tomasi, M. Persico, Chem. Rev. 94 Ž1994. 2027. w19x M. Cossi, V. Barone, R. Cammi, J. Tomasi, Chem. Phys. Lett. 255 Ž1996. 327. w20x F.J. Luque, Y. Zhang, C. Aleman, J. Gao, M. Orozco, J. Phys. Chem. 100 Ž1996. 4269. w21x J.B. Foresman, T. Keith, K.B. Wiberg, J. Soonian, M.J. Frisch, J. Phys. Chem. 100 Ž1996. 16098. w22x E. Canes, B. Menucci, J. Tomasi, J. Chem. Phys. 107 Ž1997. 3032. w23x C.M. Cortis, J.-M. Langlois, M.D. Breachy, R.A. Friesner, J. Chem. Phys. 105 Ž1996. 5472. w24x P. Claverie, J.P. Daudey, J. Langlet, B. Pullman, D. Piazzola, J. Phys. Chem. 82 Ž1978. 405. w25x C. Aleman, S.E. Galembeck, Chem. Phys. 232 Ž1998. 151. w26x C. Aleman, Chem. Phys. Lett. 302 Ž1999. 461. w27x P.C. Hariharam, J.A. Pople, Theor. Chim. Acta 28 Ž1978. 213.
162
C. Alemanr ´ Chemical Physics 244 (1999) 151–162
w28x C. Moller, M.S. Plesset, Phys. Rev. 46 Ž1934. 618. w29x S.F. Boys, F. Bernardi, Mol. Phys. 19 Ž1970. 533. w30x F.J. Luque, M.J. Negre, M. Orozco, J. Phys. Chem. 97 Ž1993. 4386. w31x R.A. Pierotti, Chem. Rev. 76 Ž1976. 717. w32x C.J.F. Boetcher, Theory of Electric Polarization, Elsevier, Amsterdam, 1973. w33x Gaussian-94, Revision B3., M.J. Frisch, G.W. Trucks, H.B. Schlegel, P.M.W. Gill, B.G. Johnson, M.A. Robb, J.R. Cheeseman, T.A. Keith, G.A. Petersson, J.A. Montgomery, K. Raghavachari, M.A. Al-Laham, V.G. Zakrzewski, J.V. Ortiz, J.B. Foresman, C.Y. Peng, P.Y. Ayala, W. Chen, M.W. Wong, J.L. Andres, E.S. Replogle, R. Gomperts, R.L. Martin, D.J. Fox, J.S. Binkley, D.J. Defrees, J. Baker, J.J.P. Stewart, M. Head-Gordon, C. Gonzalez, J.A. Pople, Gaussian, Pittsburgh PA, 1995. w34x J.J.P. Stewart, MOPAC 93 Rev. 2, Fujitsu, 1993. w35x L. Gorb, J. Leszczynski, Int. J. Quantum. Chem. 70 Ž1998. 855.
w36x C. Aleman, J.J. Navas, S. Munoz-Guerra, J. Phys. Chem. 99 Ž1995. 17653. w37x J.L. Pascual-Ahuir, E. Silla, E. Tomasi, J. Bonaccorsi, J. Comput. Chem. 8 Ž1987. 778. w38x J. Pitarch, V. Moliner, J.L. Pascual-Ahuir, E. Silla, I. Tunon, ˜ J. Phys. Chem. 100 Ž1996. 9955. w39x C.J. Cramer, D.G. Truhlar, Chem. Phys. Lett. 198 Ž1992. 74. w40x Gao, J. Biophys. Chem. 51 Ž1994. 253. w41x J. Bentzien, J. Florian, T.M. Glennon, A. Warshell, ACS Symposium Series 712, in: J. Gao, M.A. Thompson, ŽEds.., American Chemical Society, New York, 1998, pp. 16–34. w42x J. Florian, J. Sponer, A. Warshell, J. Phys. Chem. B 103 Ž1999. 884. w43x M. Orozco, C. Colominas, F.J. Luque, Chem. Phys. 209 Ž1996. 19. w44x A.H. Elcock, W.G. Richards, J. Am. Chem. Soc. 115 Ž1993. 7930. w45x P.A. Bash, U.C. Singh, R. Langridge, P.A. Kollman, Science 236 Ž1987. 564.