Bottom-up coarse-grained models describe the intermolecular structure of all-atom (AA) models with desirable accuracy and efficiency. Unfortunately, structure-based models for liquids tend to dramatically overestimate the thermodynamic pressure and, consequently, tend to vaporize under ambient conditions. By employing a volume potential to introduce additional cohesion, self-consistent pressure-matching provides a simple and robust method for accurately reproducing the pressure equation of state (EoS) for homogeneous fluids, while still preserving an accurate description of intermolecular structure. Because they depend upon the global density, though, volume potentials cannot be directly employed for inhomogeneous systems, such as liquid-vapor interfaces. In the present work, we demonstrate that volume potentials can be readily adapted as potentials of the local density. The resulting local-density potentials provide an accurate description of the structure, pressure EoS, and local density fluctuations of an AA model for liquid methanol. Moreover, we demonstrate that very slight modifications to these local-density potentials allow for a quantitative description of either local or global density fluctuations. Most importantly, we demonstrate that the resulting potentials, which were parameterized to describe a homogeneous liquid, also generate stable liquid-vapor coexistence. However, further work is necessary to more accurately reproduce the interfacial density profile.

By eliminating unnecessary atomic details, coarse-grained (CG) models facilitate investigations of phenomena that cannot be effectively studied with all-atom (AA) molecular dynamics (MD) simulations.1,2 In particular, bottom-up structure-based CG models often provide considerable efficiency gains, while still accurately modeling intermolecular structure.3–5 Unfortunately, these models often provide a poor description of thermodynamic properties and tend to demonstrate limited transferability.5 In particular, structure-based methods often generate overly repulsive potentials that dramatically overestimate the thermodynamic pressure and lack the necessary cohesion to stabilize condensed phases under ambient conditions. In practice, ad hoc constraints or corrections are often employed to ensure that the structure-based potentials generate the correct mean density at fixed external pressures.6–8 However, these ad hoc approaches can introduce discrepancies in other thermodynamic observables, such as compressibility,9 and may provide a poor description of interfacial systems.10 Consequently, several groups have sought improved transferability and thermodynamic properties by employing pair potentials that explicitly vary as a function of the local density.11–16 

Alternatively, Das and Andersen proposed a very simple and elegant approach for accurately describing the pressure of homogeneous fluids by introducing a potential, UV, that depends only upon the volume.17 Because this “volume potential” is independent of configuration, it directly impacts the thermodynamic pressure without perturbing the equilibrium configuration distribution. Subsequently, Dunn and Noid developed a self-consistent pressure-matching (PM) approach for parameterizing UV in order to quantitatively reproduce AA pressure equations of state (EoS’s).18–20 Of course, because volume potentials depend upon the global density, they cannot be directly applied for simulating liquid-vapor interfaces or other inhomogeneous systems. However, it is quite natural to consider replacing potentials of the global density with potentials of the local density. Indeed, following the work of Pagonabarraga and Frenkel,21 several recent studies have demonstrated the remarkable utility of local-density-dependent (LDD) potentials for modeling a wide range of phenomena, including lipid bilayers, shockwaves, and hydrophobic association.22–27 

In the present work, we demonstrate that the volume potential, UV, obtained from self-consistent pressure-matching can be readily adapted as a function, U ρ , of the local density. The resulting CG model quite accurately describes the pair structure, local density fluctuations, and pressure equation of state (EoS) for an AA model of liquid methanol. Moreover, very slight modifications to U ρ allow for a quantitative description of either the local density fluctuations or the pressure EoS. However, it appears that both U ρ and UV are required to simultaneously reproduce both local and global density fluctuations of homogeneous fluids. More importantly, we demonstrate that the resulting model also provides a reasonable description of liquid-vapor coexistence, although further work is required to more accurately describe the structure of the interface. Consequently, self-consistent pressure-matching may provide a particularly simple means for estimating local-density potentials, which can then be further refined via standard techniques.

Pressure-matching17,18 develops potentials of the form U PM ( R , V ) = U R ( R ) + U V ( V ) , where R is the CG configuration, V is the volume, and UR is the interaction potential, which is typically of a molecular mechanics form. Self-consistent pressure-matching determines the volume potential UV to match the simulated AA pressure EoS,18 which corresponds to minimizing a relative entropy28,29 for the constant NPT ensemble.19 In practice, we define U V ( V ) = N u ρ ( ρ 0 / ρ ¯ ) , where

u ρ ( x ) = ψ 1 x 1 + ψ 2 x 1 1 2 ,
(1)

ρ0 = N/V is the instantaneous global density, and ρ ¯ = N / V . Given this form, ψ1 and ψ2 specify corrections to the model pressure and compressibility, respectively.

In this work, we consider CG models with a potential U ( R ) = U R ( R ) + U L D D ( R ) , where

U L D D ( R ) = I U ρ ( ρ I ) ,
(2)
ρ I ρ I ( R ) = v 1 1 + J ( I ) θ ( R I J )
(3)

is the local density about site I, RIJ is the distance between sites I and J, θ (R) is an indicator function, and v is the integral of θ (R). While previous studies have emphasized the benefits of Lucy functions30,31 for modeling the local density,23,27 we adopted the functional form previously proposed by Sanyal and Shell,25 

θ ( R ) = 1 , R R 0 , c 0 + c 2 R 2 + c 4 R 4 + c 6 R 6 , R 0 R R 1 , 0 , R 1 R ,
(4)

where the constants c0, c2, c4, and c6 are explicitly defined in the supplementary material. This indicator function has a continuous first derivative. Moreover, it has finite derivatives only on the interval R 0 < R < R 1 . Importantly, the definitions for ρI and θ ensure a close relationship between the local and global densities. However, note that ρ I v 1 and, consequently, does not vanish in the gas phase.

Although U L D D is a many-body potential, it generates pair-additive forces,21 

F I , LDD = J ( I ) F ρ ( ρ I ) + F ρ ( ρ J ) r J I ,
(5)

