Coarse-grained (CG) models provide superior computational efficiency for simulating soft materials. Unfortunately, CG models with conventional pair-additive potentials demonstrate limited transferability between bulk and interfacial environments. Recently, a growing number of CG models have supplemented these pair potentials with one-body potentials of the local density (LD) around each site. These LD potentials can significantly improve the accuracy and transferability of CG models. Nevertheless, it remains challenging to accurately describe interfaces where the LD varies rapidly. In this work, we consider a new class of one-body potentials that depend upon the square of the LD gradient around each site. We investigate the impact of this square gradient (SG) potential upon both top-down dissipative particle dynamics (DPD) models and also bottom-up multiscale coarse-graining (MS-CG) models. We demonstrate that SG potentials can be used to tune the interfacial properties of DPD models without significantly altering their bulk properties. Moreover, we demonstrate that SG potentials can improve the bulk pressure–density equation of state as well as the interfacial profile of MS-CG models for acetic acid. Consequently, SG potentials may provide a useful connection between particle-based top-down models and mean-field Landau theories for phase behavior. Furthermore, SG potentials may prove useful for improving the accuracy and transferability of bottom-up CG models for interfaces and other inhomogeneous systems with significant density gradients.

## I. INTRODUCTION

By eliminating unnecessary atomic details, coarse-grained (CG) models provide the necessary efficiency for simulating phenomena that cannot be effectively studied with conventional all-atom (AA) models.^{1–9} Given a specified AA model, an “exact” CG model would describe interactions with the many-body potential of mean force (PMF), *W*, which may be defined as^{10–13}

Here, *β* = 1/*k*_{B}*T*, *u* is the AA potential, **r** is the configuration for *n* atoms, **R** is the configuration for *N* CG sites, and **M**: **r** → **R** = **M**(**r**) is a mapping that determines the CG representation of each AA configuration. In principle, the PMF can be employed to reproduce all structural and thermodynamic properties of the AA model that can be observed at the CG resolution.^{14}

In practice, though, CG models often describe intermolecular interactions with simple pair-additive potentials. For instance, structure-based bottom-up approaches often determine pair potentials that accurately reproduce the pair structure of an AA model at a single state point.^{15,16} Because *W* is a many-body function, such pair-additive approximations provide limited accuracy for describing cooperative or correlated interactions.^{17–21} Similarly, because *W* depends upon the thermodynamic state point,^{14} pair-additive approximations provide unpredictable transferability^{22–36} and often fail to reproduce thermodynamic properties, such as the pressure–density equation of state (EoS).^{37–46} While configuration-independent volume potentials can be introduced to reproduce the EoS for homogeneous systems,^{47–49} they do not help with modeling interfaces or other inhomogeneous systems.^{50–53}

Consequently, a growing number of CG models have supplemented pair-additive potentials with more complex terms to describe many-body effects.^{20,21,34,35,54–64} In particular, many recent studies have adopted one-body potentials that depend upon the local density (LD) around each site. Because the local density is defined by pair-additive contributions, these local density (LD) potentials provide similar efficiency to pair-additive potentials.^{65} These LD potentials were first introduced into top-down dissipative particle dynamics (DPD) models^{66} by Paragonabarraga and Frenkel (PF).^{65,67,68} Motivated by considerations from classical density functional theory (DFT),^{69} PF interpreted one-body LD potentials in terms of an excess free energy per particle. While traditional DPD models employed repulsive pair potentials that resulted in a quadratic equation of state,^{70} LD potentials allow many-body DPD (MDPD) models to describe much more complex phase behavior.^{71–76} Independently, Wolynes and co-workers employed similar LD potentials to model water-mediated many-body interactions in CG models for protein folding and structure prediction.^{77,78}

Subsequently, several groups have employed bottom-up approaches to parameterize LD potentials for modeling a range of complex systems, including liquid mixtures,^{79–82} polymers,^{83–85} propagating shock waves,^{86,87} and various interfaces.^{88–91} These LD potentials have improved the transferability and bulk thermodynamic properties of bottom-up models.^{84–87} However, they do not always accurately describe the structure of liquid interfaces.^{81,88,90} Classical DFT and mean-field thermodynamic theories indicate that an accurate description of liquid interfaces should account for the free energy cost associated with variations in the local density.^{69,92} Accordingly, in this work, we consider CG models with effective potentials that depend upon both the local density and also the (square of the) local density gradient around each site.

The remainder of this paper is organized as follows: Sec. II defines these square gradient (SG) potentials and analyzes their basic properties. In particular, we demonstrate that the resulting forces are pair additive, satisfy Newton’s third law, and conserve momentum, although they are not radial. Section III considers the impact of SG potentials upon top-down DPD models. We consider several functional forms for the SG potential and examine their influence upon both bulk and interfacial properties. Section IV employs the bottom-up multiscale coarse-graining (MS-CG) approach to parameterize the SG potential. We demonstrate that the SG potential improves both the bulk pressure–density equation of state and also the liquid–vapor interfacial profile of the MS-CG model. Finally, Sec. V provides concluding remarks and suggests directions for future work.

## II. SQUARE GRADIENT POTENTIALS

We consider potentials with the form

The first term is a sum of central pair potentials, *U*_{2}, that depend upon the distance, *R*_{IJ}, between each pair of sites,

The second term in Eq. (2) is a sum of local density (LD) potentials, *U*_{ρ}, that depend upon the local density, *ρ*_{I}, around each site *I*,

This local density is defined as a sum of pair-additive contributions from surrounding sites,

where $w\u0304(r)$ is a normalized weighting function that decays with distance, *r*, and vanishes at a finite cutoff. Here, we have employed the notation of PF: *w*(*r*) is an un-normalized weighting function, [*w*] = ∫d*r*4*πr*^{2}*w*(*r*) is the normalizing volume, and $w\u0304(r)=w(r)/[w]$ is the normalized weighting function.^{65} Note that including the *J* = *I* term in Eq. (5) introduces a self-density contribution, $w\u0304(0)=[w]\u22121$, that is independent of configuration and only shifts the definition of the local density. This term may impact top-down models that assume a functional form for *U*_{ρ}^{65,72} but should have no effect upon bottom-up models that parameterize flexible forms for *U*_{ρ} from AA simulation data.

The present study focuses on the third term in Eq. (2), which we refer to as a square gradient (SG) potential,

We emphasize that *U*_{∇} is not a potential function but rather a density-dependent coefficient that determines the magnitude of the SG potential. This SG potential depends upon both *ρ*_{I} and also the gradient in the local density around site *I*,

where ∇_{I} = *∂*/*∂***R**_{I} and $R\u0302IJ=RI\u2212RJ/RIJ$ is the unit vector pointing from site *J* toward site *I*. Because *R*_{IJ}, *ρ*_{I}, and $\u2207I\rho I2$ are all scalars, *U* is translationally and rotationally invariant.

Because *ρ*_{I} and ∇_{I}*ρ*_{I} are both defined by pair-additive contributions, the net force on each site *I* is also pair additive,

The total pair force may be expressed as

