Modeling Diffusive Dynamics in Adaptive Resolution Simulation of Liquid Water

Modeling Diffusive Dynamics in Adaptive Resolution Simulation of Liquid Water

Abstract

We present a dual-resolution molecular dynamics (MD) simulation of liquid water employing a recently introduced Adaptive Resolution Scheme (AdResS). The spatially adaptive molecular resolution procedure allows for changing from a coarse-grained to an all-atom representation and vice-versa on-the-fly. In order to find the most appropriate coarse-grained water model to be employed with AdResS we first study the accuracy of different coarse-grained water models in reproducing the structural properties of the all-atom system. Typically, coarse-grained molecular models have a higher diffusion constant than the corresponding all-atom models due to the reduction in degrees of freedom (DOFs) upon coarse-graining that eliminates the fluctuating forces associated with those integrated-out molecular DOFs. Here, we introduce the methodology to obtain the same diffusional dynamics across different resolutions. We show that this approach leads to the correct description of essential thermodynamic, structural and dynamical properties of liquid water at ambient conditions.

I Introduction

Water is a typical example where the interplay between different length scales determines the relevant properties of the system. Including explicit water in large biomolecular simulations is crucial Pal et al. (2002); Cheung et al. (2002) but normally not feasible due to the large number of water molecules needed to describe biomolecular function. The computational requirements of all-atom simulations in explicit water usually do not allow to investigate biologically relevant time-scalesPraprotnik and Janežič (2005). Most computational approaches renormalize the role of water into the definition of effective inter-residue interactions Cheung et al. (2002) or as a continuous field, in order to focus the computational efforts to simulate the biomolecule. However, the discrete nature of water molecules can play a role in, e.g., controlling protein functionality or the case of interactions between hydrocolloids and water, where hydrogen bonding is important, etc. To overcome this problem, a multiscale modeling of water can be very advantageous to speed up the simulation without losing physics-chemical accuracy.

Multiscale modeling is emerging as a viable way to bridge the various length and time scales involved in complex molecular systemsBroughton et al. (1999); Csanyi et al. (2004); Kremer (2004); Koumoutsakos (2005); Neri et al. (2005); Backer et al. (2005); Abrams (2005); Fabritiis et al. (2006); Lyman et al. (2006); Praprotnik et al. (2007a, b). In general, to build a multiscale model two main issues need to be addressed. The first issue consists in mapping the system into a robust reduced model while preserving the here relevant physico-chemical properties, i.e., radial distribution functions, pressure, and temperature, of the reference all-atom system. The second issue involves the definition of a robust and physically accurate procedure to smoothly join the different resolutions. We have recently proposed the Adaptive Resolution Scheme (AdResS)Praprotnik et al. (2005, 2006) to address both issues for a system of water molecules Praprotnik et al. (2007c). Typically, the reduction in the number of the system degrees of freedom (DOFs) upon coarse-graining introduces a time scale difference between the coarse-grained and explicit region. Although in certain instances the difference in time scale may be advantageous for reaching longer simulation times, in other cases having the correct diffusional dynamics is crucial. This is the case, for instance, if a dynamical property is focus of an investigation as in the translocation of biopolymers through membrane nanopores Matysiak et al. (2006).

In this paper, in order to justify our choice of the recently introduced single-site water modelPraprotnik et al. (2007c) as the optimal coarse-grained water model to be used with AdResS (when using the standard three-site TIP3P water modelJorgensen et al. (1983) as the all-atom water model), we present firstly a detailed analysis and comparison of newly developed coarse-grained water models having a different number of DOFs (from to ) with the single-site water modelPraprotnik et al. (2007c) as well as the TIP3P water modelJorgensen et al. (1983). In this way we can study the contribution of separate DOFs of the all-atom water model to the definition of its structural properties. As it turns out the single-site water model can reproduce the structural and thermodynamical properties of the all-atom system that are relevant for the multiresolution simulation equally well as the more sophisticated models. Therefore, we have chosen it as the coarse-grained model in our adaptive resolution simulations. Next, we introduce the methodology to obtain the same diffusional dynamics across the different resolutions in the hybrid atomistic/coarse-grained (ex-cg) model system composed of explicit and coarse-grained molecules as presented in Figure 1.

Figure 1: (Color online) Top figure: The coarse-grained molecule is represented at the left, and the all-atom water molecule is represented at the right. The middle hybrid molecule interpolates between the two. Bottom figure: The dashed line represents the weighting function w(x) defined in Ref. Praprotnik et al. (2005) and discussed in the text. The values w=1 and w=0 correspond to the atomistic and coarse-grained regions of the hybrid atomistic/mesoscopic system. The full line represents the interface correction weighting function s(w) defined in Figure 6. The value s=1 corresponds to the atomistic and coarse-grained regions, while s=0 when w=1/2. The explicit and interface regions are divided in three regions(A,B,C) of width 4 Å  for computing the cosine distribution in Figure 8.

