J. Phys. Chem. B 2001, 105, 5827-5832
5827
Electronic Polarization and Hydration of the Dimethyl phosphate Anion: An ab Initio Molecular Dynamics Study I-Feng Kuo and Douglas J. Tobias* Department of Chemistry, UniVersity of California, IrVine, California 92697-2025 ReceiVed: October 23, 2000; In Final Form: March 23, 2001
We report a molecular dynamics (MD) simulation study of sodium dimethlyphosphate in aqueous solution. The dimethyl phosphate anion (DMP) is a model compound for the phosphodiester linkage in phospholipids and nucleic acids. We compare results from classical MD simulations based on an empirical force field to an ab initio MD simulation in which the forces are derived from the electronic structure in the framework of density functional theory with the Car-Parrinello approach. We have characterized the molecular and electronic structure of both the DMP anion and its first solvation shell in detail. To analyze the charge distribution in individual molecules, we have used maximally localized Wannier orbitals. The radial and orientational distributions of the water molecules revealed a softening of the structure of the first solvation shell of the DMP anion in the ab initio simulation compared to the force field simulation. The ab initio simulation reveals that both the geometry and the electronic structure of the DMP anion change significantly on moving from the gas phase to aqueous solution. In contrast, the dipole moments of the water molecules are comparable to bulk water, even in the vicinity of the DMP anion and the sodium counterion.
1. Introduction The dimethyl phosphate anion (DMP, Figure 1) is a model compound for the phosphodiester linkage found in the highly charged backbone of nucleic acids and the zwitterionic headgroups of phospholipids. Theoretical investigations of the structure and energetics of DMP have employed Monte Carlo1-3 simulations based on empirical force fields.4-8 Although the early efforts were mainly focused on the two phosphodiester torsion angles with implications for the conformational flexibility of nucleic acids, attention was also directed toward the effect of hydration along with the influence of counterions.6-8 From these theoretical studies, the hydration of DMP was described in terms of three distinct hydration zones: ionic, hydrophilic, and hydrophobic regions.4,6 At the same time, the results also indicated a stabilization of different conformers due to the aqueous environments, in agreement with experimental data.7,8 Furthermore, the geometrical arrangement of water around the different regions of DMP is remarkably variable, with the most interesting results appearing in the highly ionic regions between the partial negatively charged oxygens in DMP.4,9 Around this ionic region, the water exhibits a high degree of orientational order, with a preferential orientation of a hydrogen atom pointing toward the nonester oxygen.4 Gas-phase studies have also been carried out on DMP with sophisticated quantum chemical calculations based on HartreeFock methods as well as density functional approaches.5,9,10 From these electronic structure calculations, precise relative energies of the various conformations of DMP have been determined, and the reduction of energetic barriers between different conformations in condensed phase have been evaluated with the use of polarized continuum solvation models.5 The inclusion of explicit solvent molecules in electronic structure calculations further confirmed diminished conformational preferences of DMP in aqueous solution. With the use of one explicit H2O placed between the two nonester oxygen on DMP, it was
Figure 1. Molecular structure of the DMP anion with standard atom names. In the text the nonesterified (lone) OA and OB atoms are frequently referred together as the O* atoms, and the esterified O3′ and O5′ as the O atoms.
found that even monohydrated DMP is expected to have an increase in accessible conformational space.9 With the advent of ab initio molecular dynamics (MD) simulation based on density functional theory,11 a new level of accuracy can be achieved in the theoretical prediction of the properties of aqueous solutions.12 Through the use of ab initio MD, quantities such as the radial distribution function and diffusion coefficient of water have been reproduced to within experimental accuracy, and the dipole moment of a water molecule in the liquid has been predicted.13,14 With this new tool, it was found that environmental effects can change the nature of the chemical bonds in DMP.15 Furthermore, with the trajectories obtained from ab initio MD simulations, the power spectrum was obtained for a protonated dimethyl phosphate (HDMP).16 The results of this calculation confirmed the importance of environmental effects in reshaping the potential energy surface of HDMP. A recent optimization of the sodium GpC nonahydrate crystal structure has demonstrated the potential utility of density functional electronic structure theory in highly accurate simulations of solvated biomolecules.17 One question that arises during the study of molecular fragments such as DMP is how applicable are the results obtained from simulation of molecular fragments to macromolecules such as nucleic acids or phospholipids? Recent simulation results for the radial distributions of water around the phosphate group in phospholipids18 are in agreement with previous Monte
10.1021/jp003900a CCC: $20.00 © 2001 American Chemical Society Published on Web 05/25/2001
5828 J. Phys. Chem. B, Vol. 105, No. 24, 2001
Tobias and Kuo
Carlo results for DMP. Furthermore, a three-dimensional water density distribution displays preferential solvation sites around the phosphate that are also seen in the present study of DMP.18 These and other examples demonstrate that it is reasonable to study molecular fragments such as DMP in order to gain insights into large biomolecules that are impossible to study with the same detail. In the present paper, we add a new dimension to previous studies of the DMP anion by examining for the first time the role that explicit electronic polarization plays in the molecular and electronic structure of DMP and its first hydration shell, through the use of both ab initio and nonpolarizable force field based MD simulations. 2. Computational Details A classical MD simulation was performed with an empirical force field for 2 ns under constant atmospheric pressure at a constant temperature of 300 K starting with a DMP anion and a sodium counterion in a cubic box of length 20 Å containing 221 H2O with periodic boundary conditions. Next, the size of the system was reduced by a repetitive trimming procedure, in which progressively smaller systems were produced and equilibrated at constant temperature and pressure for 2 ns at each stage. The final system generated consisted of 23 H2O, the DMP anion, and the sodium counterion in a cubic box with an edge length of 9.8 Å. Following the constant pressure and temperature equilibration, the final system was simulated at constant volume for 2 ns to generate data for comparison with the ab initio MD trajectory. The sodium ion remained well separated from the DMP anion throughout all of the simulations. The equations of motion were integrated using reversible, multiple time step, extended system algorithms19 with a time step of 3 fs. An improved version9 of the CHARMM 22 nucleic acid force field20 optimized for solute-solvent and solvent-solvent interactions was employed for DMP in conjunction with the TIP3P21 water model. The electrostatic interactions were computed using the particle mesh Ewald method.22 All of the force field based calculations were carried out using the PINY_MD program.23 The ab initio MD trajectory was generated using the CPMD program24 starting from the last spatial configuration from the force field based simulation of the system with 23 water molecules, in a cubic box with an edge length of 9.8 Å and periodic boundary conditions. The equations of motion for the electrons and ions were integrated using the method of Car and Parrinello,11 with a fictitious electron mass of 800 au, and a time step of 0.121 fs. The electronic structure was computed within the Kohn-Sham formulation of density functional theory,25,26 with the gradient corrected BLYP27,28 exchange and correlation functional, which has been shown to provide an accurate description of aqueous systems.29 Only the valence electrons were treated explicitly, and the valence-core interactions were described by norm-conserving pseudopotentials of the Troullier-Martins form30 for C, O, and H. For the sodium ion, the 2s22p6 electrons were treated explicitly, and the valencecore interactions were handled by a separable dual-space Gaussian pseudopotential.31 The Kohn-Sham orbitals were expanded in a plane-wave basis set up to an energy cutoff of 70 Rydberg. After equilibration, a 4.2 ps trajectory was used for analysis. The average ionic temperature was 300 K for the duration of the run. The ab initio simulation required approximately 47 h per ps using 12 CPUs and 533 MB of memory on a Silicon Graphics Origin 2000 computer. The optimized structure of the DMP anion in the gas phase was obtained by energy minimization using the same electronic structure description as the ab initio MD simulation. The ion
Figure 2. Radial distribution functions for water oxygen atoms around the nonesterified oxygen atoms of a dimethyl phosphate anion in a series of force field based MD simulations of different sized systems.
was placed in a cubic box of length 16.25 Å, and cluster boundary conditions were applied using the method of Martyna and Tuckerman.32 The optimization was stopped when the rootmean-squared gradient was less than 10-3. We have employed a recent advancement in the analysis of electronic structure in a periodic system using a plane wave basis that involves the transformation of the Kohn-Sham orbitals into maximally localized Wannier functions,33-36 as implemented in the CPMD program.37 In this way, the total charge distribution in a system can be partitioned into individual molecular contributions (e.g., bonding pairs and lone pairs). The coordinates of these centers of these localized orbitals, the Wannier function centers (WFCs), can then be used in electrostatic calculations, for example, to describe polarization in terms of fluctuations of the dipole moments. Moreover, it is also possible to use the WFCs to analyze changes in chemical bonding during an ab initio MD simulation, thereby providing additional insight not obtainable from force field based calculations.15 3. Results and Discussion Force Field Based Simulation Results and Finite Size Effects. A small system containing only 23 water molecules was employed to make the ab initio MD simulation tractable on our computing resources. This system size was chosen based on the results of the trimming procedure, in which we established using force field based simulations that the 23-water system faithfully reproduced the structure of the first solvation shell in the largest system considered (containing 221 water molecules). In Figure 2, we show the radial distribution of water molecules around the nonesterified oxygen atoms in the DMP anion computed from force field based simulations of systems with 189, 118, 62, 37, and 23 water molecules (corresponding box edge lengths 19.8, 16.9, 13.7, 11.5, and 9.8 Å, respectively). In the largest system g(r) f 1 at r ≈ 8 Å. The shell structure of the solvent is reflected in the peaks in the radial distribution functions. The sharp first peak with a maximum at 2.6 Å marks a tight first solvation shell containing about 5 water molecules on average. The DMP concentration in our simulations ranges from 2.1 to 17.6 M for our smallest system consisting of 23 water molecules. Even though the concentration is high, it is evident from the near coincidence of the first peak in all of the
Dimethyl phosphate Anion
J. Phys. Chem. B, Vol. 105, No. 24, 2001 5829 TABLE 1: Geometric Parameters for the Dimethyl phosphate Anion from Gas Phase Geometry Optimization and Molecular Dynamics Simulations in Aqueous Solutiona geometric parameter P-O3′ P-O5′ P-OB P-OA O3′-C O5′-C O3′-P-O5′ O3′-P-OB O5′-P-OA OA-P-OB C-O3′-P-O5′ C-O5′-P-O3′
DFT opt ab initio MD force field MD MP2 6-31G** (gas) (solution) (solution) (gas5) 1.68 1.68 1.49 1.49 1.44 1.44 99 105 105 125 75 75
1.63 (0.03) 1.63 (0.03) 1.52 (0.02) 1.52 (0.02) 1.48 (0.04) 1.48 (0.04) 104 (4) 107 (5) 105 (4) 117 (3) -
1.58 (0.03) 1.58 (0.03) 1.48 (0.02) 1.48 (0.02) 1.43 (0.03) 1.43 (0.03) 103 (3) 109 (3) 109 (3) 116 (3) -
1.68 1.68 1.51 1.51 1.43 1.43 99 106 106 126 70 70
a Distances in Å, angles in degrees. MD results are average values, with the standard deviations given in parentheses.
Figure 3. Orientational probability distributions for water molecules within the first solvation shell of the nonesterified oxygen atoms of a dimethyl phosphate anion in a series of force field based MD simulations of different sized systems.
Figure 4. Two views of water density isosurfaces in the first solvation shell of the dimethyl phosphate anion from a force field based MD simulation. The phosphate group is drawn in a ball-and-stick representation, with the nonesterified oxygens (O3′ and O5′) colored dark gray, the esterified oxygens (OA and OB) white, and the phosphorus black. The lighter density contour corresponds to bulk water density, and the darker contour to 1.5 times bulk density.
simulations that the smallest system reliably reproduces the radial distribution of the water molecules in the first solvation shell, as well as part of the second shell in the lowest concentration system. Thus, the structure of the first solvation shell of DMP is insensitive over a wide range of concentration. To quantify the orientation of water molecules in the vicinity of the phosphate group, we define θ as the angle between a PsO* vector in DMP and a OsH vector in water. The probability distribution of this angle (for waters in the first solvation shell only) is plotted in Figure 3. The sharp peak at θ ) 7° arises from tightly bound water molecules within the first solvation shell pointing one of their hydrogen atoms toward a nonesterified oxygen atom in DMP in a nearly linear hydrogen bond. The second peak at θ ) 104° arises from the OH vectors of the water hydrogens not participating in hydrogen bonds with the DMP. The bimodal distributions clearly show that water molecules do not orient themselves such that both hydrogen atoms are simultaneously interacting with an O* atom in DMP (i.e., with the water bisector aligned with the PsO* bond). The similarity of the orientational distribution from the different sized systems demonstrates that the finite size effect on this quantity is also negligible in the smallest system studied. To provide further insight into the water structure in the vicinity of DMP, we show surfaces of constant water density (density isosurfaces38) in Figure 4. It can be seen clearly that there are preferential hydration locations around the DMP anion,
mainly located around the nonesterified oxygens (O*). The rings of higher density around the PsO* bond axes are consistent with the interpretation of the bimodal angular distributions in Figure 3 in terms of water molecules having one hydrogen pointed directly at the O* atom. Compared to the highly structured solvent near the lone oxygen atoms, the solvation shell around the esterified oxygen atoms is not as well-defined due to steric shielding by methyl groups in frequently populated gauche conformations of the CsO3′sPsO5′ and CsO5′sPs O3′ torsions. It is clear from Figure 4 that the phosphate group exerts a strong, anisotropic influence on its primary solvation shell, particularly in the vicinity of the O* atoms. The density distribution in Figure 4 is very similar to that obtained for the phosphate group in MD simulations of phospholipid bilayers.18,39 The role of the anisotropic phosphate solvation in determining the properties of phospholipid bilayer/water interfaces has been discussed at length by Tobias.18 Geometry and Electronic Structure of the DMP Anion in the Gas Phase and Aqueous Solution. Optimized geometrical parameters of the gauche-gauche conformer of the DMP anion obtained using the DFT methods employed in the ab initio MD simulation are compared to the corresponding MP2 6-31G** results,5 as well as average values from the ab initio and force field based MD simulations in solution, in Table 1. As expected, our gas phase DFT results are identical to those obtained previously by Alber et al.15 using the same DFT methods. The close agreement with the MP2 results establishes the validity of the functional and pseudopotentials employed. Comparing the gas-phase geometry with the average structure from the ab initio solution simulation, it is evident that there are small but significant solvent effects on the phosphate geometry. These include a 0.03 to 0.05 Å lengthening of the PsO and OsC bonds, a 5° opening of the O3′sPsO5′ angle, and an 8° closing of the OAsPsOB angle. The CsO5′sPsO3′ dihedral angle remained in the gauche conformation throughout the ab initio solution simulation, but after 5 ps the CsO3′sPs O5′ dihedral underwent a transition from gauche to trans. In addition, the average structure from the ab initio solution simulation is similar to the crystalline ammonium DMP15 structure. The change in structure of DMP on going from solution and the solid is smaller than the change on going from gas to condensed phase (solution or solid). We expect that the presence of the sodium ion had a negligible effect on the DMP anion because the two ions were well-separated throughout the ab initio solution simulation. The closest minimum image approach of the sodium was 4.4 Å to the DMP P atom, 5.0 Å
5830 J. Phys. Chem. B, Vol. 105, No. 24, 2001
Tobias and Kuo
TABLE 2: Bond Ionicity Parameters for the Dimethyl phosphate Anion in Gas Phase and Aqueous Solutiona parameter
gas phase gaucheb
BIPO(P) BICO(C) BICH(C)
0.117 0.101 0.003
solution Gaucheb
solution transc
0.108 (0.006) 0.124 (0.007) -0.008 (0.012)
0.107 (0.006) 0.123 (0.006) -0.007 (0.012)
a Average values from MD simulations (standard deviations in parentheses). See text for the definition of the BI parameters. b BI values computed when the O-C-P-O torsions are in gauche conformations. c BI values computed when the O-C-P-O torsions are in trans conformations.
to OA, 3.4 Å to OB, 3.6 Å to O3′, 5.6 Å to O5′, 4.0 Å to CA, and 5.5 Å to CB. In addition to structural changes due to solvation, an analysis based on the WFCs has revealed that the nature of the chemical bonds of the DMP anion is different in the gas phase and solution. We discuss the solvent effect on bonding in terms of the bond ionicity parameters, BIAB(A), introduced by Alber et al.15
BIAB(A) )
(
)
dA rA dAB rA + rB
where dA is the distance between atom A and the WFC corresponding to the bonding electron pair along the AsB bond, dAB is the bond length, and rA the covalent radius of atom A. Note that a positive (negative) value of BIAB(A) implies that the AsB bond is polarized away from (toward) atom A. Alber et al. have shown that the BI parameters constitute a convenient measure of bond polarity, with larger BI magnitudes indicating more polar bonds.15 The average bond ionicity parameters, BIPO(P), BICO(C), and BICH(C), from the ab initio simulation of the DMP anion in solution are compared to the corresponding values in the optimized gas-phase structure in Table 2. The BICH(C) are small, as expected, reflecting the small polarity of the Cs H bond. It is interesting to note that the direction of the small polarization of the CsH bond changes on moving from the gas phase to solution. The relatively large, positive BI parameters (>0.1) of the PsO and CsO bonds reflect significant polarization of these bonds toward the oxygen atoms. In the gas phase, the PsO bond is more polarized than the CsO bond, whereas in solution the CsO bond is more polarized. The changes in polarization of these bonds follow the change in bond lengths on going from gas phase to solution, as an increase (decrease) in bond length is accompanied by an increase (decrease) of polarization. As noted by Alber et al. in their comparison of the DMP anion in the gas phase in the solid state,15 the bond polarities are quite sensitive to environment, and relatively insensitive to the conformation of the anion. The difference in polarity of the PsO bond between the gas phase and solution is very similar to that between the gas phase and the solid. Comparison of the Solvent Structure in the Force Field Based and ab Initio Simulations. The radial distribution of the water oxygen atoms around the nonesterified phosphate oxygen (O*) atoms from the ab initio simulation is compared to that from the force field based simulation in Figure 5. The first peak corresponding to the first solvation shell is shorter and broader, and shifted outward by about 0.1 Å in the ab initio simulation compared to the force field simulation. The average coordination number of the two O* atoms within the first solvation shell determined by integration of the first peak in the radial distribution function is 5.5 for the force field based simulation and 5.0 for the ab initio simulation. An analogous difference between ab initio and force results is also evident in
Figure 5. Comparison of radial distribution functions for water oxygen atoms around the nonesterified oxygen atoms of a dimethyl phosphate anion computed from ab initio and force field based MD simulations.
Figure 6. Comparison of orientational probability distributions for water molecules within the first solvation shell of the nonesterified oxygen atoms of a dimethyl phosphate anion computed from ab initio and force field based MD simulations.
the orientational distributions of the water molecules in the first solvation shell plotted in Figure 6 there is a slight sharpening of the peak around θ ) 7°, whereas the second peak shows a slight broadening plus a shift by 2° degrees to larger angles in the ab initio simulation. Overall, the data in Figures 5 and 6 show that the force field based simulation predicts a slightly more structured and well-defined first solvation shell around the lone oxygens on DMP compared to the ab initio simulation. Electronic Structure of Water. In the ab initio simulation the electronic structure of the water molecules fluctuates naturally in response to changes of microenvironment. These fluctuations are reflected in changes of the dipole moments of the individual water molecules, which may be computed from the ion and WFC positions, assuming that the electronic charge is concentrated in point charges located on the WFCs.13,14 In Figure 7 ,we have plotted a histogram of all the water dipole moments from the ab initio simulation of sodium DMP in solution. The histogram bears a striking similarity to the one obtained from an ab initio solution of pure water,13,14 and this similarity suggests that, although a large fraction of the water molecules are in close proximity to the DMP or sodium ions in
Dimethyl phosphate Anion
Figure 7. Histogram of water dipole moments from an ab initio MD simulation of sodium dimethyl phosphate in aqueous solution.
Figure 8. Distribution of water dipole moments as a function of the distance between the water oxygen atoms and the nonesterified oxygen atoms in the dimethyl phosphate anion during an ab initio MD simulation.
our simulation, the average electronic structure of the water molecules is not strongly perturbed by the presence of the ions. The average water dipole moment in our simulation, 2.96 D, is very similar to that in pure water, 2.95 D.13,14 To examine more closely if water molecules are polarized by the DMP anion, we have computed the dipole moment histogram as a function of the distance between the water oxygen and the DMP O* atoms (Figure 8). It is clear in Figure 8 that the distribution of dipole moments is essentially invariant to the proximity of the water to the DMP anion. Therefore, we conclude that the presence of DMP in water does not significantly affect the dipole moment of the solvent molecules, even those in the first solvation shell of the O* atoms. Solvation of the Sodium Cation. Although the focus of this paper is on the DMP anion, the solvation of the sodium cation in the ab initio simulation is also of fundamental interest.40 As mentioned above, the sodium cation was initially placed well separated from the DMP anion, and the ab initio simulation was too short to observe appreciable diffusion of either of the ions. Hence, we cannot say anything about the association (or lack thereof) of the DMP anion and the sodium cation. The radial distribution of water oxygen atoms around the sodium cation is plotted in Figure 9. There are two peaks, one
J. Phys. Chem. B, Vol. 105, No. 24, 2001 5831
Figure 9. Radial distribution of water oxygen atoms around the sodium ion in the ab initio simulation.
Figure 10. Distribution of water dipole moments as a function of the distance between the water oxygen atoms and the sodium cation during an ab initio MD simulation.
centered at 2.45 Å, and the other at roughly 4.51 Å, corresponding to the first two solvation shells, respectively. Our g(r) is similar to that from an ab initio MD simulation of an isolated sodium cation in an 11.74 Å cubic box containing 53 water molecules, but there are some quantitative differences. In particular, whereas the first peak occurs at the same distance as that of White et al., ours is slightly higher. In addition, the first minimum is significantly lower in our g(r). We calculate a coordination number of 6.0, which may be compared to the 5.2 obtained by White et al. Overall it appears that the sodium hydration shell is more ordered in our simulation, presumably due to the proximity of the bulky DMP anion. In Figure 10 we show the dipole moment histogram as a function of the distance between the sodium ion and the water oxygen atoms. It appears that the sodium ion, like the DMP anion, does not appreciably influence the electronic structure of the water molecule. 4. Conclusions We have used classical, force field based and ab initio DFT based MD simulations to reveal several novel phenomena associated with the hydration of the dimethyl phosphate anion. Because of the short length of the production ab initio simulation (4.2 ps), sampling of the dihedral angles was not sufficient to
5832 J. Phys. Chem. B, Vol. 105, No. 24, 2001 study the conformational energetics. However, we have been able to characterize the molecular and electronic structure of both the DMP anion and its first solvation shell in detail. The radial and orientational distributions of the water molecules revealed a softening of the structure of the first solvation shell in the ab initio simulation compared to the force field simulation. In addition, the ab initio simulation has shown that both the geometry and electronic structure of the DMP anion change significantly on moving from the gas phase to aqueous solution. Finally, we have found that the electronic structure of the water molecules, as reflected in their dipole moments, is hardly affected by the presence of the sodium and dimethyl phosphate ions. Thus, it appears that empirical force fields with fixed charges that do not explicitly account for electronic polarization should be able to capture the essential features of the solvation of the phosphodiester linkage. Furthermore, with solute molecules parametrized for condensed phase, empirical force field based calculations should give good structural agreement with ab initio simulations, with little dependence on the detailed nature of the condensed phase environment. We believe the results reported in this paper will prove to be useful in guiding future developments of lipid and nucleic acid force fields. Acknowledgment. We are grateful to Mark Tuckerman for assistance in setting up the ab initio MD simulation, and Frank Alber and Paolo Carloni for providing the phosphorus pseudopotential and helpful discussions. We also thank the School of Physical Sciences at UCI for a generous allocation of time on their Origin 2000 computer. References and Notes (1) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. J. Chem. Phys. 1953, 21, 1087. (2) Rao, M.; Pangali, C.; Berne, B. J. Mol. Phys. 1979, 37, 1773. (3) Owicki, J. C.; Scheraga, H. A. Chem. Phys. Lett. 1979, 47, 600. (4) Alagona, G.; Ghio, C.; Kollman, P. J. Am. Chem. Soc. 1985, 107, 2229-2239. (5) Florian, J.; Strajbl, M.; Warshel, A. J. Am. Chem. Soc. 1998, 120, 7959-7966. (6) Jayaram, B.; Mezei, M.; Beveridge, D. L. J. Comput. Chem 1987, 8, 917-942. (7) Jayaram, B.; Mezei, M.; Beveridge, D. L. J. Am. Chem. Soc. 1988, 110, 1691-1694. (8) Jayaram, B.; Ravishaker, G.; Beveridge, D. L. J. Phys. Chem. 1988, 92, 1032-1034. (9) MacKerell, A. D. J. Chim. Phys. 1997, 94, 1436-1447. (10) Liang, C. X.; Ewig, C. S.; Stouch, T. R.; Hagler, A. T. J. Am. Chem. Soc. 1993, 115, 1537-1545.
Tobias and Kuo (11) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471-2474. (12) Parrinello, M. Solid State Comm. 1997, 102, 107-120. (13) Silvestrelli, P. L.; Parrinello, M. Phys. ReV. Lett. 1999, 82, 33083311. (14) Silvestrelli, P. L.; Parrinello, M. J. Chem. Phys. 1999, 111, 35723580. (15) Alber, F.; Folkers, G.; Carloni, P. J. Phys. Chem. B 1999, 103, 6121-6126. (16) Alber, F.; Folkers, G.; Carloni, P. Theochem J. Mol. Struct. 1999, 489, 237-245. (17) Hutter, M.; Clark, T. J. Am. Chem. Soc. 1996, 118, 7574-7577. (18) Tobias, D. J. Water and Membranes: Molecular Details from MD Simulations; Bellissent-Funel, M.-C.; Ed.; IOS Press: Washington, DC, 1999; pp 293-310. (19) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117-1157. (20) Mackerell, A. D.; Wiorkiewicz-Kuczera, J.; Karplus, M. J. Am. Chem. Soc. 1995, 117, 11 946-11 975. (21) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (22) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577-8593. (23) Tuckerman, M. E.; Yarne, D. A.; Samuelson, S. O.; Hughes, A. L.; Martyna, G. J. Comput. Phys. Comm. 2000, 128, 333-376. (24) CPMD, 3.3a ed.; Hutter, J., Alavi, A., Deutsch, T., Bernasconi, M., Goedecker, S., Marx, D., Tuckerman, M., Parrinello, M., Eds.; MPI fur Festkorperforschung and IBM Zurich Research Laboratory: 1995-1999. (25) Hohenberg, P.; Kohn, W. Phys. ReV. 1964, 136, B86. (26) Kohn, W.; Sham, L. J. Phys. ReV. 1965, 140, A1133. (27) Becke, A. D. Phys. ReV. A 1988, 38, 3098-3100. (28) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785-789. (29) Sprik, M.; Hutter, J.; Parrinello, M. J. Chem. Phys. 1996, 105, 1142-1152. (30) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993-2006. (31) Hartwigsen, C.; Goedecker, S.; Hutter, J. Phys. ReV. B 1998, 58, 3641-3662. (32) Martyna, G. J.; Tuckerman, M. E. J. Chem. Phys. 1999, 110, 28102821. (33) Silvestrelli, P. L.; Marzari, N.; Vanderbilt, D.; Parrinello, M. Solid State Comm. 1998, 107, 7-11. (34) Marzari, N.; Vanderbilt, D. Phys. ReV. B 1997, 56, 12 847-12 865. (35) Boys, S. F. Quantum Theory of Atoms, Molecules, and the Solid State; A Tribute to John C. Slater; Academic Press: New York, 1966; p 253. (36) Vanderbilt, D.; King-Smith, R. D. Phys. ReV. B 1993, 48, 44424455. (37) Berghold, G.; Mundy, C. J.; Romero, A. H.; Hutter, J.; Parrinello, M. Phys. ReV. B 2000, 61, 10 040-10 048. (38) Shelley, J. C.; Sprik, M.; Klein, M. L. Langmuir 1993, 9, 916926. (39) Tu, K.; Tobias, D. J.; Blasie, J. K.; Klein, M. L. Biophys. J. 1996, 70, 595-608. (40) White, J. A.; Schwegler, E.; Galli, G.; Gygi, F. J. Chem. Phys. 2000, 113, 4668-4673.