The first contribution in Eq. (9) comes from *U*_{pair}: *F*_{2}(*R*) = −d*U*_{2}(*R*)/d*R*. The second contribution comes from *U*_{LD} and acts along $R\u0302IJ$ with a magnitude,

where *F*_{ρ}(*ρ*) = −d*U*_{ρ}(*ρ*)/d*ρ*. Thus, the pair and LD potentials both exert equal and opposite forces on each pair of particles, *I* and *J*, that act along the interparticle vector.

The third term in Eq. (9) comes from *U*_{SG} and is more complex,

where we have defined *F*_{∇}(*ρ*) = −d*U*_{∇}(*ρ*)/d*ρ*, **A**_{I} = ∇_{I}*ρ*_{I}, $f(R)=w\u0304\u2032(R)/R$, and $g(R)=w\u0304\u2033(R)\u2212f(R)$. Because each line of Eq. (11) is antisymmetric with respect to exchanging *I* and *J*, *U*_{SG} generates equal and opposite forces between each pair of particles: **F**_{∇;JI} = −**F**_{∇;IJ}. However, the pair force, **F**_{∇;IJ}, from *U*_{SG} does not act along the interparticle vector, $R\u0302IJ$. While the first two rows of Eq. (11) act along $R\u0302IJ$, the third line does not.

Because **F**_{∇;JI} = −**F**_{∇;IJ}, Eq. (9) satisfies Newton’s third law,

which ensures that the forces conserve linear momentum. Additionally, Eq. (8) implies that the excess contribution to the pressure tensor retains a pair-additive form

where $e\u0302i$ indicates the unit vector in the *i*th Cartesian direction.^{93,94}

Because *U* is rotationally invariant, the forces must also conserve angular momentum even though the pair forces are not radial. Appendix A explicitly demonstrates that *U* generates a vanishing net torque on any system of particles.

Finally, we re-iterate that, despite the complexity of Eq. (11), the force on each site can be calculated from local pair-additive contributions. Consequently, SG potentials should provide similar computational scaling to simple pair-additive potentials, albeit with a larger prefactor.

## III. DPD MODELS

### A. DPD methodology

DPD models^{66} generally assume that particles follow Newton’s equations of motion

and that the net force on each site, *I*, may be decomposed into pair-additive contributions,

The first term, $FIJC$, arises from the gradient of a conservative potential, *U*. The second and third terms in Eq. (15) describe dissipative and random forces, respectively,

Here, **V**_{IJ} = **V**_{I} − **V**_{J} is the velocity of site *I* with respect to site *J*, *θ*_{IJ} is a pairwise white noise satisfying Gaussian statistics: $\theta IJ(t)=0$ and $\theta IJ(t)\theta KL(t\u2032)=\delta IK\delta JL+\delta IL\delta JK\delta (t\u2212t\u2032)$, and

is a quadratic pairwise weighting function expressed in terms of the Heaviside function, Θ, such that the pair force vanishes between particles with *R*_{IJ} > *R*_{c;2}. According to the fluctuation–dissipation relation, the viscous and random forces are related by *σ*^{2} = 2*γk*_{B}*T*.^{95}

MDPD models^{66} adopt a conservative potential, *U* = *U*_{pair} + *U*_{LD}, which corresponds to Eq. (2) with *U*_{SG} = 0. The MDPD potential often adopts a particularly simple form.^{71–75} Specifically, the pair potential is assumed proportional to *w*_{2},

Moreover, the LD potential is assumed quadratic in the local density,

where the local density

is defined in terms of a second quadratic weighting function

that acts on a shorter length scale, *R*_{c;ρ} = 0.75*R*_{c;2}, than *w*_{2}. As before, $w\u0304\rho (r)=w\u0304\rho (r)/[w\rho ]$ and [*w*_{ρ}] = ∫d*r*4*πr*^{2}*w*_{ρ}(*r*). This definition of the local density [i.e., Eq. (21)] corresponds to Eq. (5). Given Eqs. (19)–(22), the conservative MDPD pair force may be expressed as

which corresponds to Eq. (9) with **F**_{∇;IJ}(**R**) = **0**. The parameters *A* and *B* can then be related to the model thermodynamic properties.^{73}

In this work, we consider the influence of the SG potential, *U*_{SG}, upon MDPD models. In particular, we consider three simple forms for *U*_{∇},

In the first case, we assume that the SG potential acts equivalently upon each particle irrespective of its environment. In the second and third cases, we assume that the SG potential only acts upon particles, *I*, near the interface for which *ρ*_{I} is less than some cutoff, *ρ*_{c}. In these latter two cases, *U*_{∇} is assumed to vary linearly or quadratically with the local density of particles at the interface.

### B. Simulation parameters

DPD studies often adopt dimensionless quantities. We define the length scale by the pair potential cutoff, *R*_{c;2} = 1, the energy scale by thermal energy, *k*_{B}*T* = 1, and the time scale by requiring that all particles have a mass *m* = 1. Following prior studies,^{73–75} we define the local density cutoff *R*_{c;ρ} = 0.75, the dissipative force parameter *γ* = 4.5, and the random force parameter *σ* = 3.0 to satisfy the fluctuation–dissipation relation.^{95}

We simulated all MDPD models with a version of Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)^{96} that we modified to support both LD and SG potentials. The majority of MDPD simulations considered 2400 particles in a 5 × 5 × 40 rectangular prism with periodic boundary conditions in all three dimensions. We constructed the initial configuration by randomly placing these particles in the middle 5 × 5 × 16 region of the simulation domain. We equilibrated the system for 2 × 10^{6} time steps with *δt* = 0.001. We then simulated the system for 20 × 10^{6} additional time steps with *δt* = 0.01 while sampling configurations every 500 time steps.

In order to characterize the bulk fluid in the constant NPT ensemble, we performed several additional simulations of 1296 particles in a periodic cubic box with sides of length *L* ≈ 6. These simulations were performed in the constant NPH ensemble, since the DPD random forces provide the necessary thermostatting. We employed the MTTK barostat to apply an external pressure of 1.0 with a coupling constant of 1000.0.^{97,98} Otherwise, these simulations employed equivalent conditions and parameters to the simulations described above.

The goal of the present section is to characterize the influence of the SG potential, *U*_{SG}, upon the behavior of MDPD models. Accordingly, we define a reference baseline model with *U*_{SG} = 0 such that the conservative pair force is given by Eq. (23). In particular, we choose *A* = −40.0 and *B* = 25.0 for this reference model, since prior studies have investigated this parameter combination.^{73,75} We performed seven sets of simulations in order to investigate the influence of the potential parameters upon the properties of the model. Table I summarizes the simulated parameter sets.

Set . | Fixed param. . | Var. param. . | Sample set . |
---|---|---|---|

A | B = 25 | A | −30, −35, −37.5, −42.5, −45, −50 |

B | A = −40 | B | 15, 20, 22.5, 27.5, 30, 35 |

c | A = −40, B = 25 | c | −0.03, −0.02, −0.01, 0.01, 0.032, 0.1 |