where F ρ ( x ) = d U ρ ( x ) / d x , r J I = v 1 θ ( R I J ) R ^ I J , and R ^ I J is the unit vector pointing from site J to site I. Consequently, U L D D generates forces on site I that are due to, and in the direction of, sites J for which R 0 R I J R 1 , while the magnitude of these forces is determined by the local densities, ρI and ρJ. Similarly, U L D D generates pair-additive contributions to the instantaneous pressure,

P L D D = 1 3 V I F ρ ( ρ I ) J ( I ) v 1 θ ( R I J ) R I J .
(6)

It is quite natural to optimize U ρ by minimizing a relative entropy29,32,33 as suggested by Sanyal and Shell,25 

S r e l [ U ] = d R p ( R ) ln p ( R ) / P ( R ; U ) ,
(7)

where p ( R ) is the equilibrium distribution of the CG sites that is determined by the AA model and the CG mapping, while P ( R ; U ) is the corresponding distribution for a CG model with potential U. If we define a local density operator,34 

ρ ^ ρ ( x ) ρ ^ ρ ( R ; x ) = N 1 I δ ( ρ I ( R ) x ) ,
(8)

then it follows that

δ S r e l / δ U ρ ( x ) = β N P ρ ( x ; U ) p ρ ( x ) ,
(9)

where β = 1 / k B T , while p ρ ( x ) and P ρ ( x ; U ) are averages of ρ ^ ρ ( x ) for the equilibrium distributions p ( R ) and P ( R ; U ) , respectively. Consequently, for any given UR, S r e l is minimized with respect to U ρ when the CG model reproduces the atomic local density distribution, i.e., P ρ = p ρ . In analogy to the inverse Monte Carlo method for structure-based coarse-graining,35 one can estimate the correction, Δ U ρ , for achieving this criteria,

p ρ ( x ) P ρ ( x ; U ) = d x K ( x , x ; U ) Δ U ρ ( x ) ,
(10)
K ( x , x ; U ) = β N δ ρ ^ ρ ( x ) δ ρ ^ ρ ( x ) U ,
(11)

where δ ρ ^ ρ ( x ) = ρ ^ ρ ( x ) P ρ ( x ; U ) quantifies fluctuations in the local density and the subscripted brackets indicate averages with respect to the CG distribution P ( R ; U ) . In practice, we calculate Δ U ρ by discretizing Eq. (10) and numerically solving the resulting system of linear equations.

In Sec. IV, we report constant NPT simulations of 968 methanol molecules sampling the liquid phase at 300 K and 1 bar external pressure, as well as constant NVT simulations of 2400 methanol molecules sampling liquid-vapor coexistence at 300 K. In the present section, we briefly summarize the most important details of these simulations. The supplementary material describes our calculations in much greater detail.

All simulations employed a 1 fs time step and periodic boundary conditions in all three dimensions. The AA simulations adopted the optimized potentials for liquid simulations (OPLS)-AA potential36 to model methanol interactions, while employing the particle mesh Ewald method37 to calculate electrostatic interactions. We employed the GROMACS software package38,39 to perform AA simulations. We employed a modified version of the LAMMPS software package40 to perform CG simulations. We employed standard protocols41–48 to sample the relevant ensembles.

We consider several one-site CG models for methanol. In all models, the site corresponds to the molecular center of mass. Additionally, the CG models all employ the same interaction potential, UR, which consists of the short-ranged pair potentials that are presented in Fig. S1 of the supplementary material. We calculated UR by first mapping the AA simulations of the liquid phase to the CG representation and then employing the multiscale coarse-graining (MS-CG) force-matching49,50 variational principle.51–54 

However, the CG models differ in their use of volume or local-density potentials. Given the fixed MS-CG interaction potential, we performed self-consistent pressure-matching18,19 to determine a volume potential, UV, that accurately reproduced the AA pressure EoS. Starting from this volume potential, we then developed several different local-density potentials, U L D D , as detailed in Sec. IV.

We present results for several one-site CG models of methanol. The models all employ the same short-ranged MS-CG pair potentials. The models are distinguished by their use of global-density-dependent (GDD) potentials, i.e., volume potentials, or local-density-dependent (LDD) potentials. These potentials were all optimized for constant NPT simulations of bulk liquid methanol.

Figure 1 assesses the GDD-PM model for liquid methanol. This model employs a Global-Density-Dependent (GDD) potential, i.e., a volume potential, that is obtained from self-consistent Pressure-Matching (PM). Figure 1(a) demonstrates that this model quite accurately reproduces the site-site radial distribution function (rdf) of the AA model, although it slightly underestimates the first peak and slightly overestimates the second peak of the AA rdf. Figure 1(b) and Table I demonstrate that the GDD-PM model quantitatively reproduces the simulated pressure EoS for the AA model. If the volume potential were neglected, the CG model would overestimate the internal pressure by over 103 bars and spontaneously vaporize at 1 bar external pressure. Figure 1(c) compares the density fluctuations sampled by the AA and CG models. In Fig. 1(c) and in the following, we indicate global and local densities by ρ0 and ρ1, respectively. Because the GDD-PM model is parameterized to reproduce the AA EoS, it also quantitatively reproduces the global density fluctuations of the AA model, as indicated by the dashed curves in Fig. 1(c). The solid curves demonstrate that the GDD-PM model also quite accurately describes the local density fluctuations sampled by the AA model. However, the CG distribution is slightly broader and shifted towards slightly lower densities. In particular, the error in the mean, Δ ρ 1 , reflects the slight error in the rdf, Δg(r): Δ ρ 1 = v 1 { 1 + ρ ¯ 0 R 1 d r 4 π r 2 θ ( r ) Δ g ( r ) } = 3.6 × 1 0 5 Å−3. The tendency of the CG model to sample lower local densities presumably reflects the fact that the GDD-PM model has maximum configurational entropy relative to any other model with the same pair distribution.55 Equivalently, this discrepancy suggests that the many-body potential of mean force (PMF) reflects cohesive many-body interactions that are absent from the MS-CG pair potentials. The volume potential incorporates the effects of this many-body cohesion that are necessary for accurately modeling the pressure EoS.

FIG. 1.

