Adaptive resolution simulation of liquid water

Matej Praprotnik^1,3, Silvina Matysiak², Luigi Delle Site¹, Kurt Kremer¹ and Cecilia Clementi²

¹ Max-Planck-Institut für Polymerforschung, Ackermannweg 10, D-55128 Mainz, Germany
² Department of Chemistry, Rice University, 6100 Main Street, Houston, TX 77005, USA

³ On leave from: National Institute of Chemistry, Hajdrihova 19, SI-1001 Ljubljana, Slovenia
E-mail: dellsite@mpip-mainz.mpg.de, kremer@mpip-mainz.mpg.de and cecilia@rice.edu

Received 8 June 2007, in final form 10 June 2007
Published 5 July 2007

1. Introduction

Scientific problems in physics are usually classified by typical scales of observable quantities such as energy, length or time. The appropriately determined scale together with a tailored model and/or experiment is the starting point of a successful study. The specific properties under investigation then typically define the appropriate level of resolution. Using the optimal set of degrees of freedom (DOFs) guarantees efficiency, accuracy and avoids huge amounts of unnecessary detail, which might even obscure the underlying physics. This approach dates back to the very beginning of modern physics and finds its logical continuation in systematic coarse-graining efforts for modern computational materials science and biophysics problems [1–4], where full-blown all-atom simulations are often beyond the possibilities of current and near-future computers. However, in many problems in materials science and biology different time and length scales are intrinsically interconnected, far beyond giving constant prefactors at the next coarser level of detail. Multiscale simulations are emerging as a promising tool for such problems [5–14]. Ideally, in addition to coarse-graining one would like to be able to step back to characteristic fine-grained configurations as well [15, 16]. Recently, an adaptive multiscale scheme (AdResS) has been proposed by some of us that also allows a free exchange of particles between regions of different resolution and to adjust locally the level of detail at will, while maintaining equilibrium with the fluctuating environment [17].

In this paper we address the above issue for liquid water, which continues to be a very active research field, for the obvious reason that water is the biological solvent of choice. Although many classical models have been developed for the simulation of liquid water [18–21], there is still no unique model that can reproduce all of its anomalous properties. Water, though a small molecule, intrinsically influences its surroundings on multiple scales and plays different roles at different scales. For instance, in the case of biomolecules in solution, at short length scales the hydrogen bonding with water governs the local shape and stability of folded biopolymers, whereas the `hydrophobic effect' drives the organization of the system at longer time and length scales (for instance in the formation of membranes in aqueous solution). The approach presented in this paper represents a first step towards multiscale modelling of such complex scenarios. For fluctuating systems, such as adsorption of water on a hydrophilic surface, we envisage a system in which an atomistic description is needed only near the surface, and a mesoscopic resolution can be used for water molecules further away, with molecules freely moving around and interacting with each other only with the local, more or less detailed level of resolution. Another typical example where our method will find application is that of a solvated macromolecule, e.g. a protein in water. In such a system, the solvent surrounding the macromolecule can be represented with a sufficiently high level of detail so that the specific interactions between the solvent and the solute will be correctly taken into account. On the other hand, the solvent farther away from the macromolecule, where the high resolution is not needed, can be represented on a coarse-grained level [22].

To accomplish this goal, we first propose and test a new single-site water model that reproduces the essential thermodynamic and structural features remarkably well, e.g. the tetrahedral packing of liquid water at standard external conditions, as obtained by detailed all-atom simulations^Note1 . The definition of this new single-site coarse-grained water model is needed for an adaptive resolution simulation, in order to match to the original all-atom model sufficiently well to have a smooth and fast transition between resolutions. Previously proposed simplified water models do not have the required properties to obtain such a smooth transition.

In the second step we adapt the AdResS scheme to define a robust and physically accurate procedure to smoothly join the explicit regions of atomistic and coarse-grained resolution. The resulting multiscale hybrid/mesoscopic model system [17, 23–25] is composed of explicit and coarse-grained molecules as presented in figure 1.

2. Coarse-graining of water

Recent work has focused on simplified coarse-grained models that can reproduce qualitatively the all-atom centre-of-mass radial distribution function (rdf) of water [27–31]. To reproduce further structural properties, such as the orientational preferences for nearest neighbour configurations (that is essential for hydrogen bonding), previous one-site models [28, 29] require a posteriori adjustment of the local environment [29]. The one-site model of water presented in this paper can reproduce the structural and thermodynamic properties of a widely used all-atom water model (namely, Rigid TIP3P [18]) at standard temperature and pressure with remarkable accuracy, without requiring any orientational corrections^Note2 . To construct the model we follow an iterative inverse statistical mechanics approach proposed by Lyubarstev et al  [32] (for alternative procedures see also [2] and [33]). The aim of this scheme is to numerically build a coarse-grained effective Hamiltonian that can mimic the behaviour of a set of physical observables (in this case the centre-of-mass rdf) of the system under consideration (see, e.g., Matysiak et al [34, 35] where this idea has been used to define a coarse-grained protein Hamiltonian `anchored' to experimental data). Additionally, in order to match the pressure of the coarse-grained to the all-atom model, after each iteration, a weak constant force is added to the effective force in such a way that the total effective force and potential energy are zero at the 7 Å cut-off distance [23, 33].