L1.5 | A = −40, B = 25, ρ_{c} = 1.5 | m | ±0.01, ±0.02, ±0.05, ±0.1, ±0.2 |

L2.5 | A = −40, B = 25, ρ_{c} = 2.5 | m | ±0.01, ±0.02, ±0.05, ±0.1, ±0.2 |

Q1.2 | A = −40, B = 25, ρ_{c} = 1.2 | k | ±0.1, ±0.5 |

Q1.5 | A = −40, B = 25, ρ_{c} = 1.5 | k | ±0.01, ±0.1, ±0.5 |

Set . | Fixed param. . | Var. param. . | Sample set . |
---|---|---|---|

A | B = 25 | A | −30, −35, −37.5, −42.5, −45, −50 |

B | A = −40 | B | 15, 20, 22.5, 27.5, 30, 35 |

c | A = −40, B = 25 | c | −0.03, −0.02, −0.01, 0.01, 0.032, 0.1 |

L1.5 | A = −40, B = 25, ρ_{c} = 1.5 | m | ±0.01, ±0.02, ±0.05, ±0.1, ±0.2 |

L2.5 | A = −40, B = 25, ρ_{c} = 2.5 | m | ±0.01, ±0.02, ±0.05, ±0.1, ±0.2 |

Q1.2 | A = −40, B = 25, ρ_{c} = 1.2 | k | ±0.1, ±0.5 |

Q1.5 | A = −40, B = 25, ρ_{c} = 1.5 | k | ±0.01, ±0.1, ±0.5 |

The first two sets of simulations, which we refer to as sets A and B, adopted a conventional MDPD potential with *U* = *U*_{pair} + *U*_{LD} and *U*_{SG} = 0. These sets of simulations investigated the influence of the *A* and *B* parameters, following the work of Warren.^{73} In particular, set A consisted of six simulations with *B* = 25.0, while *A* was taken from the set {−30.0, −35.0, −37.5, −42.5, −45.0, −50.0}. Conversely, set B consisted of six simulations with *A* = −40.0, while *B* was taken from the set {15.0, 20.0, 22.5, 27.5, 30.0, 35.0}.

The remaining five simulation sets adopted *A* = −40.0 and *B* = 25.0 from the baseline reference model while employing various functional forms, *U*_{∇}(*ρ*_{I}), to define *U*_{SG}. Set C consisted of six simulations adopting the constant form, *U*_{∇}(*ρ*_{I}) = *c*, while *c* was taken from the set {−0.03, −0.02, −0.01, 0.01, 0.032, 0.1}. Sets L1.5 and L2.5 both consisted of ten simulations adopting the linear (L) form, *U*_{∇}(*ρ*_{I}) = *m*(*ρ*_{c} − *ρ*_{I}) Θ (*ρ*_{c} − *ρ*_{I}), while *m* was taken from the set {±0.01, ±0.02, ±0.05, ±0.10, ±0.20}. The L1.5 and L2.5 sets differed in the choice of the LD cutoff, *ρ*_{c}. The L1.5 set adopted *ρ*_{c} = 1.5, while the L2.5 set adopted *ρ*_{c} = 2.5. Finally, the last two simulation sets Q1.2 and Q1.5 adopted the quadratic (Q) form $U\u2207(\rho I)=k\rho c\u2212\rho I2\Theta \rho c\u2212\rho I$ while adopting the cutoffs *ρ*_{c} = 1.2 and *ρ*_{c} = 1.5, respectively. The set Q1.2 consisted of four simulations with *k* taken from the set {±0.1, ±0.5}, while the set Q1.5 consisted of six simulations with *k* taken from the set {±0.01, ±0.1, ±0.5}.

### C. Results

We first consider a conventional MDPD model as a baseline for considering the impact of SG potentials. This baseline model employs a standard DPD pair potential, $U2(RIJ)=12ARc;2w2(RIJ)$, where *w*_{2}(*R*) is a quadratic weighting function that vanishes at a cutoff *R*_{c;2} = 1. The baseline model also employs a standard MDPD LD potential, $U\rho (\rho I)=\pi 30BRc;\rho 4\rho I2$, where the local density, *ρ*_{I}, around particle *I* is defined by a quadratic weighting function, *w*_{ρ}(*R*_{IJ}), that vanishes at a cutoff *R*_{c;ρ} = 0.75. In particular, the baseline reference model adopts *A* = −40.0 and *B* = 25.0, since prior studies have considered these parameters.^{73,75} Since we define distance with respect to the pair cutoff, *R*_{c;2}, all densities are reported in units of $Rc;2\u22123$ and ∇_{I}*ρ*_{I} is reported in units of $Rc;2\u22124$. Since we define energy with respect to *k*_{B}*T*, pressures are reported in units of $kBT/Rc;23$, while surface tensions are reported in units of $kBT/Rc;22$.