Comparison of AA (black curves) and GDD-PM CG (red curves) models for bulk liquid methanol. (a), (b), and (c) present the simulated rdfs, pressure EoS’s, and density distributions, respectively. In (a), the dashed red curve presents θ (r). In (c), the solid curves present distributions of the local density, ρ1, i.e., Prob(ρ1), while the dashed curves present distributions of the global density, ρ0, which are scaled by a factor of 10 for clarity, i.e., Prob(ρ0)/10.

FIG. 1.

Comparison of AA (black curves) and GDD-PM CG (red curves) models for bulk liquid methanol. (a), (b), and (c) present the simulated rdfs, pressure EoS’s, and density distributions, respectively. In (a), the dashed red curve presents θ (r). In (c), the solid curves present distributions of the local density, ρ1, i.e., Prob(ρ1), while the dashed curves present distributions of the global density, ρ0, which are scaled by a factor of 10 for clarity, i.e., Prob(ρ0)/10.

Close modal
TABLE I.

Simulated properties for liquid methanol. We present the average global density, ρ 0 , and local density, ρ 1 , in units of 10−2 molecules/Å3, the isothermal compressibility, κ = V 1 V / P T , in units of 105 bar−1, and the standard deviation in the local density, σ 1 = δ ρ 1 2 1 / 2 , in units of 10−4 molecules/Å3.

Model ρ 0 ρ 1 κ σ1
AA  1.456  1.459  11.42  7.09 
GDD-PM  1.456  1.455  11.18  7.18 
LDD-PM  1.460  1.467  11.18  7.66 
LDD-PV  1.456  1.462  11.86  7.70 
LDD-PDF  1.455  1.461  8.25  7.14 
Model ρ 0 ρ 1 κ σ1
AA  1.456  1.459  11.42  7.09 
GDD-PM  1.456  1.455  11.18  7.18 
LDD-PM  1.460  1.467  11.18  7.66 
LDD-PV  1.456  1.462  11.86  7.70 
LDD-PDF  1.455  1.461  8.25  7.14 

Next, given the interaction potential, UR, we constructed a new model by directly converting the volume potential, UV(V), to a function of the local density, Uρ(ρ1). We refer to the resulting model as the LDD-PM model since it employs a Local-Density-Dependent (LDD) potential that is determined from pressure-matching (PM). In particular, we define U ρ ( ρ 1 ) = N 1 U V ( N / ρ 1 ) = u ρ ( ρ 1 / ρ ¯ ) , which is explicitly defined by Eq. (1) and indicated by the black curve in Fig. 2. Because uρ is a monotonically decreasing function of ρ, it generates cohesive forces at all densities. In the GDD-PM model, these forces acted only upon the simulation box and had no influence upon the configurational distribution at a given volume. However, in the LDD-PM model, these forces promote local compaction around each molecule in every configuration. It should be noted that, while UV(V) is optimized for a narrow range of global densities, U ρ ( ρ 1 ) is applied to a much wider range of local densities. We also note that U ρ as ρ 1 0 , which has little significance for liquid methanol, but will preclude liquid-vapor coexistence.

FIG. 2.

Local-density potentials, U ρ . The black curve indicates uρ obtained from self-consistent pressure-matching. The blue and green curves indicate the local-density potentials employed in the LDD-PV and LDD-PDF models. The inset compares the latter potentials to uρ by presenting Δ U ρ = U ρ u ρ over the range of local densities sampled in bulk simulations of liquid methanol.

FIG. 2.

Local-density potentials, U ρ . The black curve indicates uρ obtained from self-consistent pressure-matching. The blue and green curves indicate the local-density potentials employed in the LDD-PV and LDD-PDF models. The inset compares the latter potentials to uρ by presenting Δ U ρ = U ρ u ρ over the range of local densities sampled in bulk simulations of liquid methanol.

Close modal

We define the local densities, ρ1, according to Eq. (3) in terms of the indicator function, θ (R), which is defined in Eq. (4) and indicated by the dashed curve in Fig. 1(a). Sanyal and Shell employed the same functional form for θ (R) and explicitly studied the sensitivity of the relative entropy to R1.25 In the present work, we heuristically define the transition region over the interval 6.32 Å R 9.60 Å, which corresponds to successive minima in the rdf of Fig. 1(a). Consequently, U L D D only generates forces between molecules separated by a solvation shell. Moreover, because the transition region is relatively wide and the forces are proportional to θ , the intermolecular forces due to U L D D will be relatively weak. Accordingly, we anticipate that U L D D should minimally impact the local pair structure.

Figure 3 assesses the accuracy of the resulting LDD-PM model for simulating liquid methanol in the constant NPT ensemble. In Fig. 3, the solid black curves present results for the AA model, the dashed black curves re-capitulate results for the GDD-PM model with the potential UR + UV, and the solid red curves present results for the LDD-PM model with the potential U R + U L D D .

FIG. 3.

Assessment of CG models for liquid methanol. (a) and (b) present errors in the simulated rdf and pressure EoS. (c) compares simulated local density distributions. The vertical dotted line in (b) indicates the mean density of the AA model at 1 bar external pressure, while the inset of (c) highlights the peak of the distribution. Solid and dashed black curves present results for the AA and GDD-PM models, respectively. Red, blue, and green solid curves present results for the LDD-PM, LDD-PV, and LDD-PDF models, respectively.

FIG. 3.

Assessment of CG models for liquid methanol. (a) and (b) present errors in the simulated rdf and pressure EoS. (c) compares simulated local density distributions. The vertical dotted line in (b) indicates the mean density of the AA model at 1 bar external pressure, while the inset of (c) highlights the peak of the distribution. Solid and dashed black curves present results for the AA and GDD-PM models, respectively. Red, blue, and green solid curves present results for the LDD-PM, LDD-PV, and LDD-PDF models, respectively.

Close modal