Ii Mapping of a coarse-grained model to the all-atom representation

In order to map a coarse-grained to an all-atom model, we have to construct an effective potential between coarse-grained molecules in a way that the atomistic structural properties are reproduced. We numerically build such effective Hamiltonian following an ”inverse statistical-mechanics” approach Lyubartsev and Laaksonen (1995); Matysiak and Clementi (2004, 2006). In order not to bias a priori the potential function with the choice of a particular functional form, we introduce a grid approximation, expressing a general pairwise Hamiltonian in the form:

(1)

where is the number of coarse-grained particles in the interval and is the discretized value of the potential in the interval . The averages , which can be interpreted as a discretized radial distribution function (rdf), are functions of the set of parameters . Values of the potential parameters that reproduce the rdf of the all-atom model can be obtained iteratively as follows:

  • The center-of-mass rdf is computed from all-atom simulations and used as a target function. The potential of mean forces is assumed as an initial approximation to the effective potential function:

    (2)

    where is the potential at a distance .

  • By comparing the atomistic center-of-mass rdf with the rdf obtained for the coarse-grained model at the n-th iteration (with corresponding potential parameters ) the -th correction to the effective potential can be found. Namely, the correction to the potential parameters is given by the solution of the system of linear equations Lyubartsev and Laaksonen (1995):

    (3)

    where:

    (4)
  • Ensemble averages and , are computed from MD simulations (with potential parameters ), and the correction to the potential is obtained by solving the system of linear equations (3):

    (5)
  • The corrected potential parameters are then used to perform molecular dynamics simulation with the coarse-grained model and calculate new ensemble averages for the quantities and ,.

The points above are repeated to determine a new set of corrections . The procedure is repeated until convergence is reached within the statistical simulation error.
Additionally, 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 fixed cutoff distance Reith et al. (2003); Praprotnik et al. (2006).

(6)

We use a value of Å  as the rdf for . Depending on the pressure in the current iteration being above or below the target value corresponding to the pressure of the reference all-atom system, can assume a positive or negative value.
For the atomistic representation of water we selected the rigid TIP3P model Jorgensen et al. (1983) that gives a reasonable description of the behavior of liquid water around standard (physiological) temperature and pressure, and is not too expensive computationally. The potential function of rigid TIP3P involves a rigid water with three interaction sites that coincide with the atomic positions. The model uses atom-centered point charges to represent the electrostatic, i.e., positive charges on the hydrogens and a negative charge on oxygen (). The van der Waals interaction between two water molecules is modeled by the Lennard-Jones potential with just a single interaction site per molecule centred on the oxygen atom; the van der Waals interactions involving the hydrogens are omitted. This gives the following potential:

(7)

The values of the parameters and are assigned to reproduce reasonable structural and energetics properties of liquid water Jorgensen et al. (1983).

Iii Coarse-grained models for water

Many simplified coarse-grained models have been developed to reproduce qualitatively the structural properties of waterDill et al. (2005). Existing coarse-grained models can of course only reproduce thermodynamical properties of all-atom water to a certain extendGarde and Ashbaugh (2001); Voth and Izvekov (2005); Head-Gordon and Stillinger (1993); Soper (1996); Nezbeda (2005). By using the procedure described above, we have designed water models at different level of coarse-graining (three-site, two-site, and one-site) schematically depicted in Figure 2. The gradually increasing level of coarse-graining allows us then to study the contribution of different DOFs to the structural and thermodynamical properties of the reference all-atom system. In this way we can determine the most appropriate coarse-grained model to be used in our adaptive resolution simulations.

Figure 2: (Color online) Cartoon of different water models employed in this study. From top to bottom: The rigid TIP3P water model, and three-site, two-site, one-site coarse-grained water models, respectively.

Three-site interaction model

A first level of coarse-graining can be introduced by preserving the atomic positions of each atom in the TIP3P water model and replacing the explicit electrostatic interactions (the first term in Eq. (7)) with effective short-ranged interactions. This model preserves explicitly the H-bond directionality but it is computationally less expensive than TIP3P because it removes the long-ranged Coulombic term. The calculated all-atom site-site rdfs are used as input to build the OO, OH and HH effective potentials (see Eqns. 4). In order to quantify the agreement of the rdfs (corresponding to the all-atom and coarse-grained model) we introduce the penalty function Praprotnik et al. (2005) defined as

(8)

where and are the reference site-site rdf of the atomistic and coarse-grained system, respectively, and is the Lennard-Jones constant of the water modelJorgensen et al. (1983). The extremely low values of the obtained ( of the OO rdf is , for OH is and HH is ) for the optimized effective potential for the three-site model indicates a perfect matching of the rdfs. The optimized effective potentials are as expected very similar to the effective potentials of other previously introduced 3-site coarse-grained water modelsSoper (1996); Lyubartsev and Laaksonen (2000). In addition, to check the angular properties (that are not completely defined by the rdf) we computed the distribution of the angle formed between the center-of-mass of three nearest neighbor molecules (Figure 3, top), and the distribution of the orientational order parameter (see Figure 3), as defined by Errington et al. Errigton and Debenedetti (2001):