Figure 1 presents results from constant pressure simulations of this baseline model. Panel (a) presents the pair, *w*_{2}, and LD weighting functions, *w*_{ρ}, which vary quadratically with distance, *r*, out to their respective cutoffs, *R*_{c;2} = 1 and *R*_{c;ρ} = 0.75. Panel (b) presents the simulated radial distribution function (rdf), which indicates that the pair structure of the baseline model is similar to a simple fluid. The rdf vanishes for small *r*, peaks at *r* ≈ 0.6 and then oscillates on a similar length scale until *r* ≈ 3. Since the rdf oscillates on the length scale of the weighting function, *w*_{ρ}, that is used to define the local density, the local density is likely sensitive to the local solvation structure. Panel (c) presents the simulated pressure–density equation of state as a function of the global density, *ρ* = *n*/*V*. Because the small system samples relatively large density fluctuations and the liquid is rather incompressible, the simulation samples metastable liquid states with negative internal pressure. However, the fixed external pressure, *P*_{ext} = 1, precludes the system from vaporizing. Finally, panel (d) presents the simulated distribution of global, *ρ*, and local, *ρ*_{1}, densities. While the global density distribution is sharply peaked around *ρ* ≈ 6.1, the local density distribution is much broader and peaked around the much smaller density *ρ*_{1} ≈ 2.6.

Incidentally, we observed that the narrow global density distribution includes two very close peaks that reflect oscillations in the simulated volume. While we employed the MTTK barostat^{97,98} to sample volume fluctuations, Trofimov *et al.* previously observed similar oscillations in constant pressure MDPD simulations^{99} with the Andersen barostat.^{100} However, since our primary interest is the influence of square gradient potentials upon interfacial properties, we did not investigate these oscillations further.

Figure 2 characterizes a constant NVT slab simulation of this baseline model. Panel (a) presents the average global, $\rho \u0304(z)$, and local, $\rho \u03041(z)$, densities as a function of the slab z-coordinate. The slab simulation generates coexistence between a liquid phase with a density $\rho \u0304liq\u22486$ and a vapor phase with a density $\rho \u0304vap\u22485\xd710\u22125$. The average local density of the liquid phase is $\rho \u03041;liq\u22482.5$. A rather sharp liquid–vapor interface forms near *z* ≈ ±7. Interestingly, the global density, $\rho \u0304(z)$, slightly oscillates before the interface and then rapidly plummets at the interface. In contrast, the local density, $\rho \u03041(z)$, does not oscillate but transitions rather smoothly over a considerably wider region.

Figure 2(b) presents an intensity plot for the joint distribution of z-coordinates and local densities, *ρ*_{1}, sampled by particles, i.e., *P*(*z*, *ρ*_{1}). In the bulk liquid region, the particles sample a broad distribution of local densities similar to the distribution in Fig. 1(d). In particular, note that *ρ*_{1} ≳ 1.5 in the bulk region. Near the interface, the average local density rapidly drops and the distribution sharpens.

Figure 2(c) presents an analogous intensity plot for the joint distribution of z-coordinates and local density gradients, $\u2207\rho 1$, sampled by particles, i.e., *P*(*z*, |∇*ρ*_{1}|). In the bulk liquid region, this distribution is peaked near |∇*ρ*_{1}| ≈ 3.5 and extends over a modest range of gradients. In the interfacial region, this distribution shifts to significantly larger density gradients, |∇*ρ*_{1}| ≈ 10, and also sharpens.

It is perhaps surprising that the LD gradients do not vanish in the bulk liquid region. These finite LD gradients reflect the short-ranged weighting function, *w*_{ρ}, used in defining the local density. As noted in discussing Fig. 1, this weighting function is sensitive to the many-body correlations characterizing the solvation structure.^{17} The typical LD gradients in the bulk region systematically decrease as the local density is defined over longer length scales.

Having characterized the baseline reference model, we next investigate the effects of introducing SG potentials into the MDPD models. Accordingly, Fig. 3 presents the density profiles obtained from constant NVT slab simulations for the seven families of MDPD models described in Sec. III B. Each panel in Fig. 3 presents density profiles for a single family in which a corresponding parameter was systematically varied. Figure 4 then presents the simulated bulk liquid density, *ρ*_{liq}, bulk vapor density, *ρ*_{vap}, and the liquid–vapor interfacial tension, *γ*, obtained from the corresponding simulations. The bulk liquid density was calculated from the global density profiles for −3.5 < *z* < 3.5, while the bulk vapor density was calculated from the global density profiles for |*z*| > 12.5. Throughout this work, we calculate the liquid–vapor surface tension of simulated slabs according to

where *z* indicates the direction normal to the slab.^{53} Each of the seven columns in Fig. 4 corresponds to one of the simulation families and the x-axis indicates the simulated parameter.

For reference, the first two model families illustrate the effects of varying the *A* and *B* parameters in the conventional MDPD potential. Recall that the reference baseline MDPD model corresponds to *A* = −40.0 and *B* = 25.0. The results of this baseline model are indicated by the black curves in the top two panels of Fig. 3 and by the dotted green lines in Fig. 4.

The top left panel of Fig. 3 and the first column of Fig. 4 present results for family A, in which the parameter *B* = 25.0 was fixed, while the parameter *A* varies. As *A* decreases from −30 to −50, the pair potential becomes increasingly attractive. Accordingly, the liquid density increases, the vapor density decreases over 4 orders of magnitude, the interfacial profile sharpens, and the interfacial tension increases by a factor of ∼3. The top middle panel of Fig. 3 and the second column of Fig. 4 present results for family B, in which the parameter *A* = −40.0 was fixed, while the parameter *B* varies. As *B* increases from 15 to 35, the LD potential becomes increasingly repulsive, i.e., increasingly penalizes finite local densities. Accordingly, the liquid density decreases, the vapor density increases over 2 orders of magnitude, the interfacial profile appears to broaden, and the interfacial tension decreases by a factor of ∼2.

The remaining five model families adopt the fixed parameters *A* = −40.0 and *B* = 25.0 for the conventional pair and LD potentials while adopting various functions, *U*_{∇}, for the SG potential. In particular, the top right panel of Fig. 3 and the third column of Fig. 4 present results for family c, in which *U*_{∇}(*ρ*_{I}) was assumed to be a constant, *c*, that we varied from −0.03 to 0.1. In this case, *c* < 0 favors the formation of LD gradients, while *c* > 0 penalizes these gradients. As *c* decreases from 0 to −0.03, the liquid density decreases by ∼10%, the interface broadens, and the interfacial tension decreases by a factor of roughly 3. As *c* increases from 0 to 0.1, the liquid density increases by ∼7%, and the surface tension increases by a factor of ∼2. As might be expected, the model appears more sensitive to *c* < 0 than to *c* > 0. In comparison to families A and B, the vapor density varies relatively little with the *c* parameter. More generally, we find that the vapor density is quite insensitive to the SG potentials that preserve liquid–vapor coexistence.

The remaining model families adopt SG potentials that primarily influence particles at the liquid–vapor interface. These models adopt a LD cutoff, *ρ*_{c}, such that the SG potential acts only on particles with *ρ*_{I} < *ρ*_{c}. According to Fig. 2, particles in the bulk liquid phase very rarely sample local densities below 1.5. Consequently, when *ρ*_{c} ≤ 1.5, the SG potential will very rarely act on particles in the bulk liquid phase and should have minimal influence upon the bulk thermodynamic properties. Conversely, when *ρ*_{c} > 1.5, one expects that the SG potential will exert more influence upon the bulk liquid phase.

The L1.5 and L2.5 families adopt a linear form for the SG potential with cutoffs *ρ*_{c} = 1.5 and 2.5, respectively. Specifically, *U*_{∇}(*ρ*_{I}) = *m*(*ρ*_{c} − *ρ*_{I}) Θ (*ρ*_{c} − *ρ*_{I}) such that *U*_{∇} varies linearly with local density for *ρ*_{I} < *ρ*_{c}. Given the fixed cutoff, *ρ*_{c}, *U*_{∇} depends upon a single parameter, *m*, which we varied between −0.2 and +0.2. For *m* < 0, the SG potential stabilizes local density gradients, while for *m* > 0, the SG potential penalizes local density gradients.

We present results for the L1.5 and L2.5 families in the first and second bottom panels of Fig. 3, respectively, as well as the fourth and fifth columns of Fig. 4. In both families, the bulk liquid density and interfacial tension systematically decrease with decreasing *m*; until at *m* = −0.2, the liquid–vapor interface vanishes and a single homogeneous fluid is observed. As expected, the SG potential has a greater effect upon the interfacial tension than the bulk liquid density. Moreover, since it acts on both bulk and interfacial particles particles, the L2.5 SG potential exerts considerably greater influence upon the MDPD model than the L1.5 SG potential. As observed for family c, the MDPD model appears more sensitive to parameters that stabilize local gradients and reduce the surface tension.

The third and fourth bottom panels of Fig. 3 and the last two columns of Fig. 4 present results for the Q1.2 and Q1.5 models, which employ cutoffs *ρ*_{c} = 1.2 and 1.5, respectively. The Q1.2 and Q1.5 models adopt a quadratic form for the SG potential, $U\u2207(\rho I)=k\rho c\u2212\rho I2\Theta \rho c\u2212\rho I$, such that *U*_{∇}(*ρ*_{I}) vanishes for *ρ*_{I} > *ρ*_{c} and varies quadratically with the local density for *ρ*_{I} < *ρ*_{c}. In this case, *k* < 0 stabilizes local density gradients, while *k* > 0 penalizes local density gradients. As *k* increases from −0.1 to +0.5, the liquid density and the liquid–vapor interfacial profiles appear almost constant. Over this parameter range, the interfacial tension increases by ∼5% and 20% in the Q1.2 and Q1.5 models, respectively. In contrast, for *k* = −0.5, the slab vaporizes and the system becomes a homogeneous fluid. Again, we see that the MDPD models are more sensitive to parameters that stabilize local density gradients.

Figure 4 demonstrates that SG potentials can significantly reduce the liquid–vapor surface tension. In certain parameter regimes, these SG potentials may result in a negative surface tension and drive the formation of novel phases with periodic density waves or a proliferation of interfaces. We did not observe such exotic phase behavior for the simulated parameters. Instead, in the simulated cases that the liquid slab vaporized, the system formed a homogeneous low density fluid consisting primarily of monomers in equilibrium with relatively small isotropic droplets. However, it is possible that in pathological parameter regimes, the SG potential may lead to unphysical thermodynamic instabilities.^{101,102}

According to Figs. 3 and 4, the SG potentials exert minimal impact upon the bulk liquid density in the L1.5, Q1.2, and Q1.5 models that preserve liquid–vapor coexistence. This suggests that in these cases, the SG potential can be introduced to modify the interfacial properties without altering the thermodynamic properties of the bulk liquid phase. Figure 5 presents simulated pressure–density equations of state for the Q1.2 and Q1.5 model families. These equations of state are almost independent of the SG potential for the Q1.2 model family. Thus, the Q1.2 family allows one to systematically vary the interfacial tension over a range of ∼5% without altering the bulk liquid phase. Conversely, in the case of the Q1.5 family, the equations of state shift to slightly greater densities and become very slightly steeper as *k* increases. Therefore, the Q1.5 family allows one to tune the interfacial tension by ∼20% with only slight changes to the bulk liquid density and compressibility.

## IV. BOTTOM-UP MODELS

### A. MS-CG methodology

In comparison to DPD and other top-down approaches, bottom-up approaches often employ more complex functions for describing the interactions between CG sites. Bottom-up potentials can usually be expressed by the general form

where *U*_{ζ}(*x*) is an interaction potential governing a scalar variable, *ψ*_{ζλ}(**R**), that depends upon the coordinates for the *λ* set of particles. Bottom-up approaches determine these interaction potentials based upon information from AA simulations.^{15} In particular, the MS-CG approach^{103–105} determines the potential that optimally approximates the many-body PMF by minimizing the force-matching functional,^{106,107}

where **f**_{I}(**r**) is the net force on CG site *I* in the AA configuration **r**.

The MS-CG formalism considers the CG force field to be a vector, $F={FI(R)}I=1,\u2026,N$, in an abstract vector space with **F**_{I}(**R**) = −∇_{I}*U*(**R**).^{15,106–108} In the MS-CG approach, the approximate potential in Eq. (26) specifies a basis set for a subspace of force fields. Each interaction potential, *U*_{ζ}(*x*), determines a force function, *F*_{ζ}(*x*) = −d*U*_{ζ}(*x*)/d*x*, that is represented with simple basis functions, *f*_{ζd}(*x*), and constant coefficients, *ϕ*_{ζd},

The interaction potential *U*_{ζ} determines a force on each site, *I*,

where

and $G\zeta d={G\zeta d;I(R)}I=1,\u2026,N$ defines a force field basis vector. The set of basis vectors, ${G\zeta d}$, specified by Eq. (26) then defines a basis for a subspace of CG force fields. The MS-CG variational principle determines the optimal force field in this subspace by solving a linear least squares problem for the coefficients, $\varphi \zeta d*$, that minimize *χ*^{2}.^{109}

Returning to Eq. (2), both *U*_{pair} and *U*_{LD} can be expressed in the form of Eq. (26). In particular, the present bottom-up models defined the local density, *ρ*_{I}, according to Eq. (5) and employed the Lucy function

as the un-normalized indicator function. However, the SG potential

cannot be expressed in the form of Eq. (26), and the corresponding forces cannot be expressed in the form of Eq. (29). Specifically, Eq. (11) demonstrates that these forces depend explicitly upon both *U*_{∇} and *F*_{∇}, which complicates the construction of the corresponding force field basis vectors. We emphasize again that *U*_{∇} is not a potential function but rather a density-dependent coefficient.

To construct the force field basis vectors for *U*_{SG}, we first define

According to Eq. (11), **F**_{∇;IJ} may be decomposed as

where

The function *U*_{∇}(*ρ*) is represented by simple basis functions, *f*_{∇d}(*ρ*), with constant coefficients, *φ*_{∇d},

which determines a corresponding representation of *F*_{∇}(*ρ*),

Given this representation for *U*_{∇} and *F*_{∇}, **F**_{∇;IJ} may be expressed as

where

Finally, this allows for a basis set expansion of the forces on each site *I* due to *U*_{SG},

where

defines the force field basis vector, $G\u2207d={G\u2207d;I(R)}I=1,\u2026,N$. Given this basis expansion, *χ*^{2} again becomes a quadratic function of the force field coefficients and the MS-CG potential can be determined by solving a linear least squares problem.

### B. Computational details

#### 1. All-atom simulations

In the following, we consider two all-atom (AA) simulations of acetic acid. Both simulations were performed with Gromacs 5.1.4^{110} while integrating the equations of motion with a 1 fs time step. Both simulations employed the v-rescale thermostat with a coupling constant of 0.5 ps to maintain a temperature of 300 K.^{111} Both simulations considered only neutral acetic acid molecules and modeled all interactions with the OPLS/AA force field.^{112} Dispersion interactions were truncated at 1.3 nm while using the force-switch option to transition the corresponding forces smoothly to 0 over the distance range 1.0 ≤ *r* ≤ 1.3 nm. The short-ranged component of electrostatic interactions was truncated at 1.3 nm, while the long-ranged component was treated with the particle mesh Ewald (PME) method with a grid spacing of 0.08 nm.^{113} No dispersion corrections were included in either AA simulation.

The first AA simulation modeled a homogeneous liquid in the constant NPT ensemble at a temperature of 300 K and an external pressure of 1.0 bar. This simulation considered 674 molecules in a periodic cubic box that was initially prepared with sides of *L* = 4 nm. The Parrinello–Rahman barostat was employed to apply an isotropic external pressure of 1.0 bar with an effective coupling constant of 5.5 ps.^{114} This first simulation was primarily employed for characterizing the AA pressure–density equation of state.

The second AA simulation modeled a liquid slab in the constant NVT ensemble at a temperature of 300 K. This simulation considered 1348 molecules in a fixed 4 × 4 × 20 nm^{3} box with periodic boundary conditions in all three directions. Following 5 ns of equilibration, we simulated this system for an additional 50 ns while sampling configurations and forces every 1 ps. The sampled configurations and forces were then employed to calculate the MS-CG potential.

#### 2. CG parameterization

We parameterized the bottom-up potentials according to the MS-CG method,^{103–105} as implemented in the Bottom-up Open-source Coarse-graining Software (BOCS) package.^{115} We modified the BOCS package to parameterize both LD and SG potentials. We employed cubic spline basis functions to represent the pair, *F*_{2}(*r*), and local density, *F*_{ρ}(*ρ*), force functions, as well as the SG coefficient function, *U*_{∇}(*ρ*). The pair force functions were represented with a grid spacing of 0.01 nm and were truncated at 1.4 nm. The LD force and SG coefficient functions were represented with a grid spacing of 0.1 nm^{−3} and were calculated for the entire range of local densities sampled in the AA simulations. As a technical note, these calculations included the self-term in defining the local density, i.e., the *J* = *I* contribution $w\u0304(0)$ in Eq. (5). This has no effect upon the results, since it results in corresponding shifts in the arguments of *U*_{ρ} and *U*_{∇}. For the sake of clarity, we have eliminated this self-term when presenting results since it depends upon the Lucy function cutoff, *R*_{c}.

As noted above, we defined the local density with the Lucy function in Eq. (31), which depends upon a single cut-off parameter, *R*_{c}. At present, we lack a systematic methodology for automatically determining *R*_{c}. Consequently, we parameterized a series of CG potentials for varying cutoffs and determined the optimal *R*_{c} based upon the resulting pressure–density equation of state, as detailed in Sec. IV C.

#### 3. CG simulations

We simulated all MS-CG models with a version of LAMMPS^{96} that we modified to support LD and SG potentials. These simulations employed a time step of *δt* = 1.0 fs and maintained a constant temperature *T* = 300 K via the Nosé–Hoover thermostat^{116,117} with a coupling constant of 100.0 fs and a default chain length of 3.^{118} Constant NPT simulations of the MS-CG model employed the MTTK barostat with a coupling constant of 1000.0 fs to maintain an external pressure of 0.986 atm (∼1.0 bar).^{97,98} We first equilibrated the CG simulations for 1.0 ns. We then continued these simulations for an additional 20 ns and sampled configurations every 0.5 ps.

### C. Results

This section presents results for several MS-CG models of acetic acid.^{58} Each CG model represented each molecule with a single site located at its mass center. Before considering these MS-CG models, we first briefly consider the ensemble obtained by mapping a constant NPT simulation of the AA model to the one-site CG representation. The solid black curve in Fig. 6(a) presents the resulting mapped rdf. The mapped rdf features a prominent first peak near *r* ≈ 0.5 nm with a rather broad shoulder toward larger distances and a first minimum near *r* ≈ 0.7 nm. The solid black curve in Fig. 6(b) presents the distribution of global densities, *ρ* = *N*/*V*, sampled by this AA simulation. This distribution is very sharp and peaked near *ρ* ≈ 11 nm^{−3}.