Figure 3(a) presents the errors in the simulated rdfs. As expected, U L D D had relatively little effect upon the simulated rdf. However, because U L D D introduces additional intermolecular cohesion, its inclusion slightly increases the first peak of the rdf such that the LDD-PM model reproduces the AA rdf with slightly better accuracy than the GDD-PM model. Figure 3(b) presents errors, ΔP, in the simulated pressure EoS. The EoS for the LDD-PM model is approximately parallel to the AA EoS, but shifted to slightly lower pressures by approximately 45 bars. Thus, the LDD-PM model accurately reproduces the AA compressibility, but, at 1 bar external pressure, the LDD-PM model overestimates the AA average global density by Δ ρ 0 = 4.7 × 1 0 5 Å−3, which corresponds to an error of 0.3%. Finally, Fig. 3(c) demonstrates that the LDD-PM model reproduces the AA local density distribution slightly less accurately than the GDD-PM model. Because U ρ is attractive at all densities, the local density distribution for the LDD-PM model is broader and shifted to higher densities relative to the AA or GDD-PM models. In particular, the LDD-PM model slightly overestimates the average local density of the AA model by Δ ρ 1 = 8.3 × 1 0 5 Å 3 , which is equivalent to a 0.6% error.

Because the pressure EoS for the LDD-PM model is shifted to slightly higher densities, we empirically shifted U ρ by 4.16 × 1 0 5 Å−3 towards lower densities, which has the effect of making U L D D less attractive in all configurations. Interestingly, this shift corresponds to half of the discrepancy, Δ ρ 1 , observed in Fig. 3 for the local density of the LDD-PM model. Additionally, in order to eliminate the ρ 1 0 divergence in U ρ , we quadratically extrapolated U ρ starting from the lowest global densities sampled in simulations of liquid methanol. We present the resulting potential, U ρ , as the blue curve in Fig. 2 and refer to the corresponding model as the LDD-PV model, since it employs a LDD potential that has been constructed to match the AA pressure-volume (PV) EoS. While the quadratic extrapolation is clearly visible in Fig. 2, the shift in U ρ cannot be discerned on the scale of the main figure. The inset of Fig. 2 demonstrates that, for the vast majority of local densities that are sampled at 1 bar external pressure, this shift corresponds to reducing the local-density potential by less than 0.02 kcal/mol.

The blue curves in Fig. 3 present results for the LDD-PV model. Figure 3(a) demonstrates that the LDD-PM and LDD-PV models generate indistinguishable rdfs. Figure 3(b) and Table I demonstrate that, as desired, the LDD-PV model reproduces the simulated AA pressure EoS with nearly quantitative accuracy. Figure 3(c) demonstrates that this shift in U ρ resulted in a nearly identical shift in the simulated local-density distribution. Consequently, the peak of the LDD-PV distribution now aligns better with the AA distribution, although it remains approximately 10% too low.

We also developed a LDD-PDF model by employing linear response theory to estimate the necessary correction Δ U ρ to the LDD potential in order to accurately reproduce the AA local-density probability density function (PDF). We determined Δ U ρ by numerically solving the system of equations that are obtained by discretizing Eq. (10). Since the calculated correction was quite noisy, we fit the calculated Δ U ρ to a sum of two Gaussian functions, as detailed in the supplementary material. The resulting correction is almost imperceptible on the scale of Fig. 2. While U ρ still monotonically decreases with increasing density, this correction significantly changes the curvature of U ρ for the largest and smallest simulated local densities.

Figure 3 presents results for the LDD-PDF model. As expected, Fig. 3(a) demonstrates that these slight changes to U ρ have little impact upon the rdf. Nevertheless, these changes significantly impact the simulated pressure EoS and local-density distribution. Indeed, Fig. 3(c) demonstrates that the LDD-PDF model reproduces the AA local-density distribution with nearly quantitative accuracy. However, Fig. 3(b) indicates that, while the LDD-PDF model reproduces the AA density at 1 bar external pressure, it now underestimates the compressibility of the AA model by 28%. Consequently, it appears that, if U ρ is optimized to reproduce the AA local-density distribution, then it may be necessary to also include a volume potential to accurately model the AA pressure EoS.

Finally, we investigated the transferability of the LDD models for modeling liquid-vapor interfaces. We performed simulations of AA and CG models for 2400 methanol molecules in a periodic box with sides of 3.95 nm in the x- and y-directions and 51.0 nm in the z-direction, corresponding to an 11 nm slab of liquid methanol in coexistence with 40 nm of vapor. Figure 4 presents the results of these simulations as a function of z with z = 0 being the midplane of the liquid slab. We analyzed the interface by dividing the simulated slab into sections of width Δz = 1.0 Å and performing averages over these sections. The top and bottom panels of Fig. 4 present the simulated averages of the global density, ρ0, and local density, ρ1, as a function of z. The solid and dashed black curves in Fig. 4 present results for the AA and GDD-PM models, respectively, while the solid red and green curves present results for the LDD-PM and LDD-PDF models, respectively. The results for the LDD-PV model (not shown) are very similar, but slightly worse than the results for the LDD-PDF model.

FIG. 4.

Profiles of global (a) and local density (b) for methanol liquid-vapor coexistence. The solid and dashed black curves present results for the AA and GDD-PM models, while the solid red and green curves present results for the LDD-PM and LDD-PDF models. The inset compares the interfacial global density profiles on a logarithmic scale in order to highlight the coexisting vapor phase.

FIG. 4.

Profiles of global (a) and local density (b) for methanol liquid-vapor coexistence. The solid and dashed black curves present results for the AA and GDD-PM models, while the solid red and green curves present results for the LDD-PM and LDD-PDF models. The inset compares the interfacial global density profiles on a logarithmic scale in order to highlight the coexisting vapor phase.

Close modal

Figure 4 demonstrates that the AA model generated stable liquid-vapor coexistence. As previously discussed, the GDD-PM model vaporizes and does not stabilize a liquid phase. Indeed, as noted by Wagner et al., this is to be expected for bottom-up models that are naïvely parameterized to reproduce the structure of homogeneous fluids without consideration of thermodynamic properties.27 In contrast, each LDD model does stabilize the liquid phase, although Table II indicates that the simulated liquid densities are slightly greater than the AA liquid density. Because U ρ as ρ 1 0 for the LDD-PM model, this model cannot stabilize a vapor phase. However, the LDD-PDF model, for which U ρ more reasonably describes the vapor phase, demonstrates stable liquid-vapor coexistence. The inset of Fig. 4(a) indicates that the vapor phase density is approximately an order of magnitude too small in the LDD-PDF model. This can presumably be remedied by properly tuning the ρ 1 0 limit of U ρ .

TABLE II.