(9)

where is the angle formed by the lines joining the oxygen atom of a given molecule and those of its nearest neighbors and . The parameter measures the extent to which a molecule and its four nearest neighbors adopt a tetrahedral arrangement. The excellent agreement between the rdfs and the angular distribution suggests that the explicit H-bond directionality is responsible for the local structure. On the contrary, the small deviation of the orientational order parameter distribution suggests that the electrostatic interactions are at least partly responsible for the overall tetrahedral arrangement of the system.

Two-site interaction model

A two-site simplified model of water changes the symmetry of the molecule by using two interacting sites instead of three, one corresponding to the oxygen atom and the second one providing the directionality of the dipole moment. Since the dipole moment of TIP3P is and the charge separation is q=0.82e Jorgensen et al. (1983), the effective separation of the negative and positive charge centers is Å. Thus, we define the two-site model by representing each water molecule by two spheres, which are joined by a rigid rod of length 0.5966 Å. Two-site models have been previously proposed for water Voth and Izvekov (2005), however existing models cannot exactly reproduce the radial structure of the reference all-atom water model. The two-site water model with parameters optimized by the procedure described above can reproduce remarkably well the rdfs for the different pairs of interactions: oxygen-oxygen with , oxygen-dipole with and dipole-dipole with . Figure 3 shows the angular distribution and the orientational order parameter distribution obtained for the two-site model, in comparison with the results of the all-atom and three-site models. The excellent agreement between the rdfs suggests that the inclusion of a dipole center is sufficient for the correct radial structure. The deviation on the angular distribution from the three-site model indicates that the explicit H-bond directionality has a contribution to the angular structure (in addition to the electrostatic contribution, as found in the three-site model).

One-site interaction model

In the one-site coarse-grained model, the water is represented by one spherically symmetrical site having a mass . As demonstrated in Ref. Praprotnik et al. (2007c), the one-site model can reproduce the center-of-mass rdf () and thermodynamic properties of the reference all-atom water model with remarkable accuracy. The optimized effective potential (see inset of Figure 5) has a first primary minimum at about 2.8 Å  corresponding to the first peak in the center-of-mass rdf. A second, slightly weaker and significantly broader minimum at 4.5 Å  corresponds to the second hydration shell. The combined effect of the two leads to a local packing close to that of the all-atom TIP3P water. Our effective coarse-grained potential is quite different from the previously suggested potentialsGarde and Ashbaugh (2001); Voth and Izvekov (2005); Head-Gordon and Stillinger (1993); Johnson et al. (2007): 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.

The good structural agreement between the explicit and coarse-grained models indicates that although our coarse-grained model of water is spherically symmetric and therefore does not have any explicit directionality, it approximately captures the correct local structure. Furthermore, the angular and orientational order parameter distribution coincides with what obtained for the two-site model.

Figure 3: Top figure: The center-of-mass angular distribution between three nearest neighbors for all-atom, 1-site coarse-grained, 2-site coarse-grained and 3-site coarse-grained models for water. Bottom figure: The analogous distributions of the orientational order parameter . The center-of-mass rdf of the different water models coincides exactly with the all-atom rdf (not shown).

Since the one-site model can qualitatively predict the correct angular distribution and is computationally the least expensive of all the coarse-grained models, it was our choice for the multiscale approach (presented in Praprotnik et al. (2007c)).
As a consequence of the reduced number of DOFs, there is a time scale difference in the dynamics of the coarse-grained system with respect to the atomistic one. The diffusion coefficient increases as the level of coarse-graining is increased, with the diffusion coefficient of the center of mass for the TIP3P model , for the three-site model , for the two-site model and for the one-site model .

Using the one-site coarse-grained model one can significantly speed up the simulation, with a total gain in computational time of a factor when compared to atomistic simulationsPraprotnik et al. (2007c). This is due to the reduction of the number of interactions, which are also softer than in the atomistic case, and due to an intrinsic time scale difference in the diffusive dynamics of the coarse-grained system, that is faster of about a factor 2 than the all-atom simulations counterpart (see the above values for the diffusion constant).

Iv Matching the diffusive dynamics of the coarse-grained to the all-atom model

As discussed above, the dynamics of the one-site coarse-grained model is faster than that obtained from all-atom simulations. The speed-up occurs because of the reduction in DOFs upon coarse-graining, which eliminates the fluctuating forces associated with those missing molecular DOFs, leading to the much smoother overall energy landscapeTschöp et al. (1998); Izvekov and Voth (2006). In our MD simulations we use a Langevin thermostat, and the equations of motion are in the form:

(10)