Figure 6 also presents the un-normalized Lucy weighting function, *w*(*r*), that is used to define the local density, *ρ*_{I}, around each site according to Eq. (31). This Lucy function depends upon a single cut-off parameter, *R*_{c}, beyond which *w*(*r*) = 0. The dashed curves in Fig. 6(a) present Lucy functions with *R*_{c} = 0.55, 0.70, and 0.90 nm. These distances roughly correspond to the first peak of the rdf, the first minimum of the rdf, and the second maximum of the rdf, respectively. The dashed curves in Fig. 6(b) present the sampled distribution of local densities defined by each cutoff. When the Lucy function is defined with a cutoff *R*_{c} = 0.55 nm, the local density only reflects molecules in close contact. Given such a short cutoff, the local density almost vanishes for some molecules in the homogeneous bulk liquid. As the Lucy function cutoff increases, the local density includes contributions from more distant molecules and the distribution systematically shifts to larger densities. However, even for *R*_{c} = 0.90 nm, the local density distribution is very different from the global density distribution.

We now consider several MS-CG models that employ approximate potentials, *U*, with varying complexity. The pair-additive model described all interactions with conventional pair potentials, while the LD models employed both pair and LD potentials. Finally, the SG models employed pair, LD, and SG potentials to model interactions. While Fig. 6 presents results from a constant NPT AA simulation, we calculated these MS-CG potentials based upon the configurations and forces sampled from a constant NVT AA simulation of a slab of liquid acetic acid in coexistence with its vapor phase.