Simulated properties of the liquid/vapor interface. We present the average global density for the bulk liquid phase, ρ 0 , and the bulk vapor phase, ρ 0 v , in units of molecules/Å3.

Model ρ 0 ρ 0 v
AA  1.449 × 1 0 2   9.0 × 1 0 6  
GDD-PM  …  3.0 × 1 0 3  
LDD-PM  1.460 × 1 0 2   … 
LDD-PDF  1.455 × 1 0 2   5.7 × 1 0 7  
Model ρ 0 ρ 0 v
AA  1.449 × 1 0 2   9.0 × 1 0 6  
GDD-PM  …  3.0 × 1 0 3  
LDD-PM  1.460 × 1 0 2   … 
LDD-PDF  1.455 × 1 0 2   5.7 × 1 0 7  

Figure 4 also demonstrates that the LDD models do not accurately reproduce the interfacial profile of the AA model. The atomistic global density profile transitions smoothly from the liquid to the vapor phase over an interval approximately 15 Å wide. Similarly, the AA local density profile smoothly transitions from the liquid to the vapor phase over a slightly wider interval of approximately 25 Å. In contrast, the LDD models accumulate at the liquid-vapor interface and demonstrate too sharp an interfacial profile. These deficiencies are most striking for the LDD-PM model but are also apparent for the LDD-PDF model. Interestingly, the discrepancies are more pronounced in the global density profile. Consequently, further work is required to more accurately model liquid-vapor interfacial profiles.

Motivated by the work of Das and Andersen,17 we have recently developed a self-consistent pressure-matching method that optimizes volume potentials for accurately reproducing target pressure EoS’s.18,19 However, volume potentials cannot stabilize coexistence between vapor and condensed phases. In the present work, we demonstrate that the volume potentials obtained from self-consistent pressure-matching can be readily adapted as local-density potentials. The resulting models provide a very accurate description of liquid methanol and can also be used to simulate liquid-vapor coexistence.

Very recently, the Shell and Voth groups have reported quite similar investigations of bottom-up CG models with local-density potentials.25,27 Sanyal and Shell simultaneously optimized both pair and local-density potentials by minimizing the relative entropy with respect to AA models for solvated hydrophobes.25 The resulting potentials determined accurate implicit solvent models for hydrophobic polymers and solutes with considerable transferability for longer polymers and varying solute concentrations. Conversely, Voth and co-workers optimized pair and local-density potentials by minimizing the MS-CG force-matching functional with respect to AA models for liquid-vapor interfaces.27 The MS-CG approach provided an exceptionally accurate model for methanol liquid-vapor interfaces, which we also consider here. Interestingly, although MS-CG models parameterized for bulk acetonitrile provided a very accurate model for the bulk structure, the MS-CG approach appeared less successful for acetonitrile interfaces.

The present approach combines aspects of both prior studies. In particular, we employed the MS-CG approach of Voth and co-workers to parameterize the interaction potential, UR, because the MS-CG approach is simple and direct (i.e., noniterative). Moreover, the MS-CG potential quite accurately reproduces the intermolecular packing of liquid methanol. Conversely, we employed the functional form previously adopted by Shell and co-workers for the indicator function, θ (R), since it provides a simple relation between the local and global densities and, thus, facilitates the conversion of UV to U L D D . However, motivated by simple physical intuition, we employed a relatively broad and long-ranged transition region, such that U L D D minimally impacted the intermolecular packing.

Our approach also reflects important differences. The Shell and Voth groups both employed a single variational optimization to simultaneously determine both the interaction potential and the local-density potential. In contrast, we optimized the interaction and local-density potentials separately. We employed the MS-CG variational principle to optimize the pair potentials for approximating the configuration-dependence of the PMF but employed a relative entropy variational principle to optimize UV for approximating the volume-dependence of the PMF. Perhaps even more importantly, while both prior works employed inhomogeneous environments to optimize the CG models, we first parameterized the CG model for a homogeneous liquid and then examined its transferability for inhomogeneous environments.

The present approach demonstrates several clear disadvantages. First of all, because we parameterized the CG model for the liquid phase, we obtained a relatively inaccurate description of the interfacial profile. Additionally, the separate parameterization of interaction and local-density potentials is clearly less elegant than a single self-consistent variational calculation and may potentially introduce inconsistencies between the various contributions to the total potential. These inconsistencies do not emerge with volume potentials, since UV does not alter the configuration distribution at fixed volumes but could easily arise when pair potentials are supplemented by independently optimized local-density potentials. As noted above, we minimized these inconsistencies by defining θ such that U ρ should not alter intermolecular packing. More generally, it may prove fruitful to formalize the construction of CG models in terms of mutually commuting operators corresponding to observables that are orthogonal in some sense, as suggested by Voth and co-workers.56 

An additional significant drawback of the present approach is that the volume potential is optimized over a relatively narrow range of densities, while U ρ needs to be applied over a much wider range of local densities. We speculate that this may be the largest source of discrepancies observed in both bulk and interfacial simulations. In particular, it is notable that the small corrections introduced in the LDD-PDF model primarily correspond to global densities that are rarely sampled by liquid methanol. Similarly, the errors in the interfacial profiles primarily correspond to global densities in the two-phase coexistence region that are not sampled by liquid methanol.

Nevertheless, the present approach may prove quite useful. In particular, given any interaction potential, self-consistent pressure-matching provides a particularly simple and robust method for determining a volume potential that provides the necessary cohesion to stabilize condensed phases. The present work suggests that it may be generally straightforward to convert this volume potential into a local-density potential. Moreover, by focusing on reproducing the structure and thermodynamic properties of bulk phases, the present approach avoids the possibility of introducing finite size effects that can significantly impact thermodynamic properties.57 

Moreover, notwithstanding the aforementioned limitations, it may possibly prove fruitful to separately consider structural and thermodynamic properties when constructing bottom-up CG models. Structure-based bottom-up approaches appear exceptionally useful for approximating the configuration-dependence of the many-body potential of mean force (PMF), especially the configuration-dependence governing local structure and intermolecular packing. However, thermodynamic properties may be quite sensitive to contributions in the PMF that are configuration-independent or that vary only weakly with configuration.20,58 The key insight of pressure-matching is that one can easily parameterize volume potentials, UV, to accurately model the volume-dependence of the PMF and, consequently, the pressure EoS.17,18 Notably, pressure-matching relies upon information that is not easily extracted from the configuration-dependence of the PMF and is most usefully applied in simulations with an external pressure.