Figure 2 shows the final results obtained when this procedure is applied to a box of water molecules (see section 4 for details). Perfect agreement between the all-atom and coarse-grained centre-of-mass rdfs is reached with our optimized effective potential (shown in the inset of the figure) after eight iterations. The tabulated values of the potential used in the simulations are provided in the supplementary material (available at ). The effective potential has a first primary minimum at about 2.8 Å corresponding to the first peak in the centre-of-mass rdf. The slightly weaker and significantly broader minimum at 4.5 Å corresponds to the second hydration shell. The combined effect of the two minima leads to a local packing close to that of the all-atom TIP3P water. The average number of nearest neighbours (as obtained by integrating the rdf over the first hydration shell) is 3.92 ± 0.28 in the explicit region, and 4.14 ± 0.33 in the coarse-grained region, in both case signalling a tetrahedral packing^Note3 . Our effective coarse-grained potential is qualitatively different from the previously suggested potentials [27–29]: while in previous one-site models the deepest minimum corresponds to the second hydration shell, the absolute minimum in our model is found in the first shell. In order to more thoroughly quantify the structural properties of our model, which are not completely defined by the rdf, we computed the angular distribution between the centre-of-masses of three nearest neighbour molecules and the distribution of the orientational order parameter q as defined by Errington et al  [37]: , where ψ_jk is the angle formed by the lines joining the oxygen atom of a given molecule (the centres of the spheres in the one-site coarse-grained model, respectively) and those of its nearest neighbours j and k. The parameter q measures the extent to which a molecule and its four nearest neighbours adopt a tetrahedral arrangement [37]. The good agreement between the explicit and coarse-grained water models shown in figures 2(b) and (c) indicates that although our coarse-grained model is spherically symmetric and does not have any explicit directionality, it approximately captures the orientational preferences of hydrogen bonding.

It is worth noting that the coarse-graining significantly speeds up the simulation, with a total gain in computational time of a factor ~17–20 when compared to atomistic simulations. Two factors contribute to this speed-up:

(i) 

The reduction on the number of interactions computed at each time step and the simplified form of interactions associated with the coarse-graining accounts for a factor ~8–10.

(ii) 

In addition, there is an intrinsic timescale difference in the diffusive dynamics of the coarse-grained system, which is faster by about a factor of 2 than the predictions of all-atom simulations because of the reduced number of DOFs (the diffusion coefficient of the centre-of-mass for the Rigid TIP3P model is D≈(3.7 ± 0.3) × 10^–9 m² s^–1 while for the coarse-grained model it is D≈(8.1 ± 0.4) × 10^–9 m² s^–1).

3. Multiscale water

In order to smoothly couple the atomistic and coarse-grained resolutions, we adapt the multiscale AdResS (adaptive resolution scheme) procedure of Praprotnik et al  [17, 23] outlined in section 4 to directly link a model system composed of explicit and coarse-grained molecules. Detailed comparisons between the bulk explicit simulations and the explicit regime in our hybrid set-up prove that our approach does not alter the structural properties of the water model studied.

When explicit water molecules approach a neutral fluid surface (as happens upon entering the coarse-grained liquid water region) the hydrogen bond network formed among the explicit molecules can be strongly perturbed. As a result, the explicit molecules near the resolution interface could orient and behave differently from the bulk. However, this spurious effect can be avoided by using an appropriate interface layer between the explicit and coarse-grained regime. It has been shown that the perturbation introduced by the interaction of explicit water with a neutral surface extends around 10 Å into the water layer [39]. This distance is comparable to the maximum range of electrostatic interactions in our simulations. This distance range determines the width of the interface layer to allow the molecules to equilibrate on-the-fly while moving from one resolution to the other. The results presented in this paper have been obtained by using a conservative choice for the width of the interface layer, that is, double the maximum range of interactions^Note4 .