Figure 6 demonstrates that the local density distribution in the mapped ensemble varies significantly with the Lucy function cutoff, *R*_{c}. We have previously demonstrated that the MS-CG potential and the accuracy of the MS-CG model can depend quite sensitively upon this cutoff.^{90} Accordingly, we first calculated MS-CG potentials for a range of cutoffs. We then simulated each potential in the constant NPT ensemble at 300 K and 1 bar external pressure. We estimated the equilibrium density and compressibility of each CG model from this bulk simulation. We determined the optimal MS-CG model, which corresponds to the optimal *R*_{c}, as the model that most accurately reproduced the density and compressibility of the AA model.

Figure 7 presents the results employed in optimizing *R*_{c}. The lone X presents the equilibrium density and compressibility of the AA model. The lone O presents the corresponding results for an MS-CG model that employed only a pair potential. This pair-additive MS-CG model significantly overestimates both the density and compressibility of the AA model. The series of diamonds connected by a maroon curve corresponds to a series of MS-CG models that employed both pair and LD potentials while employing the indicated cutoffs for defining the local density. Within this series of LD models, the cutoffs *R*_{c} = 0.55 and 0.90 nm most accurately reproduced the AA density and compressibility, respectively. In the following, we refer to these models as the *ρ* 0.55 and *ρ* 0.90 models, respectively. The series of plus signs connected by a dark green curve corresponds to a series of MS-CG models that employed pair, LD, and also SG potentials. Within this series of SG models, the cutoffs *R*_{c} = 0.70 and 0.72 nm quite accurately reproduced both the AA density and compressibility. In the following, we refer to these models as the ∇ 0.70 and ∇ 0.72 models, respectively.

The top, middle, and bottom panels of Fig. 8 present the calculated pair, LD, and SG potentials, respectively, for the selected MS-CG models. In particular, the black curve in the top panel presents the pair potential for the pair-additive MS-CG model. This pair potential includes two attractive minima near 0.5 and 1.0 nm that are separated by a barrier with a peak near 0.7 nm.

The brown curves in the top and central panel present the pair and LD potentials for the *ρ* 0.55 model. Since the pair and LD potentials are not unique, we have transformed the calculated potentials in order to maximize the overlap of the pair forces with the pair force of the pair-additive model.^{90} The resulting pair potential for the *ρ* 0.55 model is very similar to the pair potential for the pair-additive model. The LD potential for the *ρ* 0.55 model is centered near *ρ*_{I} ≈ 2.5 nm^{−3}, which is slightly higher than the mean local density of the AA model for this cutoff. This LD potential is quite broad and varies by ∼2 *k*_{B}*T* over almost the entire range of local densities sampled in the bulk or at the interface, 0.5 nm^{−3} < *ρ*_{I} < 5 nm^{−3}. In particular, the LD potential calculated for small local densities reflects both the vapor phase and also the liquid–vapor interface sampled by the AA slab simulation. The oscillations that appear in *U*_{ρ} for large *ρ* likely reflect statistical noise for large local densities that are very rarely sampled in the AA simulations and may also possibly reflect the boundary conditions imposed in calculating the potential.