Similarly, in parameterizing local-density potentials, it may prove beneficial to more explicitly consider the thermodynamics of phase coexistence. For instance, if one defines

Z 0 = V d R exp β U R ( R )
(12)

as the bulk-phase configuration integral for the MS-CG interaction potential, UR, then the corresponding configuration integral for LDD models may be approximated,

Z LDD = Z 0 exp β U L D D ( R ) 0 Z 0 exp β N U ρ ( ρ * ) ,
(13)

where the subscripted angular brackets denote an average with respect to the original MS-CG model (or the GDD-PM model) and ρ* indicates a characteristic local density.59,60 Consequently, U ρ should directly contribute to excess thermodynamic properties and, in particular, to the excess chemical potentials.21 While it may be convenient to employ pressure-matching to estimate U ρ ( ρ ) for a single phase, it may be quite simple to tune U ρ to accurately reproduce phase coexistence.

In closing, we note that local-density potentials appear to be a promising avenue not only for extending the transferability of bottom-up CG models but also for modeling internal degrees of freedom,23,61–64 for developing more accurate treatments of thermodynamic properties, such as phase coexistence,10,65 and even for connecting with multi-resolution66 and ultra coarse-graining67 schemes. However, CG models are sensitive to even minute changes in LDD potentials. Consequently, it is quite encouraging that self-consistent pressure-matching appears to provide a simple and robust means for obtaining first estimates of U ρ . Perhaps surprisingly, in combination with the MS-CG interaction potential, the local-density potential obtained from pressure-matching provides an improved description of the AA rdf, as well as a fairly accurate description of the AA pressure EoS and local-density fluctuations. Moreover, linear response theory then determines a local-density potential that almost quantitatively reproduces the AA local-density fluctuations, although a volume potential appears necessary to simultaneously reproduce the AA pressure EoS. Most importantly, these local-density potentials allow for a reasonable description of liquid-vapor coexistence, though further work is clearly necessary to provide a quantitative description of the interfacial profile. Thus, we conclude that this work may provide useful insights for developing more rigorous bottom-up approaches for modeling phase coexistence and inhomogeneous systems.

See supplementary material for detailed information regarding AA and CG simulations, the force- and pressure-matching calculations, the parameters for the indicator function θ , and also the numerical solution of Eq. (10).