where is the deterministic force and is a stochastic variable with . The coefficient determines the strength of the coupling to the bath [not to be confused with the friction coefficient of the system (for TIP3P )]. When is large compared to the typical time scales in the system, the stochastic dynamics reproduces dynamics obtained by molecular dynamics (MD)Kremer and Grest (1990), while if is small, then the stochastic dynamics deviates significantly from MD. Simulations performed using the the one-site model with same coupling to the bath () that is usually used for TIP3P Tironi et al. (1995) produce an increase of a factor of about in the diffusion time scale of the coarse-grained system. An accelerated dynamics can be advantageous in some cases but as mentioned previously it can be a problem if the dynamical properties themselves are under investigation. The coarse-grained dynamics can be slowed down by increasing the effective friction in the coarse-grained system. That is, by changing , the one-site model can mimic the diffusive dynamic of the all-atom system. Figure 4 shows that the diffusion coefficient D monotonically decreases as the coefficient increases. The coarse-grained model yields the same diffusion coefficient as the all-atom system (for which ) when . The additional frictional noise does not affect the short time dynamics of the system since Kremer and Grest (1990). It is worth mentioning that the structural/thermodynamics of the coarse-grained system are not effected by changes in in this range. Figure 4 shows the perfect agreement between the center-of-mass rdf for the coarse-grained and all-atom system when different values of are used for the two models. However, we have to bear in mind that by employing the Langevin thermostat, Eq. (10), we screen the hydrodynamics. Furthermore, for too large the statics is also altered because one arrives eventually at the Brownian dynamics. In order to correctly reproduce the hydrodynamics one has to resort to the Dissipative Particle Dynamics (DPD) thermostatSoddemann et al. (2003) as was done for instance in Ref.Praprotnik et al. (2007d) for the case of a macromolecule in the hybrid atomistic/mesoscale solvent.

Figure 4: The diffusion coefficient as a function of the friction constant. The triangle indicates the value obtained for the all-atom diffusion coefficient when a value is used in the simulation. The black dots indicate the diffusion coefficient obtained for the coarse-grained system for different values of . The two circles indicate the diffusion coefficients obtained for the coarse-grained system when and , respectively, are used in the simulation. The inset compares the rdfs obtained for the all-atom system and the coarse-grained system for different values of the friction coefficients.

V Adaptive resolution scheme (AdResS)

As we have proposed in Praprotnik et al. (2007c), we design an adaptive multiscale system where half of the simulation box is occupied by atomistic water molecules and the other half by the corresponding one-site coarse-grained molecules, respectively, as schematically presented in Figure 1. In order to smoothly couple the regimes of high and low level of detail, we apply the AdResS scheme Praprotnik et al. (2005, 2006); Praprotnik et al. (2007c), that allows the molecules to move freely between regimes without feeling any barrier in between. The interface region contains hybrid molecules that are composed of an all-atom molecule with an additional massless center-of-mass particle serving as an interaction site. The transition is governed by a weighting function that interpolates the molecular interaction forces between the two regimes, and assigns the identity of the molecules. In the present work we used the weighting function defined in Praprotnik et al. (2005) and described in Figure 1. The function is defined in such a way that corresponds to the atomistic region, and to the coarse-grained region, whereas the values correspond to the interface layer. This interpolation leads to intermolecular forces acting on the center-of-mass of the molecules and as:

(11)

is the total intermolecular force acting between center of mass of molecules and . is the sum of all pair atomic interactions between explicit water atoms of molecule and explicit water atoms of molecule , and is the effective pair force between the center-of-mass of two water molecules. and are the center-of-mass coordinates of molecules and . Note that the AdResS as given by Eq. (11) satisfies Newton’s Third Law. This is crucial for the diffusion of molecules across the resolution boundariesPraprotnik et al. (2007b). Since the total force as defined by Eq. (11) depends on the absolute positions of the particles and not only on their relative distances it is in general not conservative, i.e., the work done by it depends on the path taken in the transition regime. Hence the corresponding potential does not exist and the total potential energy of the system can not be definedGoldstein (1980); Ens (). Still, the intermolecular forces between molecules outside the transition regime are conservative with well defined potentials, i.e., the all-atom or effective coarse-grained potentialsPraprotnik et al. (2005).

Each time a molecule crosses the boundary between the different regimes it gains or looses (depending on whether it leaves or enters the coarse-grained region) its equilibrated rotational DOFs while retaining its linear momentum. By this choice of interactions in the interface region, the hybrid molecule interacts with molecules in the coarse-grained region on a coarse-grained level. On the other hand, the interactions of the hybrid molecules with the molecules in the explicit region are a combination of the explicit-explicit and cg-cg interactions to smoothly and efficiently equilibrate additional DOFs upon moving in the explicit regime Praprotnik et al. (2005, 2006). This change of resolution, which can be thought in terms of similarity to a phase transitionPraprotnik et al. (2007a), requires to supply or remove latent heat and thus must be employed together with i.e. a Langevin thermostat Praprotnik et al. (2005). The thermostat is coupled locally to the particle motion and provides a mean to deliver or absorb the latent heat.