The red curves in the top and central panel present the pair and LD potentials for the *ρ* 0.90 model. The pair potential for the *ρ* 0.90 model is significantly more repulsive and lacks the attractive minimum near 0.5 nm. This more repulsive pair potential is compensated by a more attractive LD potential that is steeper and centered at relatively high local densities, *ρ*_{I} ≈ 10 nm^{−3}. Figure 6(b) indicates that such high local densities, *ρ*_{I} ≥ 10 nm^{−3}, are very rarely sampled by the AA simulation.

The green and blue curves in Fig. 8 present the potentials for the ∇ 0.70 and ∇ 0.72 models, respectively. Since they employ very similar Lucy function cutoffs for defining the local density, the calculated potentials are very similar for these SG models. Interestingly, the pair potentials are most attractive for these models. The corresponding LD potentials are exceptionally broad with a minimum near *ρ*_{I} ≈ 5 nm^{−3}, which is quite close to the center of the AA local density distribution for this cutoff. The SG coefficient function, *U*_{∇}, is negative at nearly all local densities, which suggests that *U*_{SG} may slightly enhance the local structure of the bulk fluid. Moreover, *U*_{∇} slowly decreases as the local density decreases over the entire range of local densities sampled by the AA model. The SG coefficient features a well near *ρ*_{1} ≈ 0.5 nm^{−3} and an even more pronounced well near *ρ*_{1} ≈ 0. These minima presumably act to stabilize and shape the liquid–vapor interface, although it is possible that the well near *ρ*_{1} ≈ 0 is a numerical or statistical artifact for rarely sampled local densities.

Figure 9 presents the rdfs and pressure–density equations of state obtained from constant NPT bulk simulations of the optimized MS-CG models. All of the models accurately reproduce the bulk pair structure of the AA model. As indicated in Fig. 7, the MS-CG models generate rather diverse pressure–density equations of state. In particular, the MS-CG models with SG potentials most accurately reproduce the AA equation of state.

Finally, Fig. 10 and Table II present results for constant NVT simulations of a slab of liquid acetic acid in coexistence with its vapor. The solid black line in the top panel of Fig. 10 presents the density profile from the AA simulation that was employed in calculating the MS-CG potentials. The AA density profile indicates a slight oscillation near the liquid–vapor interface with a slight depression followed by a larger peak as molecules move from the liquid phase into the vapor phase. The width of the AA liquid–vapor interface is *δ*_{AA} ≈ 0.55 nm.

Model . | δ (nm)
. | γ (mN/m)
. |
---|---|---|

AA | 0.553 | 39.2 ± 5.3 |

Pair | 1.765 | 38.4 ± 0.34 |

ρ 0.55 | 1.390 | 33.4 ± 0.31 |

ρ 0.90 | 0.479 | 49.8 ± 0.37 |

∇ 0.70 | 0.603 | 55.5 ± 0.31 |

∇ 0.72 | 0.572 | 56.2 ± 0.37 |

Model . | δ (nm)
. | γ (mN/m)
. |
---|---|---|

AA | 0.553 | 39.2 ± 5.3 |

Pair | 1.765 | 38.4 ± 0.34 |

ρ 0.55 | 1.390 | 33.4 ± 0.31 |

ρ 0.90 | 0.479 | 49.8 ± 0.37 |

∇ 0.70 | 0.603 | 55.5 ± 0.31 |

∇ 0.72 | 0.572 | 56.2 ± 0.37 |

The dashed black curve demonstrates that because it was parameterized from the slab AA simulation, the pair-additive MS-CG model does allow for liquid–vapor coexistence. If the pair-additive model had been parameterized from the bulk simulation, it would dramatically overestimate the internal pressure and, consequently, vaporize under these conditions.^{50,90} Nevertheless, this model significantly overestimates both the bulk density and also the width of the interface. In contrast to the pair-additive model, the LD and SG models all reproduce the bulk density quite accurately. However, they reproduce the interfacial profile with varying accuracy.

In order to better compare the simulated interfacial profiles, the bottom panel of Fig. 10 presents the errors in the simulated CG profiles. The *ρ* 0.55 LD model very accurately reproduces the bulk density but very poorly reproduces the interfacial profile. This model generates a very broad and featureless interface that completely fails to reproduce either the density depression or peak in the AA interfacial profile. In contrast, the *ρ* 0.90 LD model slightly underestimates the bulk density but much more accurately reproduces the interfacial width and also qualitatively reproduces both the depression and the peak in the AA interfacial profile. Thus, it appears more important for LD models to reproduce the compressibility than the bulk density in order to reproduce the AA interfacial structure.

The inclusion of SG potentials allows an even more accurate description of the liquid–vapor interfacial profile. The SG models almost quantitatively reproduce both the density of the bulk phase and also the width of the interface. They also accurately reproduce the interfacial profile, although they slightly underestimate the density depression and peak at the interface. Consequently, the inclusion of SG potentials improves the description of both the bulk pressure–density equation of state and also the interfacial structure.

Table II also presents the simulated liquid–vapor surface tension for each acetic acid model. In each case, we calculated the excess pressure tensor, including the contributions of the LD and SG potentials, according to Eq. (13) and then calculated the surface tension according to Eq. (25). Rather surprisingly, although it poorly describes both the pressure–density equation of state and the interfacial profile, the pair-additive MS-CG model accurately reproduces the AA interfacial tension. By comparison, neither LD model accurately describes the AA surface tension. The *ρ* 0.55 model, which accurately describes the bulk density but overestimates both the compressibility and also the width of the interfacial profile, underestimates the surface tension by 15%. Conversely, the *ρ* 0.90 model, which underestimates the bulk density but more accurately reproduces the compressibility and interfacial profile, overestimates the surface tension by 27%. Moreover, the SG models, which most accurately reproduce the bulk pressure–density equation of state and interfacial profile, overestimate the surface tension by over 40%. Thus, the CG models that most accurately describe the internal pressure and interfacial profile provide the least accurate estimate of the AA interfacial tension.

We do not currently have an explanation for this surprising trend. These results are reminiscent of the tendency of structure-based pair potentials to dramatically overestimate the internal pressure of bulk fluids.^{40,43} Fundamentally, this error in the internal pressure results from failing to properly treat the density-dependence of the many-body PMF.^{14,32,37,38} Thus, accurately modeling the interfacial tension may require more carefully considering the dependence of the many-body PMF upon the interfacial area. Conversely, it may be possible to more accurately model the interfacial tension by parameterizing the CG potentials in a more “informative” ensemble^{119} that provides microscopic information about surface area fluctuations. Finally, we note that it is certainly possible to introduce a new “area potential” in analogy to volume potentials^{47–49} or an operator for modeling the interfacial tension.^{120–124} However, these approaches will not ensure that the CG model accurately samples interfacial fluctuations.