The results presented in figure 3 show that indeed the interface layer completely removes any orientational bias that could have been introduced in the explicit water by the interactions with the coarse-grained regime. The average orientations of three vectors are shown in the figure, namely the dipole moment of a water molecule, the vector perpendicular to the plane of the molecule and the vector joining the two hydrogen atoms of the molecule, as a function of the coordinate x spanning across the two resolutions in the simulation box. We denote the angles formed by these vectors and the normal vector pointing towards our coarse-grained water liquid as α, β and γ, respectively. The average orientation of the vectors can be quantified by the average cosine value of these angles^Note5 . The results in figures 2(d), (e), and 3 indicate that the transition regime neutralizes the interface effect of the coarse-grained liquid water, such that the molecules have fully equilibrated orientational DOFs when they enter the explicit regime.

The structural (see figure 2) and thermodynamic properties of the bulk all-atom simulations are correctly reproduced by the hybrid system with the same mean temperature (0.1% difference) and pressure (0.5% difference). The normalized density for the hybrid system is homogeneous in the coarse-grained and explicit regions with (very) small oscillations in the transition regime (see figure 4). In order to prove the free exchange of molecules between the different regimes we have computed the time evolution of a diffusion profile for molecules that were initially localized at the interface layer. Figure 4 shows that these molecules spread out asymmetrically with time. This asymmetry arises from the aforementioned difference in diffusion coefficient between the all-atom and coarse-grained regions. It is possible to obtain the same diffusional dynamics across the different resolutions by adding a local Langevin thermostat, as mentioned above (unpublished results). However, it is worth stressing that this difference in timescale can be advantageous for reaching longer simulation times in systems where multiple length- and time scales are intrinsically present.

4. Model and methods

In the present study, half of the simulation box is occupied by atomistic (Rigid TIP3P) water molecules while the other half is filled with the same number of corresponding coarse-grained molecules as schematically presented in figure 1. The two regions freely exchange molecules through a transition regime, where the molecules change their resolution and their number of DOFs accordingly. The interface region contains hybrid molecules that are composed of an all-atom molecule with an additional massless centre-of-mass particle serving as an interaction site. The transition, which needs to be smooth in order to be used in molecular dynamics simulations, is governed by a weighting function that interpolates the interaction forces between the two regimes, and assigns the identity of the particle. We used the weighting function defined in [17] in such a way that w  =  1 corresponds to the atomistic region, and w  =  0 to the coarse-grained region, whereas the values 0<w<1 correspond to the interface layer. The atomic and mesoscopic length scales are coupled via the intermolecular force acting between centres-of-mass of molecule α and β as:

where is the sum of all pair intermolecular atom interactions between explicit atoms of the molecules α and β and is the corresponding effective intermolecular force between their centres-of-mass. r_iαjβ  =  r_iα–r_jβ is the vector between atom i in molecule α and atom j in molecule β and R_αβ  =  R_α–R_β the vector between the centres-of-mass of molecules α and β, with the corresponding X coordinates X_α and X_β. Each time a molecule crosses a boundary between the different regimes it gains or loses (depending on whether it leaves or enters the coarse-grained region) its equilibrated rotational DOFs while retaining its linear momentum. To supply or remove the latent heat caused by the switch of resolution this method is employed together with a Langevin thermostat [17]. In order to remove an increased pressure in the transition regime, an interface pressure correction, as tested in [23], is employed. As discussed in an earlier publication, it is important to interpolate the forces and not the interaction potential if Newton's third law is to be satisfied [24]. The key methodological issue in applying our approach to water is how to treat the long-range electrostatic interactions. We have chosen the reaction field (RF) method, in which all molecules outside a spherical cavity of a molecular-based cut-off radius R_c  =  9 Å are treated as a dielectric continuum with a dielectric constant ε_RF [42–45]. The Coulomb force acting on a partial charge e_iα, belonging to the explicit or hybrid molecule α, at the centre of the cut-off sphere, due to a partial charge e_jβ, belonging to the explicit or hybrid molecule β, within the cavity is:

There are three main reasons that make the RF method a natural choice to be applied with our scheme (equation (1)): the first reason is that it does not impose any artificial periodicity (that would interfere with the definition of the interface layer of the system), the second one is the pairwise form of the reaction field term (equation (2)), and the third one is that, as in our scheme, the RF method works accurately only in a combination with a thermostat [45]. The last point can be better understood if the switching of the resolution is seen in analogy to a geometry-induced phase transition [24].

5. Conclusion

We have presented a general and computationally efficient multiscale model for water. We have shown that a new one-site coarse-grained model for water molecules can reproduce thermodynamic and structural properties of the Rigid TIP3P water model remarkably well, and that a smooth transition and free exchange between coarse-grained and all-atom resolutions is possible. The presented multiscale approach can be applied to any other either flexible or rigid non-polarizable classical water model, e.g. SPC or SPC/E [40–42]. We envisage that such a multiscale resolution of water will play an important role in the modelling of wet/dry interfaces and in biomolecular simulations. The gain in computational speed will allow for larger systems and more systematic studies and should at the same time extend the possibilities of such simulations to a larger user group.

Acknowledgments

We wish to thank the NSF-funded Institute for Pure and Applied Mathematics at UCLA where this work was first planned. This work has been supported in part by grants from NSF(CAREER CHE-0349303, CNS-045433 and CCF-0523908), Texas-ATP (003604-0010-2003), the Robert A Welch Foundation (Norman-Hackerman Young Investigator award and C-1570) (CC) and the Volkswagen Foundation (KK & LDS). The Rice University Cray XD1 Cluster ADA used for the calculations is supported by the NSF under grant CNS-0421109, and a partnership between Rice University, AMD and Cray.

References

[1]

Nielsen S O, Lopez C F, Srinivas G and Klein M L 2004 J. Phys.: Condens. Matter 16 R481 
IOPscience

[2]

Izvekov S and Voth G A 2006 J. Chem. Phys. 125 151101 
CrossRefPubMed

[3]

Harmandaris V A, Adhikari N P, Van der Vegt N F A and Kremer K 2006 Macromolecules 39 6708 
CrossRef

[4]

Das P, Matysiak S and Clementi C 2005 Proc. Natl Acad. Sci. USA 102 10141 
CrossRefPubMed

[5]

Delle Site L, Abrams C F, Alavi A and Kremer K 2002 Phys. Rev. Lett. 89 156103 
CrossRefPubMed

[6]

Neri M, Anselmi C, Cascella M, Maritan A and Carloni P 2005 Phys. Rev. Lett. 95 218102 
CrossRefPubMed

[7]

Koumoutsakos P 2005 Annu. Rev. Fluid Mech. 37 457 
CrossRef

[8]

Backer J A, Lowe C P, Hoefsloot H C J and Ledema P D 2005 J. Chem. Phys. 123 114905 
CrossRefPubMed

[9]

Fabritiis G D, Delgado-Buscalioni R and Coveney P V 2006 Phys. Rev. Lett. 97 134501 
CrossRefPubMed

[10]

Lyman E, Ytreberg F M and Zuckerman D M 2006 Phys. Rev. Lett. 96 028105 
CrossRefPubMed

[11]

Shi Q, Izvekov S and Voth G A 2006 J. Phys. Chem. B 110 15045 
CrossRef

[12]

Christen M and van Gunsteren W F 2006 J. Chem. Phys. 124 154106 
CrossRefPubMed

[13]

Weinan E, Engquist B, Li X, Ren W and Vanden-Eijnden E 2007 Commun. Comput. Phys. 2 367 

[14]

Shehu A, Clementi C and Kavraki L 2006 Proteins Struct. Funct. Bioinf. 65 164 
CrossRef

[15]

Heath A, Kavraki L and Clementi C 2007 Proteins Struct. Funct. Bioinf. 68 

[16]

Hess B, Leon S, Van der Vegt N and Kremer K 2006 Soft Matter 2 409 
CrossRef

[17]

Praprotnik M, Delle Site L and Kremer K 2005 J. Chem. Phys. 123 224106 
CrossRefPubMed

[18]

Jorgensen W L, Chandrasekhar J, Madura J D, Impey R W and Klein M L 1983 J. Chem. Phys. 79 926 
CrossRef

[19]

Mahoney M W and Jorgensen W L 2000 J. Chem. Phys. 112 8910 
CrossRef

[20]

Praprotnik M and Janežič D 2005 J. Chem. Phys. 122 174103 
CrossRefPubMed