The authors gratefully acknowledge financial support from the National Science Foundation (Grant No. CHE-1565631). Portions of this research were conducted with Advanced CyberInfrastructure computational resources provided by The Institute for CyberScience at The Pennsylvania State University (http://ics.psu.edu).

1.
C.
Peter
and
K.
Kremer
, “
Multiscale simulation of soft matter systems
,”
Faraday Discuss.
144
,
9
24
(
2010
).
2.
S.
Riniker
,
J. R.
Allison
, and
W. F.
van Gunsteren
, “
On developing coarse-grained models for biomolecular simulation: A review
,”
Phys. Chem. Chem. Phys.
14
,
12423
12430
(
2012
).
3.
E.
Brini
,
E. A.
Algaer
,
P.
Ganguly
,
C.
Li
,
F.
Rodriguez-Ropero
, and
N. F. A.
van der Vegt
, “
Systematic coarse-graining methods for soft matter simulations: A review
,”
Soft Matter
9
,
2108
2119
(
2013
).
4.
M. G.
Saunders
and
G. A.
Voth
, “
Coarse-graining methods for computational biology
,”
Annu. Rev. Biophys.
42
,
73
93
(
2013
).
5.
W. G.
Noid
, “
Perspective: Coarse-grained models for biomolecular systems
,”
J. Chem. Phys.
139
,
090901
(
2013
).
6.
D.
Reith
,
M.
Putz
, and
F.
Muller-Plathe
, “
Deriving effective mesoscale potentials from atomistic simulations
,”
J. Comput. Chem.
24
,
1624
1636
(
2003
).
7.
S.
Jain
,
S.
Garde
, and
S. K.
Kumar
, “
Do inverse Monte Carlo algorithms yield thermodynamically consistent interaction potentials?
,”
Ind. Eng. Chem. Res.
45
,
5614
5618
(
2006
).
8.
T.
Murtola
,
E.
Falck
,
M.
Karttunen
, and
I.
Vattulainen
, “
Coarse-grained model for phospholipid/cholesterol bilayer employing inverse Monte Carlo with thermodynamic constraints
,”
J. Chem. Phys.
126
,
075101
(
2007
).
9.
H.
Wang
,
C.
Junghans
, and
K.
Kremer
, “
Comparative atomistic and coarse-grained study of water: What do we lose by coarse-graining?
,”
Eur. Phys. J. E
28
,
221
229
(
2009
).
10.
M.
Jochum
,
D.
Andrienko
,
K.
Kremer
, and
C.
Peter
, “
Structure-based coarse-graining in liquid slabs
,”
J. Chem. Phys.
137
,
064102
(
2012
).
11.
E. C.
Allen
and
G. C.
Rutledge
, “
A novel algorithm for creating coarse-grained, density dependent implicit solvent models
,”
J. Chem. Phys.
128
,
154115
(
2008
).
12.
E. C.
Allen
and
G. C.
Rutledge
, “
Evaluating the transferability of coarse-grained, density-dependent implicit solvent models to mixtures and chains
,”
J. Chem. Phys.
130
,
034904
(
2009
).
13.
S.
Izvekov
,
P. W.
Chung
, and
B. M.
Rice
, “
The multiscale coarse-graining method: Assessing its accuracy and introducing density dependent coarse-grain potentials
,”
J. Chem. Phys.
133
,
064109
(
2010
).
14.
S.
Izvekov
,
P. W.
Chung
, and
B. M.
Rice
, “
Particle-based multiscale coarse graining with density-dependent potentials: Application to molecular crystals (hexahydro-1,3,5-trinitro-s-triazine)
,”
J. Chem. Phys.
135
,
044112
(
2011
).
15.
A. J.
Clark
,
J.
McCarty
,
I. Y.
Lyubimov
, and
M. G.
Guenza
, “
Thermodynamic consistency in variable-level coarse graining of polymeric liquids
,”
Phys. Rev. Lett.
109
,
168301
(
2012
).
16.
J.
McCarty
,
A. J.
Clark
,
J.
Copperman
, and
M. G.
Guenza
, “
An analytical coarse-graining method which preserves the free energy, structural correlations, and thermodynamic state of polymer melts from the atomistic to the mesoscale
,”
J. Chem. Phys.
140
,
204913
(
2014
).
17.
A.
Das
and
H. C.
Andersen
, “
The multiscale coarse-graining method. V. Isothermal-isobaric ensemble
,”
J. Chem. Phys.
132
,
164106
(
2010
).
18.
N. J. H.
Dunn
and
W. G.
Noid
, “
Bottom-up coarse-grained models that accurately describe the structure, pressure, and compressibility of molecular liquids
,”
J. Chem. Phys.
143
,
243148
(
2015
).
19.
N. J. H.
Dunn
and
W. G.
Noid
, “
Bottom-up coarse-grained models with predictive accuracy and transferability for both structural and thermodynamic properties of heptane-toluene mixtures
,”
J. Chem. Phys.
144
,
204124
(
2016
).
20.
N. J. H.
Dunn
,
T. T.
Foley
, and
W. G.
Noid
, “
Van der Waals perspective on coarse-graining: Progress toward solving representability and transferability problems
,”
Acc. Chem. Res.
49
,
2832
2840
(
2016
).
21.
I.
Pagonabarraga
and
D.
Frenkel
, “
Dissipative particle dynamics for interacting systems
,”
J. Chem. Phys.
115
,
5015
5026
(
2001
).
22.
M.
Homberg
and
M.
Muller
, “
Main phase transition in lipid bilayers: Phase coexistence and line tension in a soft, solvent-free, coarse-grained model
,”
J. Chem. Phys.
132
,
155104
(
2010
).
23.
J. D.
Moore
,
B. C.
Barnes
,
S.
Izvekov
,
M.
Lísal
,
M. S.
Sellers
,
D. E.
Taylor
, and
J. K.
Brennan
, “
A coarse-grain force field for RDX: Density dependent and energy conserving
,”
J. Chem. Phys.
144
,
104501
(
2016
).
24.
V.
Agrawal
,
P.
Peralta
,
Y.
Li
, and
J.
Oswald
, “
A pressure-transferable coarse-grained potential for modeling the shock Hugoniot of polyethylene
,”
J. Chem. Phys.
145
,
104903
(
2016
).
25.
T.
Sanyal
and
M. S.
Shell
, “
Coarse-grained models using local-density potentials optimized with the relative entropy: Application to implicit solvation
,”
J. Chem. Phys.
145
,
034109
(
2016
).
26.
J. F.
Dama
,
J.
Jin
, and
G. A.
Voth
, “
The theory of ultra-coarse-graining. 3. Coarse-grained sites with rapid local equilibrium of internal states
,”
J. Chem. Theory Comput.
13
,
1010
1022
(
2017
).
27.
J. W.
Wagner
,
T.
Dannenhoffer-Lafage
,
J.
Jin
, and
G. A.
Voth
, “
Extending the range and physical accuracy of coarse-grained models: Order parameter dependent interactions
,”
J. Chem. Phys.
147
,
044113
(
2017
).
28.
S.
Kullback
and
R. A.
Leibler
, “
On information and sufficiency
,”
Ann. Math. Stat.
22
,
79
86
(
1951
).
29.
M. S.
Shell
, “
The relative entropy is fundamental to multiscale and inverse thermodynamic problems
,”
J. Chem. Phys.
129
,
144108
(
2008
).
30.
L. B.
Lucy
, “
A numerical approach to the testing of the fission hypothesis
,”
Astron. J.
82
,
1013
1024
(
1977
).
31.
R. A.
Gingold
and
J. J.
Monaghan
, “
Smoothed particle hydrodynamics: Theory and application to non-spherical stars
,”
Mon. Not. R. Astron. Soc.
181
,
375
389
(
1977
).
32.
A.
Chaimovich
and
M. S.
Shell
, “
Coarse-graining errors and numerical optimization using a relative entropy framework
,”
J. Chem. Phys.
134
,
094112
(
2011
).
33.
J. F.
Rudzinski
and
W. G.
Noid
, “
Coarse-graining entropy, forces, and structures
,”
J. Chem. Phys.
135
,
214101
(
2011
).
34.
W. G.
Noid
, “
Systematic methods for structurally consistent coarse-grained models
,”
Methods Mol. Biol.
924
,
487
531
(
2013
).
35.
A. P.
Lyubartsev
and
A.
Laaksonen
, “
Calculation of effective interaction potentials from radial distribution functions: A reverse Monte Carlo approach
,”
Phys. Rev. E
52
,
3730
3737
(
1995
).
36.
W. L.
Jorgensen
,
D. S.
Maxwell
, and
J.
Tirado-Rives
, “
Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids
,”
J. Am. Chem. Soc.
118
,
11225
11236
(
1996
).
37.
T.
Darden
,
D.
York
, and
L.
Pedersen
, “
Particle mesh Ewald: An N log(N) method for Ewald sums in large systems
,”
J. Chem. Phys.
98
,
10089
10092
(
1993
).
38.
S.
Pronk
,
S.
Páll
,
R.
Schulz
,
P.
Larsson
,
P.
Bjelkmar
,
R.
Apostolov
,
M. R.
Shirts
,
J. C.
Smith
,
P. M.
Kasson
,
D.
van der Spoel
,
B.
Hess
, and
E.
Lindahl
, “
Gromacs 4.5: A high-throughput and highly parallel open source molecular simulation toolkit
,”
Bioinformatics
29
,
845
854
(
2013
).
39.
M. J.
Abraham
,
T.
Murtola
,
R.
Schulz
,
S.
Páll
,
J. C.
Smith
,
B.
Hess
, and
E.
Lindahl
, “
Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers
,”
SoftwareX
1-2
,
19
25
(
2015
).
40.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular-dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
41.
S.
Nose
, “
A molecular dynamics method for simulations in the canonical ensemble
,”
Mol. Phys.
52
,
255
268
(
1984
).
42.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
,
1695
1697
(
1985
).
43.
H. J. C.
Berendsen
,
J. P. M.
Postma
,
W. F.
van Gunsteren
,
A.
DiNola
, and
J. R.
Haak
, “
Molecular dynamics with coupling to an external bath
,”
J. Chem. Phys.
81
,
3684
3690
(
1984
).
44.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
, “
Canonical sampling through velocity rescaling
,”
J. Chem. Phys.
126
,
014101
(
2007
).
45.
M.
Parrinello
and
A.
Rahman
, “
Crystal structure and pair potentials: A molecular-dynamics study
,”
Phys. Rev. Lett.
45
,
1196
1199
(
1980
).
46.
G. J.
Martyna
,
M. L.
Klein
, and
M.
Tuckerman
, “
Nosé–Hoover chains: The canonical ensemble via continuous dynamics
,”
J. Chem. Phys.
97
,
2635
2643
(
1992
).
47.
G. J.
Martyna
,
D. J.
Tobias
, and
M. L.
Klein
, “
Constant pressure molecular dynamics algorithms
,”
J. Chem. Phys.
101
,
4177
4189
(
1994
).
48.
G. J.
Martyna
,
M. E.
Tuckerman
,
D. J.
Tobias
, and
M. L.
Klein
, “
Explicit reversible integrators for extended systems dynamics
,”
Mol. Phys.
87
,
1117
1157
(
1996
).
49.
F.
Ercolessi
and
J. B.
Adams
, “
Interatomic potentials from first-principles calculations: The force-matching method
,”
Europhys. Lett.
26
,
583
(
1994
).
50.
A. J.
Chorin
, “
Conditional expectations and renormalization
,”
Multiscale Model. Simul.
1
,
105
118
(
2003
).
51.
S.
Izvekov
and
G. A.
Voth
, “
A multiscale coarse-graining method for biomolecular systems
,”
J. Phys. Chem. B
109
,
2469
2473
(
2005
).
52.
S.
Izvekov
and
G. A.
Voth
, “
Multiscale coarse graining of liquid-state systems
,”
J. Chem. Phys.
123
,
134105
(
2005
).
53.
W. G.
Noid
,
J.-W.
Chu
,
G. S.
Ayton
,
V.
Krishna
,
S.
Izvekov
,
G. A.
Voth
,
A.
Das
, and
H. C.
Andersen
, “
The multiscale coarse-graining method. I. A rigorous bridge between atomistic and coarse-grained models
,”
J. Chem. Phys.
128
,
244114
(
2008
).
54.
W. G.
Noid
,
P.
Liu
,
Y. T.
Wang
,
J.-W.
Chu
,
G. S.
Ayton
,
S.
Izvekov
,
H. C.
Andersen
, and
G. A.
Voth
, “
The multiscale coarse-graining method. II. Numerical implementation for molecular coarse-grained models
,”
J. Chem. Phys.
128
,
244115
(
2008
).
55.
A. P.
Lyubartsev
and
A.
Laaksonen
, “
Osmotic and activity coefficients from effective potentials for hydrated ions
,”
Phys. Rev. E
55
,
5689
5696
(
1997
).
56.
J. W.
Wagner
,
J. F.
Dama
,
A. E. P.
Durumeric
, and
G. A.
Voth
, “
On the representability problem and the physical meaning of coarse-grained models
,”
J. Chem. Phys.
145
,
044108
(
2016
).
57.
A.
Lyubartsev
,
A.
Mirzoev
,
L. J.
Chen
, and
A.
Laaksonen
, “
Systematic coarse-graining of molecular models by the Newton inversion method
,”
Faraday Discuss.
144
,
43
56
(
2010
).
58.
C. N.
Likos
, “
Effective interactions in soft condensed matter physics
,”
Phys. Rep.
348
,
267
439
(
2001
).
59.
J. G.
Kirkwood
, “
Statistical mechanics of fluid mixtures
,”
J. Chem. Phys.
3
,
300
313
(
1935
).
60.
B.
Widom
, “
Some topics in the theory of fluids
,”
J. Chem. Phys.
39
,
2808
2812
(
1963
).
61.
P.
Espanol
, “
Dissipative particle dynamics with energy conservation
,”
Europhys. Lett.
40
,
631
(
1997
).
62.
K.-H.
Lin
,
B. L.
Holian
,
T. C.
Germann
, and
A.
Strachan
, “
Mesodynamics with implicit degrees of freedom
,”
J. Chem. Phys.
141
,
064107
(
2014
).
63.
T. T.
Foley
,
M. S.
Shell
, and
W. G.
Noid
, “
The impact of resolution upon entropy and information in coarse-grained models
,”
J. Chem. Phys.
143
,
243104
(
2015
).
64.
G.
Faure
,
R.
Delgado-Buscalioni
, and
P.
Español
, “
The entropy of a complex molecule
,”
J. Chem. Phys.
146
,
224106
(
2017
).
65.
B.
Mukherjee
,
L.
Delle Site
,
K.
Kremer
, and
C.
Peter
, “
Derivation of coarse grained models for multiscale simulation of liquid crystalline phase transitions
,”
J. Phys. Chem. B
116
,
8474
8484
(
2012
).
66.
R.
Potestio
,
S.
Fritsch
,
P.
Español
,
R.
Delgado-Buscalioni
,
K.
Kremer
,
R.
Everaers
, and
D.
Donadio
, “
Hamiltonian adaptive resolution simulation for molecular liquids
,”
Phys. Rev. Lett.
110
,
108301
(
2013
).
67.
J. F.
Dama
,
A. V.
Sinitskiy
,
M.
McCullagh
,
J.
Weare
,
B.
Roux
,
A. R.
Dinner
, and
G. A.
Voth
, “
The theory of ultra-coarse-graining. 1. General principles
,”
J. Chem. Theory Comput.
9
,
2466
2480
(
2013
).

Supplementary Material