## V. CONCLUSIONS

The present work investigated the use of SG potentials for simulating CG models of interfaces. By explicitly modeling the free energy cost associated with density gradients, these SG potentials may enable more efficient simulations of complex phase behavior with lower resolution models. We first systematically examined the influence of SG potentials upon top-down DPD models. By considering various functional forms for *U*_{∇}, we demonstrated that SG potentials can be used to tune interfacial properties while only slightly perturbing the bulk properties of DPD models. In addition, we employed the bottom-up MS-CG variational principle to determine an optimal SG potential based upon AA simulations of acetic acid. When using an appropriate length scale to define the local density, the SG potential significantly improved the fidelity of the MS-CG model for describing both the bulk pressure–density equation of state and also the structure of the liquid–vapor interface. Thus, SG potentials appear to hold considerable promise for improving the accuracy, efficiency, and transferability of CG models for interfacial phenomena.

The present work indicates several directions for improving bottom-up CG models of interfaces and inhomogeneous systems. In particular, the large discrepancies in the surface tension of the LD and SG models clearly demand attention. Future studies should explore the possibility of parameterizing these potentials from more informative simulations that sample larger interfacial fluctuations.^{119} Moreover, future studies should more rigorously consider the influence of interfaces upon the many-body PMF.^{14,91} Furthermore, the present work further highlights the sensitivity of MS-CG models to the length scale, *R*_{c}, employed in defining the local density. If the local density is defined over a sufficiently large length scale, then LD potentials can be introduced to modify thermodynamic properties without impacting the bulk structure.^{89} Unfortunately, the resulting models then tend to demonstrate artifacts at interfaces where the density rapidly varies.^{90} Conversely, if the local density is defined on shorter length scales, then the local density can become correlated with the pair structure. In this case, LD and SG potentials can provide a very accurate description of the bulk and interfacial structure as well the pressure–density equation of state.^{91} Unfortunately, in this case, the accuracy of the MS-CG model can depend quite sensitively upon *R*_{c}.^{90} Clearly, it would be very helpful to develop theory and computational methods to understand and optimize this cutoff.^{83} Additionally, future work should further consider CG potentials that depend upon other order parameters and their gradients.

More generally, the present work also suggests several additional directions for modeling interfaces and phase behavior. For instance, future studies should develop an appropriate stability criterion^{101,102} for simulating many-body potentials that depend upon both the local density and its gradient. Moreover, future studies should more deeply mine the relationship between classical DFT and bottom-up CG models. Similarly, future studies should explore the relationship between Eq. (2) and the effective Hamiltonian employed in mean-field Landau models for inhomogeneous systems.^{92} In particular, it may be possible to design LD and SG potentials to simulate the formation of novel phases with a large number of interfaces. Thus, it is hoped that these SG potentials may prove fruitful for gaining new insight into complex interfaces and phase behavior.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge financial support from the National Science Foundation (Grant Nos. CHE-1565631 and CHE-1856337). Parts of this research were conducted with Advanced CyberInfrastructure computational resources provided by the Institute for Computational and Data Sciences at the Pennsylvania State University (http://icds.psu.edu). In addition, parts of this research were conducted with XSEDE resources awarded by Grant No. TG-CHE170062. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation (Grant No. ACI-1548562).^{125} W.G.N. also gratefully acknowledges helpful discussions with Markus Deserno.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: CONSERVATION OF MOMENTUM

In this appendix, we explicitly demonstrate that Eqs. (10) and (11) conserve both linear and angular momentum. It is convenient to first note that **F**_{IJ} may be expressed as

where

correspond to the three rows contributing to **F**_{∇;IJ} in Eq. (11). Note that the scalars *F*_{2}(*R*_{IJ}), *F*_{ρ;IJ}(**R**), *S*_{1;IJ}(**R**), and *S*_{2;IJ}(**R**) are all symmetric with respect to *I* ↔ *J*. Conversely, **A**_{∇;IJ} is a vector that is antisymmetric with respect to *I* ↔ *J*. Consequently, Eq. (A1) implies that **F**_{IJ}(**R**) = −**F**_{JI}(**R**). This ensures that the net force on the system of particles vanishes and that linear momentum is conserved.

We next consider conservation of angular momentum. Because the total force is pairwise additive, the total torque may be decomposed into pair-additive contributions as

where

and **R**_{IJ} = **R**_{I} − **R**_{J}. Because **N**_{IJ} = **N**_{JI},

To proceed further, we employ Eq. (7) to evaluate **A**_{∇;IJ},

The pair torque may then be expressed as

where

Note that Eq. (A9) preserves the *I* ↔ *J* symmetry, but **A**_{I|JK} is antisymmetric with respect to *J* ↔ *K*, i.e., **A**_{I|JK} = −**A**_{I|KJ}. Substituting Eq. (A9) into Eq. (A7), one obtains

where the sum is evaluated over all triples of three distinct particles. For any fixed *I*, the sum over *J* and *K* vanishes due to the antisymmetry of **A**_{I|JK}, i.e., ∑_{J≠K(≠I)}**A**_{I|JK} = **0**. Consequently, the total torque vanishes, **N** = **0**, and the angular momentum of the system is preserved.

### APPENDIX B: VALIDATING THE MS-CG IMPLEMENTATION

In order to validate our numerical implementation of the SG potential, we consider a toy model with a known potential of the form given by Eq. (2). Given the simulated trajectory for this toy model, we attempted to recover the original potential by determining the coefficients that minimized *χ*^{2}. Since the simulated force field lies within the basis set included in the MS-CG calculation, *χ*^{2} should attain a minimum value of 0 when the calculated CG force field perfectly matches the simulated forces, i.e., when the variational calculation recovers the simulated force field.

Figure 11 presents the results of this calculation. Specifically, the left panel presents the pair potential, *U*_{2}(*r*), the center panel presents the LD potential, *U*_{ρ}(*ρ*), and the right panel presents the coefficient of the SG potential, *U*_{∇}(*ρ*). The dashed green curves present the simulated functions, while the solid black curves present the functions obtained by minimizing *χ*^{2}. The MS-CG calculation quantitatively recovers *U*_{∇} for the range of local densities that are significantly sampled by the simulation. However, the calculated pair and LD potentials differ noticeably from the simulated potentials. These discrepancies arise because the forces resulting from *U*_{2} and *U*_{ρ} are not independent, as we have previously discussed in Ref. 90. Accordingly, we applied the transformation described in Ref. 90 to the calculated pair and LD potentials. The solid red curves in Fig. 11 present the resulting potentials. Aside from the discrepancies at local densities that are not sampled in the simulation, the transformed potentials quantitatively match the original simulated potentials. Thus, the MS-CG variational calculation quantitatively recovers the simulated pair, LD, and SG forces. This consistency validates our implementation of the SG potential into both the LAMMPS simulation engine and the BOCS force-matching code.