[21]

Guillot B 2002 J. Mol. Liq. 101 219 
CrossRef

[22]

Praprotnik M, Delle Site L and Kremer K 2007 J. Chem. Phys. 126 134902 
CrossRefPubMed

[23]

Praprotnik M, Delle Site L and Kremer K 2006 Phys. Rev. E 73 066701 
CrossRef

[24]

Praprotnik M, Kremer K and Delle Site L 2007 Phys. Rev. E 75 017701 
CrossRef

[25]

Praprotnik M, Kremer K and Delle Site L 2007 J. Phys. A: Math. Theor. 40 F281 
IOPscience

[26]

Limbach H J, Arnold A, Mann B A and Holm C 2006 Comput. Phys. Commun. 174 704 
CrossRef

[27]

Garde S and Ashbaugh H S 2001 J. Chem. Phys. 115 977 
CrossRef

[28]

Voth G A and Izvekov S 2005 J. Chem. Phys. 123 134105 
CrossRefPubMed

[29]

Head-Gordon T and Stillinger F H 1993 J. Chem. Phys. 98 3313 
CrossRef

[30]

Soper A K 1996 Chem. Phys. 202 295 
CrossRef

[31]

Nezbeda I 2005 Mol. Phys. 103 59 
CrossRef

[32]

Lyubartsev A P and Laaksonen A 1995 Phys. Rev. E 52 3730 
CrossRef

[33]

Reith D, Putz M and Müller-Plathe F 2003 J. Comput. Chem. 24 1624 
CrossRefPubMed

[34]

Matysiak S and Clementi C 2004 J. Mol. Biol. 343 235 
CrossRefPubMed

[35]

Matysiak S and Clementi C 2006 J. Mol. Biol. 363 297 
CrossRefPubMed

[36]

Lazaridis T and Karplus M 1996 J. Chem. Phys. 105 4294 
CrossRef

[37]

Errington J R and Debenedetti P G 2001 Nature 409 318 
CrossRefPubMed

[38]

Kremer K and Grest G S 1990 J. Chem. Phys. 92 5057 
CrossRef

[39]

Jedlovszky P, Vincze A and Horvai G 2004 Phys. Chem. Chem. Phys. 6 1874 
CrossRef

[40]

Berendsen H J C, Postma J P M, van Gunsteren W F and Hermans J 1981 Intermolecular Forces ed B Pullman (Dordrecht: Reidel) p 331 

[41]

Berendsen H J C, Grigera J R and Straatsma T P 1987 J. Phys. Chem. 91 6269 
CrossRef

[42]

Praprotnik M, Janežič D and Mavri J 2004 J. Phys. Chem. A 108 11056 
CrossRef

[43]

Neumann M 1983 Mol. Phys. 50 841 
CrossRef

[44]

Neumann M 1985 J. Chem. Phys. 82 5663 
CrossRef

[45]

Tironi I G, Sperb R, Smith P E and van Gunsteren W F 1995 J. Chem. Phys. 102 5451 
CrossRef

Notes

Note1  The coarse-grained model, as it is presented here, has been designed to perform well specifically for the selected temperature and pressure (standard external conditions). However, one can adjust the model for different thermodynamic conditions by a re-parametrization of the effective potential between molecules.

Note2  Since a water molecule is electrically neutral the interaction site has a zero electric charge. Furthermore, in contrast to a Stockmayer fluid, the presented one-site model does not have a dipole (or any higher electric multipole). All the electrostatic interactions, including the reaction field correction, are therefore implicitly taken into account via the effective intermolecular potential.

Note3  Our single-site model is therefore significantly different from a Lennard-Jones liquid, which is closely packed with 12 nearest neighbours (see for instance [36]). Indeed, results from extensive tests show that if the functional form of a Lennard-Jones potential is used in the definition of the water coarse-grained model, the iterative inverse statistical mechanics procedure does not converge, and the rdf of the explicit model cannot be reproduced (results not shown).

Note4  However, a smaller interface layer can be used without affecting the fluxes between the layers or the results in the two resolution regimes (results not shown).

Note5  The random orientation corresponds to for the molecular normal vector, and for the vector joining the two hydrogens, while the random orientation of the dipole moment corresponds to . The different values for the random orientations of the vectors arise from the fact that opposite directions are physically indistinguishable both for the molecular normal vector and the vector joining the two hydrogens, while the vector of the dipole moment has a specific directionality [39].