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.

## I. INTRODUCTION

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, *U*_{V}, 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 *U*_{V} 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, *U*_{V}, obtained from self-consistent pressure-matching can be readily adapted as a function, $ U \rho $, 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 \rho $ allow for a quantitative description of either the local density fluctuations or the pressure EoS. However, it appears that both $ U \rho $ and *U*_{V} 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.

## II. THEORY

Pressure-matching^{17,18} develops potentials of the form $ U PM ( R , V ) \u2009=\u2009 U R ( R ) \u2009+\u2009 U V ( V ) $, where $R$ is the CG configuration, *V* is the volume, and *U*_{R} is the interaction potential, which is typically of a molecular mechanics form. Self-consistent pressure-matching determines the volume potential *U*_{V} to match the simulated AA pressure EoS,^{18} which corresponds to minimizing a relative entropy^{28,29} for the constant NPT ensemble.^{19} In practice, we define $ U V ( V ) =N u \rho ( \rho 0 / \rho \xaf ) $, where

*ρ*_{0} = *N*/*V* is the instantaneous global density, and $ \rho \xaf =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

is the local density about site *I*, *R*_{IJ} is the distance between sites *I* and *J*, $\theta $(*R*) is an indicator function, and $v$ is the integral of $\theta $(*R*). While previous studies have emphasized the benefits of Lucy functions^{30,31} for modeling the local density,^{23,27} we adopted the functional form previously proposed by Sanyal and Shell,^{25}

where the constants *c*_{0}, *c*_{2}, *c*_{4}, and *c*_{6} 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 $\theta $ ensure a close relationship between the local and global densities. However, note that $ \rho I \u2265 v \u2212 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}

where $ F \rho ( x ) =\u2212d U \rho ( x ) /dx$, $ r J I = v \u2212 1 \theta \u2032 ( 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 \u2264 R I J \u2264 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,

It is quite natural to optimize $ U \rho $ by minimizing a relative entropy^{29,32,33} as suggested by Sanyal and Shell,^{25}

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}

then it follows that

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

where $\delta \rho ^ \rho ( x ) = \rho ^ \rho ( x ) \u2212 P \rho ( 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 $\Delta U \rho $ by discretizing Eq. (10) and numerically solving the resulting system of linear equations.

## III. METHODS

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 potential^{36} to model methanol interactions, while employing the particle mesh Ewald method^{37} to calculate electrostatic interactions. We employed the GROMACS software package^{38,39} to perform AA simulations. We employed a modified version of the LAMMPS software package^{40} to perform CG simulations. We employed standard protocols^{41–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, *U*_{R}, which consists of the short-ranged pair potentials that are presented in Fig. S1 of the supplementary material. We calculated *U*_{R} by first mapping the AA simulations of the liquid phase to the CG representation and then employing the multiscale coarse-graining (MS-CG) force-matching^{49,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-matching^{18,19} to determine a volume potential, *U*_{V}, 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.

## IV. RESULTS

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.

### A. GDD model for 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 10^{3} 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, $\Delta \rho 1 $, reflects the slight error in the rdf, Δ*g*(*r*): $\Delta \u27e8 \rho 1 \u27e9 = v \u2212 1 { 1 + \rho \xaf \u222b 0 R 1 d r 4 \pi r 2 \theta ( r ) \Delta g ( r ) } =\u22123.6\xd71 0 \u2212 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.

Model . | $ \rho 0 $ . | $ \rho 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 . | $ \rho 0 $ . | $ \rho 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 |

### B. LDD models for liquid methanol

Next, given the interaction potential, *U*_{R}, we constructed a new model by directly converting the volume potential, *U*_{V}(*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 \rho ( \rho 1 ) = N \u2212 1 U V ( N / \rho 1 ) = u \rho ( \rho 1 / \rho \xaf ) $, 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 *U*_{V}(*V*) is optimized for a narrow range of global densities, $ U \rho ( \rho 1 ) $ is applied to a much wider range of local densities. We also note that $ U \rho \u2192\u221e$ as $ \rho 1 \u21920$, which has little significance for liquid methanol, but will preclude liquid-vapor coexistence.

We define the local densities, *ρ*_{1}, according to Eq. (3) in terms of the indicator function, $\theta $(*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 $\theta $(*R*) and explicitly studied the sensitivity of the relative entropy to *R*_{1}.^{25} In the present work, we heuristically define the transition region over the interval 6.32 Å $\u2264R\u22649.60\u2009$ Å, 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 $ \theta \u2032 $, 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 *U*_{R} + *U*_{V}, and the solid red curves present results for the LDD-PM model with the potential $ U R + U L D D $.

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 $\Delta \rho 0 =4.7\xd71 0 \u2212 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 \rho $ 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 $\Delta \rho 1 =8.3\xd71 0 \u2212 5 \u2009 \xc5 \u2212 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 \rho $ by $4.16\xd71 0 \u2212 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, $\Delta \rho 1 $, observed in Fig. 3 for the local density of the LDD-PM model. Additionally, in order to eliminate the $ \rho 1 \u21920$ divergence in $ U \rho $, we quadratically extrapolated $ U \rho $ starting from the lowest global densities sampled in simulations of liquid methanol. We present the resulting potential, $ U \rho $, 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 \rho $ 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 \rho $ 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 $\Delta U \rho $ to the LDD potential in order to accurately reproduce the AA local-density probability density function (PDF). We determined $\Delta U \rho $ 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 $\Delta U \rho $ 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 \rho $ still monotonically decreases with increasing density, this correction significantly changes the curvature of $ U \rho $ 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 \rho $ 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 \rho $ 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.

### C. Transferability of LDD models to interfaces

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.

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 \rho \u2192\u221e$ as $ \rho 1 \u21920$ for the LDD-PM model, this model cannot stabilize a vapor phase. However, the LDD-PDF model, for which $ U \rho $ 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 $ \rho 1 \u21920$ limit of $ U \rho $.

Model . | $ \rho 0 \u2113 $ . | $ \rho 0 v $ . |
---|---|---|

AA | 1.449 $\xd7\u20091 0 \u2212 2 $ | $9.0\xd71 0 \u2212 6 $ |

GDD-PM | … | $3.0\xd71 0 \u2212 3 $ |

LDD-PM | 1.460 $\xd7\u20091 0 \u2212 2 $ | … |

LDD-PDF | 1.455 $\xd7\u20091 0 \u2212 2 $ | $5.7\xd71 0 \u2212 7 $ |

Model . | $ \rho 0 \u2113 $ . | $ \rho 0 v $ . |
---|---|---|

AA | 1.449 $\xd7\u20091 0 \u2212 2 $ | $9.0\xd71 0 \u2212 6 $ |

GDD-PM | … | $3.0\xd71 0 \u2212 3 $ |

LDD-PM | 1.460 $\xd7\u20091 0 \u2212 2 $ | … |

LDD-PDF | 1.455 $\xd7\u20091 0 \u2212 2 $ | $5.7\xd71 0 \u2212 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.

## V. DISCUSSION

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, *U*_{R}, 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, $\theta $(*R*), since it provides a simple relation between the local and global densities and, thus, facilitates the conversion of *U*_{V} 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 *U*_{V} 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 *U*_{V} 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 $\theta $ such that $ U \rho $ 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 \rho $ 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, *U*_{V}, 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

as the bulk-phase configuration integral for the MS-CG interaction potential, *U*_{R}, then the corresponding configuration integral for LDD models may be approximated,

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 \rho $ 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 \rho ( \rho ) $ for a single phase, it may be quite simple to tune $ U \rho $ to accurately reproduce phase coexistence.

## VI. CONCLUSIONS

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-resolution^{66} and ultra coarse-graining^{67} 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 \rho $. 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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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).