As in Ref. Praprotnik et al. (2007c) the reaction field (RF) method is used, in which all molecules outside a spherical cavity of a molecular based cutoff radius Å  are treated as a dielectric continuum with a dielectric constant Neumann (1983, 1985, 1986); Tironi et al. (1995); Im et al. (2001); Praprotnik et al. (2004). The Coulomb force acting on a partial charge , belonging to the explicit or hybrid molecule , at the center of the cutoff sphere, due to a partial charge , belonging to the explicit or hybrid molecule , within the cavity is:

(12)

To suppress the unphysical pressure fluctuations emerging as artifacts of the scheme given in Eq. 11 we employ an interface pressure correction Praprotnik et al. (2006) within the transition zone. The latter requires a re-parametrization of the effective potential in the system composed of exclusively hybrid molecules (with a constant of w=1/2). We redefine the center-of-mass intermolecular forces as

(13)

where the function is defined as

(14)

so that s[0]=1, s[1]=1, and s[1/4]=0, as shown in Figure 6. The interface correction force is the effective pair force between molecules and defined by the effective pair potential, which is obtained by mapping the hybrid model system composed of exclusively hybrid molecules with to the explicit model system. The corrected effective potential shown in Figure 5 is obtained by mapping the system containing only hybrid molecules. The minima of the effective potential become deeper, and a higher barrier separates the first and second hydration shell when compared to the center-of-mass effective potential . The center-of-mass rdf of the hydrid system exactly reproduces the structure of the reference system as shown in Figure 5.
The mixing function in Eq. 14 can correctly reproduce the thermodynamic properties of the () system with the same mean temperature (0.1 % of difference) and pressure (0.5 % of difference).

Figure 5: The center-of-mass rdfs for explicit, ex-cg(w=1/2), and coarse-grained systems. The inset shows the corrected effective pair potential for hybrid molecules [dashed line] and the reference effective potential for the coarse-grained molecules [full line].
Figure 6: The interface correction weighting function s(w). The value s=1 corresponds to the atomistic and coarse-grained regions and s=0 when w=1/2.

Detailed comparisons between the bulk explicit simulations and the explicit regime in our hybrid setup prove that this approach does not alter the structural properties of the water model studied. In particular, Figure 7 shows that the structural properties of the explicit regime in the multiscale system are exactly the same as in bulk explicit simulations. No change is detected in the orientational preferences of water molecules near the interface region. It is worth mentioning that similar results are obtained in the interface region (results not shown).

Figure 7: The center-of-mass (fp=), OH (fp=) and HH (fp=) rdfs for the explicit region in the hybrid system [dots], and bulk [line] systems.

Figure 8 shows the probability density function of the orientational parameters (O-H vector and the normal of the interface surface) and (H-H vector and the normal of the interface surface) in three consecutive layers of width 4 Å  next to the interface region along the X-axis (see Figure 1 for the definition of the regions).

Figure 8: (Color online) Cosine distribution of the angle formed by the H-H vector (top panel) and O-H vector (bottom panel) of the water molecules with the interface normal vector pointing toward the explicit region in three consecutives regions of width 4 Å  each (see Figure 1 for the definition of the regions).

The uniform cosine distribution of the water molecules in layers inside and next to the interface indicates that the orientational DOFs are fully equilibrated in the hybrid region, as the molecules do not have to reorient upon crossing the interface region and there is no change of behaviour at the interfacePraprotnik et al. (2007c).

Vi Position dependent thermostat

As mentioned above, the structural/thermodynamics of the coarse-grained system are not changed by increasing the background friction in the Langevin thermostat. To obtain the same diffusional dynamics across different resolutions, the coefficient in the Langevin thermostat is changed on-the-fly depending on the number of DOFs of the molecules. As shown in Figure 9, when a constant coefficient is used, two regimes are observed in systems composed of only hybrid molecules:

  • For the value of the diffusion coefficient D is essentially constant;

  • For the value of the diffusion coefficient D drops steeply with .

Figure 9: Top figure: The dashed curve indicates the dependency of the friction coefficient as a function of the particle identity when a position dependent thermostat is used. The full line shows the constant value of the friction coefficient when a regular thermostat is used. Bottom figure: The dots indicate the diffusion of the molecules when the regular thermostat is used. The triangles indicate the diffusion of the molecules when the position dependent thermostat is used.

When the molecular identity is greater than the hybrid molecules start to equilibrate their orientational structure and their diffusive dynamics is slowed down, as shown in Figures 9 and 10.

Figure 10: (Color online) The H-H,O-H and O-O rdfs for different particles identities. The water molecules start to correct their structural properties when .

