Adaptive resolution molecular dynamics simulation through coupling to an internal particle reservoir
Sebastian Fritsch, Simon Poblete, Christoph Junghans, Giovanni Ciccotti, Luigi Delle Site, Kurt Kremer
aa r X i v : . [ phy s i c s . c o m p - ph ] N ov Adaptive resolution molecular dynamics simulation through coupling to an internalparticle reservoir
Sebastian Fritsch, ∗ Simon Poblete, ∗ Christoph Junghans, Giovanni Ciccotti, Luigi Delle Site, † and Kurt Kremer ‡ Max Planck Institut f¨ur Polymerforschung, Ackermannweg 10, D-55128 Mainz, Germany EMPS Building, Room 302/B UCD, School of Physics, University College Dublin,Belfield, Dublin 4, Ireland and Dipartimento di Fisica and CNISM,Universit´a “La Sapienza”, Piazzale Aldo Moro 2, 00185 Rome, Italy
For simulation studies of (macro) molecular liquids it would be of significant interest to be ableto adjust or increase the level of resolution within one region of space, while allowing for the freeexchange of molecules between open regions of different resolution or representation. We generalizethe adaptive resolution idea and suggest an interpretation in terms of an effective generalized grandcanonical approach. The method is applied to liquid water at ambient conditions. PACS numbers:02.70.Ns, 61.20.Ja, 80.82, 05.10.-a, 05.20.Gg
Many (macro-) molecular systems display phenomenaand properties that are inherently multiscale in nature,and thus particularly challenging for both theoretical andexperimental investigations. Here computer simulationsrepresent a powerful tool of investigation, though trulylarge all-atom simulations of complex molecular systemsare neither possible, because of CPU requirements, norin many cases very useful, as they often produce an ex-cess of data. To the aim of linking large scale genericproperties and local specific chemistry, a variety of scalebridging simulation techniques, ranging from hierarchi-cal parameterizations of models covering different levelsof resolution, see e.g. [1–7], to interfaced layers of differ-ent resolutions (see e.g. [8–18]) have been developed.As most techniques are sequential – the whole systemis treated on one level of resolution at a time – switch-ing between resolutions may be performed only for thewhole system. In many cases however it would be desir-able to adjust the level of resolution on the fly only ina smaller well defined region of space, i.e. consideringmore details, according to the problem of interest, whilekeeping the larger surrounding on a coarser, computa-tionally much more efficient level. This idea has beensuccessfully employed for problems of crack propagationin hard matter [19]. However for soft matter or liquidsinherent molecular number fluctuations pose special chal-lenges. They require a simulation setup, which allowsmolecules or parts of large molecules to cross bound-aries between areas of different representation withoutany barrier, keeping the overall thermodynamic equilib-rium intact. While crossing between regions of differ-ent resolution, they either become more coarse-grainedor resume more atomistic detail that is they acquire orlose degrees of freedom (DOFs). In this context, adap-tive resolution MD simulations would represent a naturalGrand Canonical set up: An atomistic region/cavity iscoupled to a large, ideally infinite coarse-grained regionacting as a particle/molecule reservoir. In general forsuch a coupling scheme the only requirement is that the coarse grained surrounding supplies/removes moleculesand heat in a way that in the atomistic region µ ,V andT remain constant. The adaptive resolution method(AdResS, Adaptive Resolution Scheme), which allowsone to couple two systems with different resolution withinone MD simulation [20–24], describes a first step in thisdirection. Molecules change their resolution and thustheir number of DOFs when moving through space andpassing from one region through a transition zone intothe other region. The method has successfully been ex-tended into the quantum [25] and the continuum [26, 27]region. A different approach based on the interpolationof Lagrangians can be found in Ref. [28] In the follow-ing we will generalize the AdResS concept and show thatthe atomistic/ fine grained subsystems can be viewed asa system in contact with a (huge) coarse-grained par-ticle/molecule reservoir, thus effectively representing a µ VT ensemble on the level of the high resolution sub-system, while the overall systems remains to be an NVTensemble. Actually, by linking the coarse grained regionin addition to a continuum, as shown in Refs. [26, 27],this is easily extended to a µ VT ensemble on the level ofthe whole system. We now first shortly review the adap-tive resolution method AdResS and the concept of ther-modynamic force which assures overall thermodynamicequilibrium. We then extend this to a more general cou-pling scheme and apply it to the example of an adaptiveresolution simulation of liquid water, where the coarse-grained water model matches the compressibility of theall-atom model, while the pressure is higher by three or-ders of magnitude.In the AdResS method molecules smoothly changetheir level of resolution as a function of position by mov-ing through a transition region as illustrated in Fig.1.In order to accomplish the process described above, twodifferent representations are interpolated by a transi-tion function w ( x ) with w = 1 in region A , the finegrained more microscopic or all-atom one, and w = 0in the coarse-grained region B , respectively. Further Atomistic Hybrid Coarse-grained a x w ( x ) b FIG. 1: Pictorial representation of the adaptive simulationbox and local molecular representation. On the left, is thelow resolution (coarse-grained) region, indicated by B , and onthe right, is the high resolution (atomistic) region A. In themiddle is the transition (hybrid) region H with the switchingfunction w ( x ) (curve in grey), let us assume that a particle in B is the coarse-grainedrepresentation of a collection of atoms in region A andthe position of the coarse particle is the center of massof the corresponding atomistic molecule. With thisin mind two molecules α, β in region A (fine-grained, w ( X α ) = w ( X β ) = 1, X α , X β being the x-coordinatesof the centers of mass of the two molecules, cf Fig. 1)interact via the standard forces between all individualatoms of the molecules; in region B (coarse-grained, w ( X α ) = w ( X β ) = 0) two coarser molecules interactwith each other via central forces; just as in a bulk sys-tem of all atom molecules or coarse grained molecules,respectively. In all other cases of mixed resolution in-teractions, these two forces are interpolated as illustratedin Fig. 1 so that the total force F αβ between two differentmolecules consist of two parts. The all-atom contributionreads: F Aαβ = w ( X α ) w ( X β ) X i,j F Aα i β j α = β (1)While the coarse-grained parts is: F Bαβ = (1 − w ( X α ) w ( X β )) F Bαβ α = β (2) F Aα i β j is the force between atom i of molecule α andatom j of molecule β and F Bαβ is the coarse-grained cen-ter of mass – center of mass interaction. Note that inthe transition region the intra-molecular atom-atom in-teractions are left unaffected. For technical details ofthe implementation we refer to the supporting mate-rial. Altogether we simulate a system at constant vol-ume, partitioned in two different (ideally macroscopic)regions (coarse-grained and fine-grained) by a transitionregion of finite width. This region facilitates the exchangeof molecules and ensures equilibrium between the coarseand fine grained regions, moreover, by construction (as explained below), should not affect the density fluctua-tions (compressibility). Beyond the crucial technical taskof assuring the thermodynamically correct coupling be-tween the two large regions, the transition region doesnot have, a priori, any thermodynamic/physical meaning.An important characteristic of the adaptive approach isthat while the overall number of molecules is conserved,this is not the case for the number of atoms in the in-dividual regions. Eventually we aim for a flat densityprofile throughout the whole simulation box, includingthe transition region, so that we have a thermodynamicstate point well defined without (kinetic) barriers for theexchange of molecules.In the same spirit, the temperature has to be well de-fined and equal throughout the system. However, whilethis is simple to achieve in the ”bulk” regions, the tran-sition region requires some care since in this region achange in the number of degrees of freedom occurs. Byswitching from coarse-grained to atomistic or vice versa,degrees of freedom (DOFs) are reintroduced or removedand the related equilibration is accomplished by a ”local”thermostating procedure (for details see the supportingmaterial). The approach presented above, paves the wayto an open system simulation scheme within a generalizedGrand Canonical framework.A common method to derive coarse grained interac-tions is based on matching the center of mass radialdistribution functions between the molecules in the twosystems. As a result the compressibilities κ remainunchanged [29] by construction while the pressures inthe fine and coarse-grained representation can be differ-ent [30]. In terms of the grand potential − pV this reads(for identical volumes) p A ( µ A , T ) V = p B ( µ B , T ) V ; κ A = κ B , (3)where µ A and µ B are the respective chemical poten-tials. This will result in density variations in the sys-tem. Instead we would like to have a constant density ρ throughout the whole simulation box, while keeping iden-tical compressibilities. Different pressures at the bound-aries of the transition region will produce a drift forceon the molecules, leading to an unphysical inhomoge-neous system. Different compressibilities would affectthe molecule number fluctuations in the - still finite -regions. Thus we keep κ A = κ B and propose to exactlycompensate this pressure generated drift force by a ther-modynamic force F th ( x ) acting on the molecules and Eq.3 can then be rewritten as (see Fig. 1) p A ( µ A , T ) + ρ M α Z ba F th ( x )d x ! V = p B ( µ B , T ) V (4)with ρ being the reference, uniform, molecular density.Eq.4 tells that the subsystem A ( B ) is in equilibriumwith a reservoir µ B ( µ A ) respectively, notwithstandingthe fact that the pressures, p A and p B , and local chemi-cal potentials, µ A and µ B , may be different. To overcomethese variations, the related thermodynamic work, whichexactly compensates any related drift force between thedifferent regions, is provided/absorbed by F th . Let theall-atom system, at a given density ρ , serve as refer-ence system, then the force on a molecule with mass M α balancing −∇ p ( x ), reads F th ( x ) = M α ρ ∇ p ( x ) (5)where p ( x ) is the local pressure at the target density ρ . F th is applied for a < x < b , where a ( b ) is theleft (right) boundary of the transition region. In Eq.5, p ( x = a ) = p A and p ( x = b ) = p B . Now the density willremain unchanged. With this generalization the originalequation for the coarse-grained forces in the transitionregion, Eq. 2 is changed so that the force ˆ F Bα acting onthe center of mass of molecule α reads:ˆ F Bα = X β F Bαβ + F th ( X α ) (6)The above line of arguments does not require the pres-sures in the two (”macroscopic”) regions, p A and p B ,to be the same, allowing to couple rather different sys-tems. Formally this simply means that the thermody-namic force as well as the thermostat can perform workon the molecules to adjust the virial pressure and thethermal energy, respectively, while passing from one re-gion to the other. This procedure allows one to couplealmost arbitrary systems and keep them in equilibriumwith respect to each other. For typical coarse-grainingprocedures like (iterative) Boltzmann inversion or Re-verse Monte Carlo which are mostly based on radial dis-tribution functions and reproduce structural rather thanthermodynamic aspects very well, this is of special ad-vantage.Unfortunately the optimized pressure profile p ( x ) is notdirectly accessible as it has to be measured under theconstraint of enforcing the equilibrium density ρ in thesystem. Thus we developed an iterative approach basedon the density profile. Within a linear approximation afirst estimate of the local pressure p ( x, ρ ( x )) for an en-forced overall constant density ρ is: p ( x, ρ ( x ))) ≈ p A + 1 ρ κ T ( ρ − ρ ( x )) (7)where κ T is the (constant) isothermal compressibility,and ρ ( x ) is the, non uniform, density to which the systemwould adjust itself to without applying an external force.In order to enforce the uniform density ρ , Eq. 7 andEq. 5 allow to obtain an initial thermodynamic force,that can be refined iteratively through F i +1th ( x ) = F i th ( x ) − M α ρ κ T ∇ ρ i ( x ) , (8) until the system evolves to the target uniform density.Here for i = 1, we have ρ ( x ) = ρ ( x ) and F ( x ) = 0.The prefactor 1 /ρ κ T can be interpreted as a variationof local chemical potential, by the identity [29] (cid:18) ∂µ∂ρ (cid:19) V,T = 1 ρ κ T , (9)a fact explored already before in this context [23]. In thefollowing we will demonstrate the power of this conceptfor the illustrative case of liquid water, where the com-pressibility of the SPC/E water and a structure basedcoarse-grained model are the same, while the pressure inthe coarse-grained system is 6000 bar[31].Coarse graining liquid water with a focus on the liquidstructure usually leads to models for which the pressureis significantly higher while the compressibility remainsessentially unchanged. Thus pressure correction termsare often used, which add linear terms to the coarse-grained potential. The compressibility is then affected as κ T ∝ R r ( g ( r ) − r and small deviations at high r val-ues can have dramatic effects [31]. However, for adaptiveresolution methods, but also solvent-solute systems, it isimportant that the compressibility in the coarse-grainedand all-atom regions is the same. This minimizes finitesize effects for our (in most cases on purpose as smallas possible) atomistic region, imposed by the - still fi-nite - CG region. For many applications in computa-tional biology, but also for studies of water itself, theSPC/E model is employed [32]. A standard structurebased coarse-grained model perfectly fits the radial dis-tribution functions, but does not comply with the tetra-hedral packing [31, 33]. It has -within the error bars- thesame compressibility but a very high pressure of about6000 bar, making it an ideal test case for the thermody-namic force principle described above. If not properly ac-counted for, this pressure difference between the differentmodels would lead to density distortions throughout thesimulation box and create a transient flow of moleculesfrom the coarse-grained to the atomistic region. Using κ AT at the outset of the iteration scheme of Eq. 8 asdescribed before, after only four iterations for the ther-modynamic force we arrive at a flat density profile, asillustrated in Fig. 2.The thermodynamic force creates a stationary statewith the correct target density, removing both the effectof the pressure difference of the models on the distribu-tion of molecules in either region and the artifacts fromthe change of resolution in the transition region. Theadvantages are also clearly displayed, when looking atlocal density fluctuations. While in the thermodynamiclimit the molecule number fluctuations are related to thecompressibility through ρk B T κ T = h N i−h N i h N i , we hereanalyze an effective compressibility in terms of moleculenumber fluctuations in slabs along the axis of resolutionchange. The increased compressibility for a pressure cor- A H B(c) x [nm] < N > − < N > < N > (b) p [ b a r ] (a) ρ / ρ AdResS p cg = p at AdResS κ cg = κ at full atomisticlocal p xx ρ /M α R xa F th d x ′ FIG. 2: (a) Density profile for the water simulation for the un-corrected AdResS simulation of water and several iterations ofthe thermodynamic force with A being the all-atom SPC/Eregion and B the coarse-grained region. The final densityprofile shows deviation of less than 0.4 %. (The arrow indi-cates the development of the density with the iteration.) (b)Compensation of the local pressure change between regionsby the thermodynamic force. (c) Molecule number fluctua-tions in slabs (width ∆ x = 1nm) along the x-axis for pressurecorrected and uncorrected coarse-grained models. rected coarse-grained SPC/E water model leads to an in-crease in the local density fluctuations. Within our setupwe do not alter the effective compressibility anywhere inthe system as shown in Fig. 2c. Actually the moleculenumber fluctuations in subvolumes of region A , are thesame as in similar subvolumes of a full all-atom simula-tion. Hence the particle distribution in region A showsthe expected Gaussian distribution as it is the case in thefull all-atom simulation, shown in Fig. 3. Leveling outany pressure undulations in the transition region, whilekeeping the compressibilities the same also means thatthere is no barrier for molecule exchange between regions(see supporting information). Coupling two systems withdifferent pressures means that the thermodynamic forceperforms work on the molecules, when they cross from N p ( N ) AdResSAll-atomGaussian approximation
FIG. 3: Particle number distribution in a slab of width l =4 nm and for an atomistic zone in adaptive resolution with thesame width. Both distributions match the expected Gaussianbehavior with mean µ = h N i and variance σ = ρk B T κ T h N i one region to the other. In practice the integral over thethermodynamic force over the transition region does nothave to be zero, which would be the case if the pressureson either side were the same. Thus for the present setupthe thermodynamic force has to exactly compensate thelocal pressure variation due to the constant density. Thecomparison of the final local pressure (see supporting in-formation for details) and the thermodynamic force asderived in form of an external local pressure field is shownin Fig. 2b. They perfectly match, illustrating the over-all consistency of the approach. Thus by applying thethermodynamic force we arrive at a robust and efficientsimulation setup, which enables us within a single MDsimulation to keep rather different systems with quitedifferent intrinsic pressures but identical compressibili-ties in equilibrium with each other.The fact that the number of molecules at a given res-olution fluctuates suggests to view each subsystem as areservoir of molecules for the other. Strictly speaking,to interpret such a setup as a Grand Canonical ensem-ble would require macroscopic reservoirs. In our case itis obvious that the number of molecules, and thus thereservoir, is finite. In fact, we have shown that the av-erage and the fluctuating number of molecules in eachsubsystem and even in smaller subregions throughoutthe whole simulation box are the same as analyzed ina full atomistic system. For the general case where theall-atom as well as the coarse-grained system are at thesame (or even well defined different) densities the sim-ulation setup leads to a situation, where a difference inthe excess chemical potential µ (the temperature partfor the internal degrees of freedom is taken care of bythe thermostat) between the regions, is compensated bythe work performed by the thermodynamic force. Thuson the level of the (large) subsystems we run a ( µ, V, T )MD simulation, while we have a ( N, V, T ) ensemble glob-ally. The extreme example of liquid SPC/E water attemperature 300K and pressure 1 bar, where the coarse-grained model at the same temperature and density hasa pressure of more than ≈ Acknowledgment
We thank B. D¨unweg and M.Praprotnik for valuable discussions. We thank M. De-serno and D. Donadio for comments on the manuscript.Part of the work has been performed with the MultiscaleMaterials Modeling Initiative of the Max Planck Society.S. P. thanks the DAAD-Conicyt grant for financial sup-port. G.C. wish to acknowledge financial support fromSFI Grant 08-IN.1-I1869 and from the Istituto Italianodi Tecnologia (IIT) under the seed project grant 259,SIMBEDD. C.J. was financially supported by SFB 625. ∗ Shared First Authorship: S. Fritsch & S. Poblete † Electronic address: [email protected] ‡ Electronic address: [email protected][1] L. Delle Site, C. F. Abrams, A. Alavi, and K. Kremer,Phys.Rev.Lett. , 156103 (2002).[2] L. Delle Site, S. Leon, and K. Kremer,J. Am.Chem.Soc. , 2944 (2004). [3] L. M. Ghiringhelli, B. Hess, N. F. A. van der Vegt, andL. Delle Site, Am. Chem. Soc. , 13460 (2008).[4] K. Reuter, D. Frenkel, and M. Scheffler, Phys. Rev.Lett. , 116105 (2004).[5] J. Rogal, K. Reuter, and M. Scheffler, Phys. Rev. B ,155410 (2008).[6] K. Tarmyshov and F. M¨uller-Plathe, J. Chem. Phys. ,074702(2007).[7] In Coarse Graining of Condensed Phase and Biomolecu-lar Systems , edited by Gregory A. Voth, (Chapman andHall/CRC Press, Taylor and Francis Group,2008).[8] G. Lu, E. B. Tadmor, and E. Kaxiras, Phys. Rev. B ,024108 (2006).[9] J. Rottler, S. Barsky, and M. O. Robbins, Phys. Rev.Lett. , 148304 (2002).[10] G. Csanyi, T. Albaret, M. C. Payne, and A. D. Vita,Phys. Rev. Lett. , 175503 (2004).[11] A. Laio, J. VandeVondele, and U. R¨othlisberger,J. Chem. Phys. , 6941 (2002).[12] D. E. Jiang and E. A. Carter, Acta Materialia , 4801(2004).[13] G. Lu and E. Kaxiras, Phys. Rev. Lett. , 155501(2005).[14] E. Villa, A. Balaeff, L. Mahadevan and K. Schulten, Mul-tiscale Model. Simul. , 527 (2004).[15] M. Neri, C. Anselmi, M. Cascella, A. Maritan and P.Carloni, Phys. Rev. Lett. , 218102 (2005).[16] R. Delgado-Buscalioni and P. Coveney, Phys. Rev. E ,046704 (2003).[17] P. Koumoutsakos, Annu. Rev. Fluid Mech. , 457(2005).[18] E. Lyman, F. Ytreberg and D. Zuckerman, Phys. Rev.Lett. , 028105 (2006).[19] R. E. Rudd and L. Q. Broughton,Phys. Stat. Solid. B , 251 (2000).[20] M. Praprotnik, L. Delle Site, and K. Kremer, J. Chem.Phys. , 224106 (2005).[21] M. Praprotnik, L. Delle Site, and K. Kremer, Phys.Rev. E , 066701 (2006).[22] M. Praprotnik, L. Delle Site, and K. Kremer, Annu. Rev.Phys. Chem. , 545 (2008).[23] S. Poblete, M. Praprotnik, K. Kremer and L. Delle Site,J. Chem. Phys. , 114101 (2010).[24] C. Junghans and S. Poblete, J. Comp. Phys. Comm. ,1449 (2010)[25] A. Poma and L. Delle Site, Phys. Rev. Lett , 250201(2010).[26] R. Delgado-Buscalioni, K. Kremer, and M. Praprotnik,J. Chem. Phys. , 114110 (2008).[27] R. Delgado-Buscalioni, K. Kremer, and M. Praprotnik,J. Chem. Phys. , 244107 (2009).[28] A. Heyden and D. G. Truhlar, J. Chem. theory Comput. , 217, (2008)[29] In Theory of Simple Liquids , J.P. Hansen and I. R.McDonald, Academic Press, 2006.[30] S. Izvekov and G. Voth J. Chem. Phys. , 134105(2005).[31] H. Wang, C. Junghans and K. Kremer, Euro. Phys. J. E , 221 (2009).[32] H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma,J. Phys. Chem. , 6269 (1987).[33] J. R. Errington and P. G. Debenedetti, Nature , 318(2001). r X i v : . [ phy s i c s . c o m p - ph ] N ov Supporting Material: Adaptive resolution molecular dynamics simulation throughcoupling to an internal particle reservoir
Sebastian Fritsch, ∗ Simon Poblete, ∗ Christoph Junghans, Giovanni Ciccotti, Luigi Delle Site, † and Kurt Kremer ‡ Max Planck Institut f¨ur Polymerforschung, Ackermannweg 10, D-55128 Mainz, Germany EMPS Building, Room 302/B UCD, School of Physics, University College Dublin,Belfield, Dublin 4, Ireland and Dipartimento di Fisica and CNISM,Universit´a “La Sapienza”, Piazzale Aldo Moro 2, 00185 Rome, Italy
I. FORCE CALCULATION
In Eqn. 1 and Eqn. 2 of the original paper the inter-actions between molecules α and β subject to the inter-polation of forces is written in a rather general form as F Bαβ = (1 − w ( X α ) w ( X β )) F Bαβ (1)for the center of mass force and F Aαβ = w ( X α ) w ( X β ) X i,j F Aα i β j α = β (2)for the atom level interactions. From there one arrives atthe interpolation formula for the total force between twomolecules used in previous publications: F αβ = w ( X α ) w ( X β ) F Aαβ + (1 − w ( X α ) w ( X β )) F Bαβ (3)From this it is possible to see that, (i) the force F Aαβ in the hybrid region is actually a function of all theatomic coordinates { r i } i ∈ α and { r j } j ∈ β and hence neednot to be in the direction of the connection vector R αβ of the centers of mass and can lead to e.g. rotationof the molecule. (ii) Since the intra molecular forcesdo not contribute to the total force on a molecule, wedo not have to include the intramolecular interactionsin the interpolation scheme; they are kept acting onthe individual atoms throughout the atomistic andtransition regions just as they are.In general Eqns. 1,2 describe four different cases forthe forces between two different molecules, which applyto the whole simulation box:(i) For w ( X α ) = w ( X β ) = 0 two coarse-grained moleculesinteract via the effective coarse-grained interaction F αβ = F Bαβ , (4)leading only to a translational motion of the moleculeswhich are now reduced to point particles. In this regime ∗ Shared First Authorship: S. Fritsch & S. Poblete † Electronic address: [email protected] ‡ Electronic address: [email protected] there is no intermolecular or intramolecular atomisticcontribution.(ii) For w ( X α ) = w ( X β ) = 1 both molecules are in theatomistic region. Thus atom i of molecule α interactswith all atoms j of molecule β F α i ,β = X j ∈ β F Aα i β j . (5)The situation is that of a standard all-atom simulationwith forces acting on each atom.(iii) For 0 < w ( X α ) < w ( X β ) = 0 at least molecule α is in the transition regime and β is not in the coarse-grained region. Then both parts of the interpolation con-tribute to the motion of the molecule, hence the coarse-grained part is also contributing to the force between theatoms. Following Eqns. 1,2 the force on atom α i origi-nating from all atoms of molecule β reads: F α i ,β = w ( X α ) w ( X β ) X j ∈ β F Aα i β j +(1 − w ( X α ) w ( X β )) m i M α F B αβ . (6)The first part describes the standard atom atom inter-actions, while the second contains the distribution ofthe coarse grained interaction onto the individual atoms.By this the coarse grained interactions are naturally in-corporated into the all atom resolution, which we keepthroughout the whole transition region.(iv) For w ( X α ) = 0 and w ( X β ) = 0 only the second partof Eqn. 6, the center of mass contribution, survives andhas to be distributed to the atoms of molecule α , whichgives F α i ,β = m i M α F B αβ (7)for atom α i .With this scheme all forces between atoms and moleculesare well defined throughout the whole simulation box.As mentioned in the main text, the thermostat needssome attention. Kinetic energy and the related tem-perature are both well defined in the coarse-grained andthe atomistic region. More attention is required in thetransition region. By switching from coarse-grained toatomistic or vice versa, degrees of freedom (DOFs) (in-cluding their related thermal energies) are introduced orlost. The related contribution to the thermal energy willbe furnished by a standard local thermostat as has beendiscussed in Refs. [7–10]. Thus the molecules/particlesare ”handed over” properly thermalized at the bound-aries of the transition region. In summary, one can ap-ply any standard local thermostats like, e.g, a Langevinthermostat, which automatically take care of the ”latentheat”[8, 9] related to the fading DOFs. AdResS simu-lations from the point of view of the whole system cor-respond to simulations of a N V T ensemble, where N counts the number of molecules, but not the number ofDOFs (which in our case are not strictly but only in av-erage proportional). In the current implementation weuse a Langevin thermostat, which acts without any in-terpolation on the atoms. The friction is proportionalto the masses of the atoms and thus by themostattingeach atom, the center of mass of every molecule obeysthe fluctuation-dissipation-theorem as well [7]. A purecenter of mass thermostat would also be sufficient, butslightly more difficult to implement in the integrator.The consequences for the definition of kinetic energy etc.that follow from ”partial” degrees of freedom have beenextensively discussed in [8].The above discussion considers molecules which are re-placed by their centers of mass in the coarse-grained re-gion. However a more general mapping from explicitto coarse-grained, like a atomistic polymer mapped to achain of spheres (so called ”coarse-grained” beads) con-nected by springs, is possible. Then the interpolation isapplied to the coarse-grained beads. II. PRESSURE CALCULATION
We can numerically analyze the local pressure in thestationary state and its behavior when an external fieldis applied making use of the method proposed by Todd,Evans and Daivis [1]. The time averaged local pressuretensor is given by¯ p µν ( x ) = 12 A ∆ x (cid:28) X { α : | X α − x | < ∆ x } ( v α · ˆ e µ )( v α · ˆ e ν ) M α (cid:29) t (8)+ 12 A (cid:28) X α F α ˆe µ sgn( X α · ˆe x − x ) (cid:29) t (9)where A is a cross-sectional area perpendicular to the x -direction. The first term corresponds to the ideal gascontribution calculated in a small slab centered at x withthickness 2∆ x , and the second represents the contribu-tion of the forces. Both sums run over molecules (whichare labeled with α as before), where M α , v α and X α aremass, velocity and position. Consequently, the ideal gascontribution of the center of mass is fully included in thesum regardless of the weighting function of the molecule,since those degrees of freedom fully contribute to the av-erage in any representation. The internal degrees of free-dom, however, are not considered here. Therefore, theresolution of each molecule is implicit in the force con-tribution that depends on the weighting function. In the case of two-body forces, the second sum is reduced tothe count of the forces that cross the area A . Particlesfar away from this plane (at distances bigger than themaximal cutoff of the non-bonded interactions) do notcontribute to the sum. III. SIMULATION PROTOCOL
The simulated system consisted of 8507 watermolecules in a box of dimensions 16 nm × × κ atT = 4 . × − m /N [4].The coarse-grained potential was obtained by mappingeach water molecule in the center of mass to a singlecoarse-grained bead by an iterative Boltzmann inversionbased on radial distribution functions. 100 steps of 0.2nseach were carried using the VOTCA software package [5].As mentioned before, such a coarse-graining does not pre-serve the state point, yielding a coarse-grained pressureof approximately 6200 bar, while the compressibility isnot altered. That is the local slope of the equation ofstate is conserved by such a procedure, while the abso-lute value is shifted. Additional corrections to reproducethe atomistic pressure can be applied [6], but have beenshown to alter the compressibility [4] by a factor of 5,which one has to avoid in order to obtain the correct par-ticle number fluctuation in the different regions of space.To obtain the highest accuracy in the final density pro-file the thermodynamic force was extended by half thecut-off length (0.6 nm) into the atomistic zone (becauseof the mixed atomistic coarse grained interaction of twomolecules, when one molecule is in the transition region).Without this extension the density deviated at the mostby about 1 % from the target density, showing a drop inthe atomistic zone. This is due to the fact that close tothe coarse-grained zone, the atomistic molecules interactwith the hybrid molecules an thus have a slightly alteredstate point. IV. DIFFUSION PROFILE
To proof that molecules indeed diffuse back and forthbetween both regions without regard to the local levelof resolution, we can measure diffusion profiles along theaxis of resolution change. At t = t all particles in a slabof width ∆ x = 1nm are selected which are then ”traced”(in terms of their density profile) over time (fig. 1). Thisis done for a set of particles in both regions. As shownby the symmetry of the decay for both cases, the tran-sition region does not introduce any further transitionbarrier. On the coarse-grained side the particle distribu- BHA x [nm] ρ / ρ [ ] FIG. 1: Interdiffusion between the two regions. This showsthe density profiles over time for particles originally in a cubicsubvolume with width ∆ x = 1nm at t = t . The profiles areasymmetric as diffusion constants are different in each repre-sentation. By this it is clear that both domains constantlyexchange particles. tion decays significantly faster, illustrating the enhanceddiffusion constant in this region. [1] B. D. Todd, D. J. Evans and P. J. Daivis (1995) Pressuretensor for inhomogeneous fluids , Phys. Rev. E .[2] B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl(2008) GROMACS 4: Algorithms for Highly Efficient,Load-Balanced, and Scalable Molecular Simulation , J.Chem. Theory Comput. , 435.[3] I. G. Tironi, R. Sperb, P. E. Smith and W. F. vanGunsteren (1995) A generalized reaction field method formolecular dynamics simulations , J. Chem. Phys. ,5451.[4] H. Wang, C. Junghans and K. Kremer (2009)
Compar-ative atomistic and coarse-grained study of water: Whatdo we lose by coarse-graining? , Euro. Phys. J. E , 221.[5] V. R¨uhle, C. Junghans, A. Lukyanov, K. Kremer and D.Andrienko (2009) Versatile Object-Oriented Toolkit forCoarse-Graining Applications , J. Chem. Theory Comput. , 3211. [6] D. Reith, M. P¨utz and F. M¨uller-Plathe (2003) Derivingeffective mesoscale potentials from atomistic simulations ,J. Comput. Chem. , 1624.[7] C. Junghans and S. Poblete (2010) A reference implemen-tation of the adaptive resolution scheme in ESPResSo , J.Comp. Phys. Comm. , 1449.[8] M. Praprotnik, K. Kremer, L. Delle Site (2007)
Frac-tional dimensions of phase space variables: a tool forvarying the degrees of freedom of a system in a multi-scale treatment , J. Physics A , F281.[9] M. Praprotnik, K. Kremer, L. Delle Site (2007) Adaptivemolecular resolution via a continuous change of the phasespace dimensionality , Phys. Rev. E. , 017701[10] M. Praprotnik, L. Delle Site and K. Kremer (2008) Mul-tiscale Simulation of Soft Matter: From Scale Bridging toAdaptive Resolution , Annu. Rev. Phys. Chem.59