Non-conformal coarse-grained potentials for water
NNon-conformal coarse-grained potentials for water
Tonalli Rodríguez-López,
1, 2, a) Yuriy Khalak, b) and Mikko Karttunen
3, 4, 5, c) Departamento de Física, Universidad Autónoma Metropolitana, Iztapalapa, Apdo 55 534, México DF 09340,México. Department of Chemistry & Waterloo Institute of Nanotechnology, University of Waterloo,200 University Avenue West, Waterloo, Ontario, Canada N2L 3G1 Department of Mathematics and Computer Science & Institute for Complex Molecular Systems,Eindhoven University of Technology, P.O. Box 513, MetaForum, 5600 MB Eindhoven,the Netherlands Department of Chemistry, Western University, 1151 Richmond Street, London, Ontario,Canada N6A 5B7 Department of Applied Mathematics, Western University, 1151 Richmond Street, London, Ontario,Canada N6A 5B7 (Dated: 4 September 2018)
Water is a notoriously difficult substance to model both accurately and efficiently. Here, we focus on de-scriptions with a single coarse-grained particle per molecule using the so-called Approximate Non-Conformal(ANC) and generalized Stockmayer potentials as the starting points. They are fitted using the radial dis-tribution function and the liquid-gas density profile of the atomistic SPC/E model by downhill simplexoptimization. We compare the results with monatomic water (mW), ELBA, as well as with direct IterativeBoltzmann Inversion (IBI) of SPC/E. The results show that symmetrical potentials result in non-transferablemodels, that is, they need to be reparametrized for new state-points. This indicates that transferabilitymay require more complex models. Furthermore, the results also show that the addition of a point dipole isnot sufficient to make the potentials accurate and transferable to different temperatures (300 K-500 K) andpressures without an appropriate choice of properties as targets during model optimization.
I. INTRODUCTION
Water is arguably the most important and studied sub-stance in the world. It covers approximately 70% of theEarth’s surface and is vital for all known forms of life. Inparticular, at the microbiological scale, membranes thatprovide the container and selective protection for all cells,as well as proteins, peptides, and drugs, all require anaqueous solution to function.
Despite being the most studied substance, our under-standing of water is far from comprehensive. Its physicalbehavior and interactions with biological molecules keep providing surprises. Developing computational andtheoretical models for water has proven to be extremelydemanding and remains a topic of intense research .The focus of this work is on coarse-grained (CG) single-particle models of water. The general goal of coarse-graining is to reduce the number of degrees of free-dom allowing for significantly longer simulation timesand larger systems yet retaining the essential physicalproperties . There is no single approach to coarse-graining and a large number of CG models with uniquestrengths and weaknesses exist . Apart from implicitsolvent models, each of the approaches has their ownmodel(s) for water . Here, we operate at a modestlevel of coarse-graining to be able to make direct com-parisons with the underlying atomistic description. a) E-mail: [email protected] b) E-mail: [email protected] c) E-mail: [email protected]
The reduction in degrees of freedom and averaging overproperties such as electrostatic interactions and hydrogenbonding can sometimes lead to unwanted and surprisingeffects. Observed effects include freezing of water athigher temperatures (this occurs with the original wa-ter of the very popular MARTINI model; water freezesat ± K but the polarizable MARTINI water modellater corrected this to lower temperatures between 280and 285 K ), freezing and ordering of CG water closeto surfaces or the tradeoff between reproducing thestructure factor and tetrahedral packing with single wa-ter models . Each of the coarse-grained (as well asatomistic) water models have their strengths and weak-nesses, and it is important to pay attention that the cho-sen model is approriate for the intended range of phys-ical parameters such as temperature or pressure. Darre et al. , and Hadley and McCabe provide recent reviews ofcoarse-grained water models .To overcome some of these issues, we introduce threecustomized models for CG water based on the Approx-imate Non-Conformal (ANC) potential . These newmodels, called ANC-1G, ANC-2G and GSM, will be de-tailed below. We characterize them and compare theirbehavior to some of the well-known CG water models,namely monatomic Water (mW) and the electrostatic-based (ELBA) model . As an all-atom referencemodel for water, we use the extended simple point charge(SPC/E) model and apply the iterative Boltzmann In-version (IBI) method to it to obtain a CG reference.Majority of modern water (and biomolecular) models,whether coarse-grained or atomistic, rely on the Lennard-Jones potential. The motivation for diverting from the a r X i v : . [ c ond - m a t . s o f t ] A ug usual Lennard-Jones potential-based approach in favor ofa non-conformal model is that liquid water has alreadylong ago been shown to be poorly described by van derWaals’ law of corresponding states resulting in modelsthat only qualitatively reproduce density and compress-ibility dependence on pressure . As shown by Pitzeralready in 1939, the law of corresponding states holds ex-actly only for fluids with conformal potentials, that is, po-tentials that remain invariant under energy and distancescaling like Lennard-Jones. To describe interactions morerealistically beyond Lennard-Jones, extensions and gen-eralizations of non-conformal models involving more pa-rameters especially in the context of colloids have beendeveloped, see e.g. Refs. 37 and 38. II. MODELS AND METHODSA. Models
1. Atomistic reference model
The reference model, SPC/E water , is characterizedby three point masses with oxygen-hydrogen distance of . Å and the HOH angle equal to the tetrahedral an-gle . ◦ , with charges on the oxygen and hydrogenequal to − . e and . e , respectively, and withthe Lennard-Jones parameters of oxygen-oxygen interac-tion set to (cid:15) LJ = 0 . kcal/mol and σ LJ = 3 . Å. Wechose SPC/E since it is commonly used and reproducesproperties such as isothermal compressibility and criticalpoint temperature better than most other popular three-site models . While a four-site model like TIP4P/2005 does yield more accurate liquid densities , it is lessfrequently used in large-scale simulations.
2. New CG parameterizations
Using the SPC/E model as a reference, weparametrized three new coarse-grained models, ANC-1G,ANC-2G and GSM. In them, water molecules are rep-resented as single particles that interact through eitherisotropic (ANC-1G and ANC-2G) or anisotropic (GSM)interactions. All three are based on the original ANCpair potential of del Rio et al. In the procedure below, we fit the potential twice us-ing the the radial distribution function as the first targetand the liquid-gas density profile as the second. This al-lows us to determine which target property provides abetter model. Structure based fitting as the sole targetis used in methods such as the iterative Boltzmann In-version (IBI) . The hope is that adding a second targetwill help to obtain a better model; density was chosenbecause it is a commonly used metric in the developmentof atomistic water models. We would like to note thatthis approach is different from using multiple simultane-ous targets as is used with, for example, the mW and ELBA models. Fitting to multiple properties at onceproduces models that, at least in principle, approximatemore properties at once, but do not necessarily exactlymatch any. Thus, fitting separate models to individualproperties and fitting a single model to multiple proper-ties simultaneously are two different philosophies. Theresults here provide a comparison between the two ap-proaches.The ANC potential retains the length and energyparameters from Lennard-Jones, and enforces non-conformality through a third parameter called softness .The ANC pair potential is given by φ anc ≡ φ anc ( r ij ; δ m , (cid:15), s ) = (cid:15) (cid:34)(cid:18) δ m ζ ij (cid:19) − (cid:18) δ m ζ ij (cid:19) (cid:35) , (1) where ζ = δ + (cid:2) r − δ (cid:3) / s , (2)and r ij denotes the relative distance between the i th and j th particle, ( δ m ) length scale, ( (cid:15) ) the depth of the poten-tial well, and non-conformality is enforced by the softnessof the potential ( s ) . Softness measures the ratio betweenof the slope and that of the reference potential, in thiscase the Lennard-Jones. The reference potential is recov-ered at s = 1 . It has been shown to provide an excellentapproximation for many gases . Furthermore, it hasan analytical expression for the second virial coefficient. As a new step, we add Gaussian functions to theANC potential. This is motivated by the fact thatstructure-based inversion methods, such as Inverse MonteCarlo, Inverse Boltzmann and Force Matching, show thatcoarse-graining from SPC/E and TIP4P atomistic watermodels yields effective potentials with two minima.
A straightforward way to implement this is to add aGaussian function to the ANC potential. The Gaussiancontribution can be given by φ g ( r ij ; h, p, q ) = h exp (cid:34) − ( r ij − p ) q (cid:35) (3)with peak height h , peak position p , and a standard de-viation q . By adding one and two Gaussian functionsto the ANC potential (Eq. 1), we obtained the ANC-1Gand ANC-2G potentials, respectively. While some pre-vious studies have employed as many as four Gaussianfunctions , we determined that using two is sufficient inorder not to introduce excess free parameters.The ANC, ANC-1G and ANC-2G are still symmetric.Another way to extend the ANC potential to capturemore molecular properties is to introduce orientationaldependence. In atomistic water models, this is usuallyachieved by placing opposite atomic partial charges onhydrogens and oxygens to induce intermolecular hydro-gen bonding, and at the coarse-grained level the Polariz-able Martini , mW and the Mercedes-Benz modelshave orientational dependence. In the latter two, the ap-proach is inspired by the models for silicon which alsoexhibit tetrahedral order. Both the mW and Mercedes-Benz model also display the density anomaly of water.The polarizable Martini model, albeit at a higher level ofcoarse-graining, uses a fluctuating dipole thus resemblingthe GSM approach here. We use the GSM potential,which combines a point dipole from the Stockmayer potential with the ANC potential (Eq. 1) without anyGaussian functions: φ gsm ( r ij , Ω ij ; δ m , (cid:15), µ ) = φ anc + φ dd ( r ij , Ω ij ; µ ) , (4)where φ dd is the dipole-dipole interaction given as φ dd ( r ij , Ω ij ; µ ) = 14 πε (cid:34) (cid:126)µ i · (cid:126)µ j r ij − (cid:126)r ij · (cid:126)µ i ) ( (cid:126)r ij · (cid:126)µ j ) r ij (cid:35) , with µ being the dipole strength, ε the permittivity ofvacuum, and Ω ij stands for the set of angles that de-fines the relative orientation between molecules i and j .In this case, the Stockmayer potential and dipolar hardspheres are obtained from Eq. (4) with s = 1 and s → ,respectively.
3. Other CG models
Since we want to assess transferrability to differentstate points, we also studied some of the most pop-ular coarse-grained water models to compare with thenew ANC and GSM potentials, namely, monatomic wa-ter (mW) water and the electrostatic-based (ELBA)model .The mW model is an adaptation of the Stillinger-Weber silicon potential that favors a tetrahedral coor-dination of molecules. The mW model was parameter-ized to reproduce the experimental melting temperatureof hexagonal ice as well as the density and enthalpy of va-porization of liquid water at ambient conditions . In themW model, each water molecule is mapped to one beadthat interacts through both a two-body and a three-bodypotential, described by E = (cid:88) i (cid:88) j>i φ ( r ij ) + (cid:88) (cid:88) (cid:88) φ ( r ij , r ik , θ ijk ) , (5)where φ ( r ij ) = Aε (cid:34) B (cid:18) σr ij (cid:19) − (cid:35) exp (cid:18) γσr ij − aσ (cid:19) and φ ( r ij , r ik , θ ijk ) = λε [cos θ ijk − cos θ ] × exp (cid:18) γσr ij − aσ (cid:19) exp (cid:18) γσr ik − aσ (cid:19) , where A = 7 . and B = 0 . . The restof the parameters are θ = 109 . ◦ , (cid:15) = 6 . kcal/mol, σ = 2 . Å, a = 1 . , λ = 23 . , and γ = 1 . . Lastly, we also included the ELBA model origi-nally developed as a solvent for lipid membranes withthe aim of reducing the computational cost. The ELBAmodel describes a water molecule as a single CG bead em-bedded with a point dipole. The potential energy for apair of CG beads is the sum of Lennard-Jones and dipoleinteractions, both terms in a shifted-force form and aregiven by φ sf − lj =4 (cid:15) (cid:40)(cid:34)(cid:18) σr ij (cid:19) − (cid:18) σr ij (cid:19) (cid:35) + (cid:34) (cid:18) σr c (cid:19) − (cid:18) σr c (cid:19) (cid:35) × (cid:18) r ij r c (cid:19) − (cid:18) σr c (cid:19) + 4 (cid:18) σr c (cid:19) (cid:41) (6)and φ sf − dd = 14 πε (cid:34) − (cid:18) r ij r c (cid:19) + 3 (cid:18) r ij r c (cid:19) (cid:35) × (cid:34) ( (cid:126)µ i · (cid:126)µ j ) r ij − (cid:126)r ij · (cid:126)µ i ) ( (cid:126)r ij · (cid:126)µ j ) r ij (cid:35) where the set of parameters are (cid:15) = 0 .
55 kcal mol − , σ = 3 . Å, µ = 2 . , and r c = 12 Å. This model wasoriginally parametrized to reproduce bulk density anddiffusion coefficient of liquid water at ◦ C and 1 atm . B. Model Optimization
The model parameters for both the ANC and GSMmodels were optimized to best reproduce the behaviorof the SPC/E reference system. Parametrizations weremade using both bulk and slab systems, using the ra-dial distribution function and the density profile, respec-tively, as the target properties. The ANC parametriza-tions, however, became unstable at the intermediate op-timization parameters for the slab geometry. The gen-eral procedure for potential optimization has been previ-ously described elsewhere and implemented in the Ver-satile Object-Oriented Toolkit for Coarse-graining Ap-plications (VOTCA) package which was used here.Here, the focus is on the selection of the potential pa-rameters during the potential update step.Due to limited sampling, results of MD simulationscontain noise, which can lead to significant errors in nu-merically determined gradients of the target properties.Therefore, gradient based function optimization meth-ods are a poor choice for optimizing model parameters.To circumvent this limitation, we used the downhill sim-plex algorithm to optimize all model parameters. It isa deterministic optimization procedure designed to mini-mize an objective (or penalty) function of n +1 variables, y ( x , x , . . . , x n +1 ) . This method requires only functionevaluations, not derivatives. The idea is to employ amoving simplex in the n -dimensional parameter spaceto surround the optimal point and then shrink the sim-plex until its dimensions reach a specified error toler-ance. A simplex is a geometrical figure consisting, in n dimensions, of n + 1 vertices x , x , . . . , x n +1 connectedby straight lines and bounded by polygonal faces. It istransformed successively using basic operations such asreflection, expansion, contraction, or reduction in orderto move towards the global minimum, i.e., to minimizethe objective function . In the context coarse-graining,the use of simplex optimization was pioneered in theMüller-Plathe group for coarse-graining of polymers .It is currently also used in speeding up optimization ofparemeters for complex biomolecular force fields such asAmber and others .To obtain a reference for the CG pair potential andto verify consistency, we also employed the IBI method with the atomistic structure of the bulk system as thetarget property. Both the simplex procedure and IBIwere performed with the VOTCA package. C. Simulation details
For the reference system, each SPC/E water moleculewas mapped onto its center of mass to represent a one-bead coarse-grained water molecule. Each SPC/E sys-tem consisted of , water molecules. A cubic boxof side l b = 40 . Å was used for the bulk system, anda rectangular box obtained by expanding the cubic boxby a factor of two in the z -direction ( l s = 2 l b ) for theslab system, generating a slab of liquid water surroundedby vapor on both sides (Fig. 1). Equations of motionswere integrated with the velocity-Verlet algorithm anda time step of . Both reference and CG simulationswere simulated in the canonical ensemble (constant num-ber of molecules N , constant volume V and constanttemperature T ) using the Nosé-Hoover thermostat with T = 300 K and a relaxation time of
200 fs , exceptfor CG liquid-gas equilibrium simulations, where a relax-ation time of
500 fs was used. In the case of the GSMmodel, the NVT/sphere integrator was used for dipolesin order to update the orientation of the dipole momentand to apply the Nosé-Hoover thermostat to its rota-tional degrees of freedom. This was necessary to pre-vent large amounts of energy from accumulating in rota-tional degrees of freedom while translational degrees offreedom experience a much smaller effective temperature.All MD simulations were performed with the LAMMPSpackage modified to permit simulations using the ANCpotential.Atomistic reference simulations were run for
10 ns ,bond lengths and angles were constrained using theSHAKE algorithm , long-range electrostatic interac-tions were handled with the particle-particle particle-mesh method , and the van der Waals interactionswere cut off at Å. For the case of point dipoles, FIG. 1: Configuraton of the slab system. Bulk water isplaced into the center and spreads out along the Z -axis,forming regions of smaller density, sometimesevaporating into the surrounding vacuum. Figureproduced with VMD .long-range interactions were handled with the ewald/dispsolver. About − simplex iterations were performed foreach of CG pair potentials studied. The number of it-erations depends on the selection of the initial set ofparameters and can not be determined a priori . Forfaster (wall time) convergence of the parameters, sam-pling times were kept short, similar to other works usingthe VOTCA package for coarse-graining . During ev-ery iteration, starting from the same initial configuration,a damped dynamics energy minimization method, calledquick-min , was performed, followed by CG simulationsfor
100 ps with the first
50 ps ignored during analysis.Radial distribution functions, density profiles, and pairpotentials (Figs. 2 and 3) were obtained from these sim-ulations.
D. Measurement of Densities in Liquid-Gas Equilibria
To measure liquid and gas densities as a function oftemperature, we performed simulations using theslab geometry (Fig. 1). The systems were divided intoequally sized sections along the z -axis for analysis andthe average density of each section ρ ( z ) was computed.A periodic function ρ fit ( z ) = ρ g + (cid:88) n = − f ( z mod ( l b ) , nl b ) (7)composed of a sum of logistic functions f ( z, nl b ) = ρ l − ρ g e − k ( z − z + nl b ) − ρ l − ρ g e − k ( z − z + nl b ) (8)was fitted to reproduce the densities from the simulationsusing the non-linear least squares curve fitting algorithmof SciPy . Each pair of logistic functions (Eq. 8) rep-resents a region of liquid water of density ρ l surroundedby water vapor of density ρ g located in a periodic cell n with cell width l b along the z-axis. The transition re-gions between the two states are characterized by a de-cay parameter k and are centered on z and z . Use of aperiodic function with contributions from multiple peri-odic images of the same liquid region (Eq. 7) is requiredto account for the possibility of a liquid-gas transitionoccurring near the periodic boundaries and to ensure asmooth curvature of the fitted function in the region. III. RESULTS AND DISCUSSION
In this work we employed two kinds of interactions: 1)isotropic, Eq. (1), 2) anisotropic via inclusion of a point-dipole interaction, Eq. (4). In both cases, parameterswere obtained from the radial distribution functions anddensity profiles. We present first the results obtainedfrom bulk systems, and then those obtained from slabgeometry. To test transferability, the parameters werealso employed in simulations over a temperature rangefrom
300 K to
550 K . A. Bulk water system
The radial distribution functions g ( r ) and the pair po-tentials u ( r ) obtained from the simplex procedure areshown in Figs. 2a and b, respectively, for both isotropicand anisotropic interactions. The IBI results, the atom-istic reference, as well as the results for mW and ELBAmodels are shown for comparison.Figure 2a shows that using the ANC potential with asingle Gaussian function (ANC-1G) yields an additionalunphysical peak in g ( r ) near . Å. We were able tocorrect this by applying a second Gaussian function tothe potential (ANC-2G). The resulting g ( r ) follows veryclosely the g ( r ) from IBI and the underlying atomisticreference (SPC/E) model. This is similar to previousfindings showing that a spherically symmetricalsingle particle potential needs to have two wells to re-cover atomistic water structure. As Fig. 2b shows, theIBI procedure produces two potential wells, where theirmaxima correspond to radii of the first and second hydra-tion shells. Both of the ANC models are purely repulsive,as the added Gaussian functions fill up the minimum ofthe base ANC potential, Eq. (1). Meanwhile the IBI effec-tive potential does have an attractive region. Therefore,it is possible that better parameters for ANC-1G andANC-2G exist but are unreachable by the simplex opti-mization procedure or require thermodynamic propertiesas targets in addition to structural ones.The g ( r ) of the GSM potential, on the other hand, ex-hibited a significantly shifted narrow second peak at . Åand a very broad and flat first minimum at . Å to . Å.The first peak, however, closely follows the reference sys-tem. These discrepancies indicate that adding a pointdipole to a single particle water model is not sufficient to make it reproduce atomistic structure. Furthermore,visual inspection of the trajectory revealed that parti-cles described by the GSM potential were not percolatingvery far.As the above shows, addition of a simple dipole mo-ment is insufficient to capture the effects of tetrahedralhydrogen bonding of water. Addition of a quadrupolemoment may be required, but that was not done in thecurrent study. The mW water model mimicks this ef-fect with a three-body term instead. With the exceptionof a slightly wider first peak, the mW water reproducesthe reference g ( r ) very well. Although it does not havea second minimum in its effective potential (but has athree-body term), the repulsive force drops significantlynear . Å, the region of the first minimum in the g ( r ) of the IBI and the reference atomistic potentials. This issimilar to the behavior of ANC-1G and ANC-2G in thesame region, Fig. 2b. B. Slab system
As an alternative means of model parametrization, wetried to fit to the density of a slab configuration at
300 K .As discussed by Ismail et al. , density profile can be usedas an order parameter in a system containing liquid-vaporinterfaces. Unlike Ismail et al. who compared surfacetensions in different atomistic water models and stud-ied capillary waves, we restrict ourselves to the densityprofile only and use it to obtain the vapor-liquid phase di-agram (the next section). However, this approach failedfor the ANC-1G and ANC-2G potentials, as during opti-mization the simulation code could not handle the inter-mediate parameters. This is appears to be due to the re-pulsive nature of those potentials. In addition, as shownby Ismail et al. , even the atomistic models underesti-mate surface tension. In that light, the failure of theANC-1G and ANC-2G potentials is not a surprise. Im-provements are a topic of a future study.As shown in Fig. 3, the density profiles of GSM, mWand the ELBA model are in good agreement with theatomistic model. As an example, for GSM the den-sity matches well the atomistic model ( . / cm vs. . / cm for SPC/E). C. Transferability
Significant deviations in liquid density were observedwhen the GSM-slab (parametrized to density) and mWpotentials were used in slab systems at higher tempera-tures ( T ≥
350 K ), Fig. 4. Water described by the GSM-slab potential evaporates much easier than the referenceand has a much lower triple point. The exact position ofthe triple point proved to be impossible to determine withcurrent methods. Meanwhile, the mW model remainedliquid and stayed at a nearly constant density regard-less of temperature. This can likely be explained by the r [ Å ] g ( r ) (a) AtomisticIBIGSMmWELBAANC-1GANC-2G r [ Å ] u [ k c a l / m o l] (b) IBIGSMmWELBAANC-1GANC-2G
FIG. 2: Radial distribution (a) and pair potential (b) ofANC-1G, ANC-2G, and GSM potentials followingSimplex optimization of the bulk water system usingthe radial distribution function as the target property.Results from iterative Boltzmann inversion (IBI)method and the atomistic reference are also shown.These are further compared to previously publishedmodels . The pair potential shown for GSM doesnot include the dipole contributions, while the pairpotential for mW is obtained from an IBI fit to the mWmodel and shows an equivalent potential lacking 3-bodyinteractions. In general all dashed pair potentials areresults of IBI.short range nature of its potential, where a molecule ingas phase is too far away from any other molecules tocontribute to the total potential energy of the system.Furthermore, when ANC-1G and ANC-2G potentialsparametrized to bulk g ( r ) were simulated in slab systems,both evaporated and spread evenly through the wholevolume at all attempted temperatures. This is due to z [ Å ] g / c m AtomisticGSMmWELBA
FIG. 3: Density z-profile of the slab system at
300 K .The center of the slab is at z = 4 . Å.the repulisive nature of the potentials as discussed in theprevious section. While they are suitable for use at thetemperature and mean density they were parametrizedat, these potentials, as well as the GSM and mW po-tentials, appear not to be transferable to other thermo-dynamic conditions. These results suggest that a trans-ferable water model needs to be parametrized not onlyagainst multiple properties (density or g ( r ) , diffusion co-efficient, tetrahedral order parameter, etc.), but possiblyalso at different state points (temperature, pressure) si-multaneously. One possibility is to interpolate betweensuch state points. Addressing this in detail is, however,beyond the current study. Density [ g/ cm ] T [ K ] AtomisticGSMmWELBA
FIG. 4: Vapor-liquid phase diagram for the GSMpotential derived from slab simulations, SPC/E ,mW and ELBA models. Points near the criticalpoint are not shown due to difficulty of finding densitiesin that region. IV. CONCLUSIONS
An ideal water model should be both fast and accu-rate. However, as we coarse-grain the system to make thesimulations faster, we also lose degrees of freedom, com-plexity of the potential, and (at least some) descriptivepower of the model. This is exemplified by the ANC-1Gand GSM parametrizations in the bulk system, wherethe functional forms of the potentials prevent accuratereproduction of g ( r ) even though g ( r ) was the targetproperty for parameter optimization. This necessitatesa more complex two-body potential, often one withoutan explicit functional form, like the one obtained fromIBI.While such potentials are unique for each statepoint , they change with thermodynamic variableslike temperature and pressure. This is the transferabilityproblem : potentials parametrized for one state pointare often poor representations of other state points. Thisis exemplified in the slab systems by the density pro-file of the of GSM and the complete evaporation of theANC-1G and ANC-2G potentials. All existing CG mod-els suffer from this to some degree and often involve acompromise between quantitatively reproducing the gen-eral trend over many states and reproducing any givenstate exactly.Furthermore, even with a two-body effective poten-tial exactly reproducing the g ( r ) at a state point, itis impossible to recapture the thermodynamic prop-erties of an underlying multi-body potential , an ef-fect known as representability. In water, the multi-body contribution arises from interactions involving non-uniform charge distribution, which are ignored in single-particle CG models. Representability problems can alsobe observed in anisotropic systems being modeled withisotropic potentials . Even though it is anisotropic, thepoint dipole introduced in the GSM is still a two-bodypotential. It is not enough to obtain accurate thermo-dynamic properties of water, so the potential needs toinclude either higher orders of the multipole expansionor a true multi-body contribution, likely one that explic-itly incorporates hydrogen bonding. Therefore, in sub-sequent works, we will expand the GSM potential to in-clude higher order moments and explore a parametriza-tion method aimed to simultaneously fit physical prop-erties at multiple state points. The best-fit parametersfound here are given below in Appendix A. ACKNOWLEDGMENTS
Appendix A: Best-fit parameters
The best-fit parameters for the ANC-1, ANC-2G, andGSM models are presented in Tables I and II.TABLE I: Values of the best-fit parameters obtained forANC-1G and ANC-2G models (Eqs. 1 and 3) with g ( r ) as the target property ANC-1G ANC-2G GSM (cid:15) [kcal / mol] δ m [ Å ] s h [kcal / mol] p [ Å ] q [ Å ] h [kcal / mol] – 0.102 – p [ Å ] – 5.53 – q [ Å ] – 0.368 – µ [D] – – 3.58 TABLE II: Values of the best-fit parameters obtainedfor GSM model (Eq. 4 ) with density z -profile as thetarget property GSM (cid:15) [kcal / mol] δ m [ Å ] s µ [D] M. Chaplin, Nat. Rev. Mol. Cell Biol. , 861 (2006). P. Ball, Nature , 291 (2008). T. D. Kühne and R. Z. Khaliullin, Nature Comm. , 1450 (2013). L. B. Skinner, C. J. Benmore, J. C. Neuefeind, and J. B. Parise,J. Chem. Phys. , 214507 (2014). J. C. Palmer, F. Martelli, Y. Liu, R. Car, A. Z. Panagiotopoulos,and P. G. Debenedetti, Nature , 385 (2014). J. T. Titantah and M. Karttunen, Soft Matter , 7977 (2015). D. Thirumalai, G. Reddy, and J. E. Straub, Acc. Chem. Res. , 83 (2012), 1107.4820. M.-C. Bellissent-Funel, A. Hassanali, M. Havenith, R. Hench-man, P. Pohl, F. Sterpone, D. van der Spoel, Y. Xu, and A. E.Garcia, Chem. Rev. , 7673 (2016). B. Guillot, J. Mol. Liq. , 219 (2002). Z. Wu, Q. Cui, and A. Yethiraj, J. Phys. Chem. B , 10524(2010). Y. Mao and Y. Zhang, Chem. Phys. Lett. , 37 (2012). S. Izadi, R. Anandakrishnan, and A. V. Onufriev, J. Phys. Chem.Lett. , 3863 (2014). C. J. Tainter, L. Shi, and J. L. Skinner, J. Chem. Theory Com-put. , 2268 (2015). I. Pethes and L. Pusztai, J. Chem. Phys. , 064506 (2017). Z. Steinczinger, P. Jóvári, and L. Pusztai, J. Mol. Liq. , 19(2017). T. Murtola, A. Bunker, I. Vattulainen, M. Deserno, and M. Kart-tunen, Phys. Chem. Chem. Phys. , 1869 (2009). S. Riniker, J. R. Allison, and W. F. van Gunsteren, Phys. Chem.Chem. Phys. , 12423 (2012). H. I. Ingólfsson, C. a. Lopez, J. J. Uusitalo, D. H. de Jong, S. M.Gopal, X. Periole, and S. J. Marrink, Wiley Interdiscip. Rev.Comput. Mol. Sci. , 225 (2013). T. Lazaridis and R. Versace, Isr. J. Chem. , 1074 (2014). K. R. Hadley and C. McCabe, Mol Simul , 671 (2012). L. Darré, M. R. Machado, and S. Pantano, Wiley Interdiscip.Rev. Comput. Mol. Sci. , 921 (2012). S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, andA. H. D. Vries, J. Phys. Chem. B , 7812 (2007). S. O. Yesylevskyy, L. V. Schäfer, D. Sengupta, and S. J. Marrink,PLoS Comput Biol , e1000810 (2010). S. V. Bennun, A. N. Dickey, C. Xing, and R. Faller, Fluid PhaseEq. , 18 (2007). C. Xing and R. Faller, J. Phys. Chem. B , 7086 (2008). M. Habibi, C. Denniston, and M. Karttunen, Europhys. Lett. , 28005 (2014). H. Wang, C. Junghans, and K. Kremer, Eur. Phys. J. E , 221(2009). F. del Río, J. E. Ramos, and I. A. McLure, J. Phys. Chem. B , 10568 (1998). V. Molinero and E. B. Moore, J. Phys. Chem. B , 4008 (2009). M. Orsi and J. W. Essex, PLoS One , e28637 (2011). M. Orsi and J. W. Essex, Faraday Discuss. , 249 (2013). H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys.Chem. , 6269 (1987). D. Reith, M. Pütz, and F. Müller-Plathe, J. Comput. Chem. ,1624 (2003). J. D. van der Waals, in KNAW, Proceedings, Vol. 15 (1913) pp.1912–1913. K. M. Watson, Ind. Eng. Chem. , 398 (1943). W.-G. Dong and J. H. Lienhard, Can. J. Chem. Eng. , 158–161(1986). M. G. Noro and D. Frenkel, J. Chem. Phys. , 2941 (2000). G. Foffi and F. Sciortino, J. Phys. Chem. B , 9702 (2007). C. Vega and J. L. F. Abascal, Phys. Chem. Chem. Phys. ,19663 (2011). J. L. F. Abascal and C. Vega, J. Chem. Phys. , 234505 (2005). I. A. McLure, J. E. Ramos, and F. del Río, J. Phys. Chem. B , 7019 (1999). F. del Río, J. E. Ramos, and I. A. McLure, Phys. Chem. Chem.Phys. , 4937 (1999). J. E. Ramos, F. del Río, and I. A. McLure, Phys. Chem. Chem.Phys. , 2731 (2000). J. E. Ramos, F. del Río, and I. A. McLure, Phys. Chem. Chem.Phys. , 2634 (2001). A. González-Calderón and A. Rocha-Ichante, J. Chem. Phys. , 034305 (2015). A. P. Lyubartsev, M. Karttunen, I. Vattulainen, and A. Laak-sonen, Soft Mater. , 121 (2003). V. Rühle, C. Junghans, A. Lukyanov, K. Kremer, and D. An-drienko, J. Chem. Theory Comput. , 3211 (2009). N. M. Barraz, E. Salcedo, and M. C. Barbosa, J. Chem. Phys. (2009), 10.1063/1.3213615. C. L. Dias, T. Ala-Nissila, M. Grant, and M. Karttunen, J.Chem. Phys. , 54505 (2009). E. Ávalos, F. del Río, and S. Lago, J. Phys. Chem. B , 508(2005). W. H. Stockmayer, J. Chem. Phys. , 398 (1941). F. H. Stillinger and T. A. Weber, Phys. Rev. B , 5262–5271(1985). M. Orsi, Mol. Phys. , 1 (2014). S. Y. Mashayak, M. N. Jochum, K. Koschke, N. R. Aluru,V. Rühle, and C. Junghans, PLoS One , e0131754 (2015). J. A. Nelder and R. Mead, Comput. J. , 308 (1965). R. Faller, H. Schmitz, O. Biermann, and F. Müller-Plathe, J.Comput. Chem. , 1009 (1999). H. Meyer, O. Biermann, R. Faller, D. Reith, and F. Müller-Plathe, J. Chem. Phys. , 6264 (2000). D. Reith, H. Meyer, and F. Müller-Plathe, Macromolecules ,2335 (2001). R. M. Betz and R. C. Walker, J. Comput. Chem. , 79 (2014). C. G. Mayne, J. Saam, K. Schulten, E. Tajkhorshid, and J. C.Gumbart, J. Comput. Chem. , 2757 (2013). W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson,J. Chem. Phys. , 637 (1982). S. Nosé, J. Chem. Phys. , 511 (1984). W. G. Hoover, Phys. Rev. A , 1695 (1985). S. Plimpton, J. Comput. Phys. , 1 (1995). W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graph. ,33 (1996). J.-P. Ryckaert, G. Ciccotti, and H. J. Berendsen, J. Comput.Phys. , 327 (1977). R. W. Hockney and J. W. Eastwood,Computer Simulation Using Particles (CRC Press, New York,1988). E. Pollock and J. Glosli, Comput. Phys. Commun. , 93 (1996). A. Toukmaji, C. Sagui, J. Board, and T. Darden, J. Chem. Phys. , 10913 (2000). D. Sheppard, R. Terrell, and G. Henkelman, J. Chem. Phys. , 134106 (2008). E. Jones, T. Oliphant, P. Peterson, et al., “SciPy: Open sourcescientific tools for Python,” (2001–), [Online; accessed 2016-11-16]. A. E. Ismail, G. S. Grest, and M. J. Stevens, J. Chem. Phys. , 014702 (2006). R. Henderson, Phys. Lett. A , 197–198 (1974). M. E. Johnson, T. Head-Gordon, and A. A. Louis, J. Chem.Phys. (2007). A. Louis, Mol. Phys.109