In order to obtain the same diffusive dynamics for different levels of resolution we propose the following functional form for the friction coefficient in the Langevin thermostat:

(15)

This choice provides a simple interpolation between the two limit values of and . The parameters and are and , respectively. As shown in Figure 9 systems containing only hybrid molecules with different particle identities exhibit the same diffusional dynamics when has the functional form proposed in Eq. 15.

In our simulation we employ a local Langevin thermostat with a position dependent friction Carmeli and Nitzan (1983); Moix and Hernandez (2005); Grigolini (1988); Straus and Voth (1992) to match the diffusion constants of the coarse-grained and all-atom regimes. The Langevin equation with a position dependent coefficient can be written as Carmeli and Nitzan (1983):

(16)

where is:

(17)
(18)

The functional form of depends on the weighting function and it is shown in Figure 9 (see Eq.15).

Figure 11: Top figure: Normalized density profile in the x-direction of the mixed system when the position dependent thermostat is used. Bottom figure: Normalized density profile in the x-direction of the mixed system when a constant friction constant is used.

There is no difference in the structural and thermodynamic properties of the system when a position dependent is used instead of a constant coefficient . The density is homogeneous in both the coarse-grained and explicit regions with (very) small oscillations in the transition regime, cf. Figure 11 .

Figure 12: (a) Time evolution of a diffusion profile when the position dependent thermostat is used for molecules that are initially (at time ps) located inside two neighboring slabs at opposite sides of the mid interface layer. The diffusion profile is averaged over different time origins. Vertical lines denote the boundaries of the interface layer. Bottom figure: The diffusion profile for the all-atom side is averaged over different time origins. Middle figure: The same as before but for molecules that are initially localized inside the slab on the coarse-grained side of the interface region. Top figure: Difference between the corresponding diffusion-profile of the coarse-grained and the all-atom regions. (b) Time evolution of a diffusion profile with the regular thermostat using the same criteria as given in (a).

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 within two slabs of width Å  neighboring the interface layer. The molecules initially localized within the two slabs spread out symmetrically with time when a position dependent friction is used in the local thermostat. On the contrary, the molecules spread out asymmetrically with time as shown in Figure 12 when a fixed coefficient is used. This asymmetry arises from the above-mentioned difference in diffusion coefficient D between the all-atom and coarse-grained regions (see Ref. Praprotnik et al. (2007c)). Figure 13 shows the behaviour of the diffusion coefficient D as a function of the x position using particles that are within a slab of Å  of the mixed system for the two different type of thermostats. As Figure 13 illustrates, the molecules are slowed down when they cross the interface from the coarse-grained region to the explicit region if a constant friction is used for the thermostat. The change in the diffusion coefficient D is localized in the interface region while D is approximately constant inside each region, coarse-grained [] and explicit []. When a position dependent friction is used for the local thermostat the diffusive dynamics of the molecules is the same as for the all-atom system across all regions, as shown in Figure 13.

Figure 13: Diffusion coefficient as a function of the x position using particles that are within a slab of Å  of the hybrid system for different type of thermostats. Each point in the plot is centered in each slab. The average is over different time origins. The horizontal line is the diffusion coefficient of a system containing only all-atom water molecules. The dots indicate the diffusion of the molecules across the hybrid system when the regular thermostat is used. The triangles indicate the diffusion of the molecules across the hybrid system when the position dependent thermostat is used.

Vii Conclusions

In this paper we extended the multiscale model for water we have recently proposed Praprotnik et al. (2007c) to obtain the same diffusional dynamics across the different resolutions. We also studied the accuracy of different coarse-grained water models in reproducing the structural properties of the all-atom system. The results of this study show that for our purposes the single-site model performs equally well as the more sophisticated (three-site and two-site) coarse-grained models and is hence the most appropriate for the presented adaptive resolution simulations. We envisage that such adaptive resolution simulations of water will play an important role in the modeling of biomolecular systems where the coupling of different time and length scales is crucial to understand their physico-chemical properties.

Aknowledgments

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, Texas-ATP, the Robert A. Welch Foundation (C.C.) and the Volkswagen foundation (K.K & L. D. S.). The Rice Cray XD1 Cluster ADA used for the calculations is supported by NSF, Intel, and Hewlett Packard.

Viii Appendix

viii.1 Multiscale simulation protocol

All the results presented in the paper were obtained by performing nVT simulations using the ESPResSo Limbach et al. (2006) (for the one-site and multi-scale model) and AMBER 8 Pearlman et al. (1995) (for the 2-site and 3-site models) simulation packages with a Langevin thermostat, with a friction constant  ps when a regular thermostat is used and a time step of  ps. All the models considered a system of of 1464 water molecules at  K and  g/cm. The density was obtained from an all-atom NPT simulation with  atm using a Reaction-field method for the electrostatics and a cut-off method for the pressure. The results presented for the 2-site and 3-site models were obtained for a density  g/cm. This value of the density was obtained from an all-atom NPT simulation with  atm using an Ewald-method for the electrostatics and a long-ranged correction for the pressure after the cut-off. Periodic boundary conditions and minimum image convention were applied in all directions. The bonds and angle of the water molecules were constrained by using the RATTLE procedure. After warm-up and equilibration, several trajectories of  ns each were collected for each different model. All simulations were performed with a force capping to prevent possible force singularities that could emerge because of overlaps with neighboring molecules when a given molecule enters the interface layer from the coarse-grained side. The size of the system is  Å in the x direction and  Å in both the y and z directions. An interface layer of  Å between the coarse-grained and all-atom models is set along the x direction. To treat all the molecules of the system equally regardless of their level of detail, we define the pressure of the system according to their lower level of detail, that is by the molecular pressure Praprotnik et al. (2005, 2006); Berendsen et al. (1984); Akkermans and Ciccotti (2004). The temperature was evaluated using the fractional analog of the equipartition theorem:

(19)

where is the average kinetic energy per fractional quadratic DOF with the weight Praprotnik et al. (2007a). Via Eq. (19) the temperature is also rigorously defined in the transition regime in which the rotational DOFs are partially ’switched on/off’.

For each of the models considered, the following expression:

(20)

is fitted within line thickness to the tabulated effective potential.

viii.2 Three-site interaction model

The coarse-graining procedure described in the paper converges after 12 iterations, yielding the effective potentials shown in Figure 14.

Figure 14: Effective pair potentials for the 3-site water model. The dashed, dashed-dot, and full line curves indicate the O-O, H-H, and O-H interactions, respectively.

viii.3 Two-site interaction model

The coarse-graining procedure converges after 25 iterations; the resulting effective potential is shown in Figure 15.

Figure 15: Effective pair potentials for the 2-site water model. The dashed, dashed-dot, and full line curves indicate the O-O, dipole-dipole, and O-dipole interactions, respectively.

viii.4 One-site interaction model

The coarse-graining procedure converges after 8 iterations. The resulting potential for the effective interaction between the center-of-mass of two molecules is shown in Figure 16.

Figure 16: The effective interaction between the center-of-mass of two molecules as obtained for the single-site water model.

References

  1. S. Pal, J. Peon, and A. Zewail, Proc. Natl Acad. Sci. USA 99, 15297 (2002).
  2. M. Cheung, A. Garcia, and J. Onuchic, Proc. Natl Acad. Sci. USA 99, 685 (2002).
  3. M. Praprotnik and D. Janežič, J. Chem. Phys. 122, 174103 (2005).
  4. J. Q. Broughton, F. F. Abraham, N. Bernstein, and E. Kaxiras, Phys. Rev. B 60, 2391 (1999).
  5. G. Csanyi, T. Albaret, M. C. Payne, and A. DeVita, Phys. Rev. Lett. 93, 175503 (2004).
  6. K. Kremer, in Multiscale Modelling and Simulation, edited by S.Attinger and P. Koumoutsakis (Springer, 2004), Lecture Notes on Computational Science and Engineering.
  7. P. Koumoutsakos, Annu. Rev. Fluid Mech. 37, 457 (2005).
  8. M. Neri, C. Anselmi, M. Cascella, A. Maritan, and P. Carloni, Phys. Rev. Lett. 95, 218102 (2005).
  9. J. A. Backer, C. P. Lowe, and H. C. J. Ledema, J. Chem. Phys. 123, 114905 (2005).
  10. C. F. Abrams, J. Chem. Phys. 123, 234101 (2005).
  11. G. D. Fabritiis, R. Delgado-Buscalioni, and P. V. Coveney, Phys. Rev. Lett. 97, 134501 (2006).
  12. E. Lyman, F. M. Ytreberg, and D. M. Zuckerman, Phys. Rev. Lett. 96, 028105 (2006).
  13. M. Praprotnik, K. Kremer, and L. Delle Site, Phys. Rev. E 75, 017701 (2007a).
  14. M. Praprotnik, K. Kremer, and L. Delle Site, J. Phys. A: Math. Theor. 40, F281 (2007b).
  15. M. Praprotnik, L. Delle Site, and K. Kremer, J. Chem. Phys. 123, 224106 (2005).
  16. M. Praprotnik, L. Delle Site, and K. Kremer, Phys. Rev. E 73, 066701 (2006).
  17. M. Praprotnik, S. Matysiak, L. Delle Site, K. Kremer, and C. Clementi, J. Phys.: Condens. Matter 19, 292201 (2007c).
  18. S. Matysiak, A. Montesi, A. B. Kolomeisky, M. Pasquali, and C. Clementi, Phys. Rev. Lett. 96, 118103 (2006).
  19. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J Chem Phys 79, 926 (1983).
  20. A. P. Lyubartsev and A. Laaksonen, Phys Rev E 52, 3730 (1995).
  21. S. Matysiak and C. Clementi, J. Mol. Biol. 343, 235 (2004).
  22. S. Matysiak and C. Clementi, J. Mol. Biol. 363, 297 (2006).
  23. D. Reith, M. Putz, and F. Müller-Plathe, J Comput Chem 24, 1624 (2003).
  24. K. A. Dill, T. M. Truskett, V. Vlachy, and B. Hribar-Lee, Annu. Rev. Biophys. Biomol. Struct. 34, 173 (2005).
  25. S. Garde and H. S. Ashbaugh, J. Chem. Phys. 115, 977 (2001).
  26. G. A. Voth and S. Izvekov, J. Chem. Phys. 123, 134105 (2005).
  27. T. Head-Gordon and F. H. Stillinger, J. Chem. Phys. 98, 3313 (1993).
  28. A. K. Soper, Chem Phys 202, 295 (1996).
  29. I. Nezbeda, Molec Phys 103, 59 (2005).
  30. A. P. Lyubartsev and A. Laaksonen, Chem. Phys. Lett. 325, 15 (2000).
  31. M. E. Johnson, T. Head-Gordon, and A. A. Louis, J. Chem. Phys. 126, 144509 (2007).
  32. J. R. Errigton and P. G. Debenedetti, Nature 409, 318 (2001).
  33. W. Tschöp, K. Kremer, J. Batoulis, T. Bürger, and O. Hahn, Acta Polym. 49, 61 (1998).
  34. S. Izvekov and G. A. Voth, J. Chem. Phys. 125, 151101 (2006).
  35. K. Kremer and G. S. Grest, J. Chem. Phys. 92, 5057 (1990).
  36. I. G. Tironi, R. Sperb, P. E. Smith, and W. F. van Gunsteren, J. Chem. Phys. 102, 5451 (1995).
  37. T. Soddemann, B. Dünweg, and K. Kremer, Phys. Rev. E 68, 046702 (2003).
  38. M. Praprotnik, L. Delle Site, and K. Kremer, J. Chem. Phys. 126, 134902 (2007d).
  39. H. Goldstein, Classical Mechanics (Addison-Wesley Publishing Company, 1980), 2nd ed.
  40. In principle, the coupling of different scales could also be done on the respective potentials instead of on forces. However, we showed in Refs. Praprotnik et al. (2007a, b); Delle Site (2007) that, although in this case the system is conservative with a defined total potential energy, this leads to a violation of Newton’s Third Law and consequently to a non-conservation of the linear momentum. Very recently, an alternative “energy-conserving” approach was proposed in Ref. Ensing et al. (2007), where the authors proposed a scheme that couples the atomistic and coarse-grained potentials instead of forces. However, the proposed scheme is ultimately the same as AdResS (with a slightly modified weighting function), i.e., Eq.(11) is employed for the force evaluation (satisfying the Newton’s Third Law), as the “spurious” forces associated with the potential coupling are neglected. By doing this the computed energies do not correspond to forces employed in integrating the equations of motion.
  41. L. Delle Site, Phys. Rev. E (in press, see also arXiv:0709.2579v1).
  42. B. Ensing, S. O. Nielsen, P. B. Moore, M. L. Klein, and M. Parrinello, J. Chem. Theory Comput. 3, 1100 (2007).
  43. M. Neumann, Mol. Phys. 50, 841 (1983).
  44. M. Neumann, J. Chem. Phys. 82, 5663 (1985).
  45. M. Neumann, J. Chem. Phys. 85, 1567 (1986).
  46. W. Im, S. Berneche, and B. Roux, J. Chem. Phys. 114, 2924 (2001).
  47. M. Praprotnik, D. Janežič, and J. Mavri, J. Phys. Chem. A 108, 11056 (2004).
  48. B. Carmeli and A. Nitzan, Chem. Phys. Lett. 102, 517 (1983).
  49. M. J. Moix and R. Hernandez, J. Chem. Phys. 122, 114111 (2005).
  50. P. Grigolini, J. Chem. Phys. 89, 4300 (1988).
  51. J. B. Straus and G. A. Voth, J. Chem. Phys. 96, 5460 (1992)
  52. H. J. Limbach, A. Arnold, B. A. Mann, and C. Holm, Comp. Phys. Comm. 174, 704 (2006); http://www.espresso.mpg.de
  53. H. C. Anderson, J. Comput. Phys. 52, 24 (1983).
  54. H. Berendsen, J. Postma, W. V. Gunsteren, A. D. Nola, and J. Haak, J. Chem. Phys. 81, 3684 (1984).
  55. R. L. C. Akkermans and G. Ciccotti, J. Phys. Chem. B 108, 6866 (2004).
  56. D. A. Pearlman, D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatam, D. M. Ferguson, U. Chandra Singh, P. Weiner, and P. A. Kollman, AMBER, V. 4.1 (1995).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
63977
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description