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.

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 as10–13 

exp[βW(R)]=V(nN)Vndrexp[βu(r)]δ(M(r)R).
(1)

Here, β = 1/kBT, u is the AA potential, r is the configuration for n atoms, R is the configuration for N CG sites, and M: rR = 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 transferability22–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) models66 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.

We consider potentials with the form

U(R)=Upair(R)+ULD(R)+USG(R).
(2)

The first term is a sum of central pair potentials, U2, that depend upon the distance, RIJ, between each pair of sites,

Upair(R)=(I,J)U2(RIJ).
(3)

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,

ULD(R)=IUρ(ρI).
(4)

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

ρI=J(I)w̄(RIJ),
(5)

where w̄(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] = ∫dr4πr2w(r) is the normalizing volume, and w̄(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̄(0)=[w]1, 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,

USG(R)=IU(ρI)IρI2.
(6)

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,

IρI=J(I)w̄(RIJ)R̂IJ,
(7)

where ∇I = /RI and R̂IJ=RIRJ/RIJ is the unit vector pointing from site J toward site I. Because RIJ, ρI, and Iρ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,

FI(R)=IU(R)=J(I)FIJ(R).
(8)

The total pair force may be expressed as

FIJ(R)=F2(RIJ)+Fρ;IJ(R)R̂IJ+F;IJ(R).
(9)

The first contribution in Eq. (9) comes from Upair: F2(R) = −dU2(R)/dR. The second contribution comes from ULD and acts along R̂IJ with a magnitude,

Fρ;IJ(R)=Fρ(ρI)+Fρ(ρJ)w̄(RIJ),
(10)

where Fρ(ρ) = −dUρ(ρ)/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 USG and is more complex,

F;IJ(R)=w̄(RIJ)F(ρI)AI2+F(ρJ)AJ2R̂IJ2g(RIJ)U(ρI)AIU(ρJ)AJR̂IJR̂IJ2f(RIJ)U(ρI)AIU(ρJ)AJ,
(11)

where we have defined F(ρ) = −dU(ρ)/dρ, AI = ∇IρI, f(R)=w̄(R)/R, and g(R)=w̄(R)f(R). Because each line of Eq. (11) is antisymmetric with respect to exchanging I and J, USG generates equal and opposite forces between each pair of particles: F∇;JI = −F∇;IJ. However, the pair force, F∇;IJ, from USG does not act along the interparticle vector, R̂IJ. While the first two rows of Eq. (11) act along R̂IJ, the third line does not.

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

FIJ(R)=FJI(R),
(12)

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

Pxs; ij(R)1VIêiRIFIêj=1VIJ(>I)êiRIJFIJêj,
(13)

where êi indicates the unit vector in the ith 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.

DPD models66 generally assume that particles follow Newton’s equations of motion

ddtRI=VI,ddtVI=FI/mI
(14)

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

FI=J(I)FIJC+FIJD+FIJR.
(15)

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,

FIJD=γw2(RIJ)R̂IJVIJR̂IJ,
(16)
FIJR=σw2(RIJ)θIJR̂IJ.
(17)

Here, VIJ = VIVJ is the velocity of site I with respect to site J, θIJ is a pairwise white noise satisfying Gaussian statistics: θIJ(t)=0 and θIJ(t)θKL(t)=δIKδJL+δILδJKδ(tt), and

w2(R)=1R/Rc;22Θ(Rc;2R)
(18)

is a quadratic pairwise weighting function expressed in terms of the Heaviside function, Θ, such that the pair force vanishes between particles with RIJ > Rc;2. According to the fluctuation–dissipation relation, the viscous and random forces are related by σ2 = 2γkBT.95 

MDPD models66 adopt a conservative potential, U = Upair + ULD, which corresponds to Eq. (2) with USG = 0. The MDPD potential often adopts a particularly simple form.71–75 Specifically, the pair potential is assumed proportional to w2,

U2(R)=12ARc;2w2(R).
(19)

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

Uρ(ρI)=π30BRc;ρ4ρI2,
(20)

where the local density

ρI=J(I)w̄ρ(RIJ)
(21)

is defined in terms of a second quadratic weighting function

wρ(R)=1R/Rc;ρ2Θ(Rc;ρR)
(22)

that acts on a shorter length scale, Rc;ρ = 0.75Rc;2, than w2. As before, w̄ρ(r)=w̄ρ(r)/[wρ] and [wρ] = ∫dr4πr2wρ(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

FIJC(R)=A1RIJ/Rc;2Θ(Rc;2RIJ)+BρI+ρJ1RIJ/Rc;ρΘ(Rc;ρRIJ)R̂IJ,
(23)

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, USG, upon MDPD models. In particular, we consider three simple forms for U,

U(ρI)=c,constant,mρcρIΘρcρI,linear,kρcρI2ΘρcρI,quadratic.
(24)

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.

DPD studies often adopt dimensionless quantities. We define the length scale by the pair potential cutoff, Rc;2 = 1, the energy scale by thermal energy, kBT = 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 Rc;ρ = 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 × 106 time steps with δt = 0.001. We then simulated the system for 20 × 106 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, USG, upon the behavior of MDPD models. Accordingly, we define a reference baseline model with USG = 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.

TABLE I.

A summary of the parameters employed in the MDPD simulations.

SetFixed param.Var. param.Sample set
B = 25 A −30, −35, −37.5, −42.5, −45, −50 
A = −40 B 15, 20, 22.5, 27.5, 30, 35 
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 
SetFixed param.Var. param.Sample set
B = 25 A −30, −35, −37.5, −42.5, −45, −50 
A = −40 B 15, 20, 22.5, 27.5, 30, 35 
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 = Upair + ULD and USG = 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 USG. 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(ρI)=kρcρI2Θρcρ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}.

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 w2(R) is a quadratic weighting function that vanishes at a cutoff Rc;2 = 1. The baseline model also employs a standard MDPD LD potential, Uρ(ρI)=π30BRc;ρ4ρI2, where the local density, ρI, around particle I is defined by a quadratic weighting function, wρ(RIJ), that vanishes at a cutoff Rc;ρ = 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, Rc;2, all densities are reported in units of Rc;23 and ∇IρI is reported in units of Rc;24. Since we define energy with respect to kBT, 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, w2, and LD weighting functions, wρ, which vary quadratically with distance, r, out to their respective cutoffs, Rc;2 = 1 and Rc;ρ = 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, Pext = 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.

FIG. 1.

Analysis of constant pressure simulations for the baseline reference MDPD model. Panel (a) presents the un-normalized weighting functions w2 and wρ as solid and dashed curves, respectively. Panels (b) and (c) present the simulated rdf and pressure–density equation of state, respectively. Panel (d) presents the simulated distribution of global, ρ = n/V, and local, ρ1, densities as the solid and dashed curves, respectively. Note that the global density distribution has been scaled to match the height of the local density distribution.

FIG. 1.

Analysis of constant pressure simulations for the baseline reference MDPD model. Panel (a) presents the un-normalized weighting functions w2 and wρ as solid and dashed curves, respectively. Panels (b) and (c) present the simulated rdf and pressure–density equation of state, respectively. Panel (d) presents the simulated distribution of global, ρ = n/V, and local, ρ1, densities as the solid and dashed curves, respectively. Note that the global density distribution has been scaled to match the height of the local density distribution.

Close modal

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 barostat97,98 to sample volume fluctuations, Trofimov et al. previously observed similar oscillations in constant pressure MDPD simulations99 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, ρ̄(z), and local, ρ̄1(z), densities as a function of the slab z-coordinate. The slab simulation generates coexistence between a liquid phase with a density ρ̄liq6 and a vapor phase with a density ρ̄vap5×105. The average local density of the liquid phase is ρ̄1;liq2.5. A rather sharp liquid–vapor interface forms near z ≈ ±7. Interestingly, the global density, ρ̄(z), slightly oscillates before the interface and then rapidly plummets at the interface. In contrast, the local density, ρ̄1(z), does not oscillate but transitions rather smoothly over a considerably wider region.

FIG. 2.

Analysis of global, ρ, and local, ρ1 densities sampled in a constant NVT slab simulation of the baseline reference MDPD model. Panel (a) presents the global (solid curve) and local (dashed curve) density profiles as a function of the z coordinate through the slab. Panel (b) presents an intensity plot quantifying the joint distribution, P(z, ρ1), of z coordinates and local densities, ρ1, sampled by particles. Panel (c) presents an analogous intensity plot quantifying the joint distribution, P(z, |∇ρ1|), for the particle z coordinate and the magnitude of the gradient in the local density, |∇ρ1|.

FIG. 2.

Analysis of global, ρ, and local, ρ1 densities sampled in a constant NVT slab simulation of the baseline reference MDPD model. Panel (a) presents the global (solid curve) and local (dashed curve) density profiles as a function of the z coordinate through the slab. Panel (b) presents an intensity plot quantifying the joint distribution, P(z, ρ1), of z coordinates and local densities, ρ1, sampled by particles. Panel (c) presents an analogous intensity plot quantifying the joint distribution, P(z, |∇ρ1|), for the particle z coordinate and the magnitude of the gradient in the local density, |∇ρ1|.

Close modal

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, ρ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

γ=12LzPzz12Pxx+Pyy,
(25)

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.

FIG. 3.

Density profiles from constant NVT slab simulations for families of MDPD models. The top row presents results for the A, B, and c model families from left to right. The bottom row presents results for the L1.5, L2.5, Q1.2, and Q1.5 model families from left to right. In each case, the adjacent legend indicates the values of the varying parameter.

FIG. 3.

Density profiles from constant NVT slab simulations for families of MDPD models. The top row presents results for the A, B, and c model families from left to right. The bottom row presents results for the L1.5, L2.5, Q1.2, and Q1.5 model families from left to right. In each case, the adjacent legend indicates the values of the varying parameter.

Close modal
FIG. 4.

Analysis of constant NVT slab simulations for families of MDPD models. The top, middle, and bottom rows present the bulk liquid density, the bulk vapor density, and the interfacial tension, respectively. The columns present results for the A, B, c, L1.5, L2.5, Q1.2, and Q1.5 model families from left to right. The x-axis for each column indicates the value of the varying parameter. The error bars indicate estimates of the standard error. In each panel, the horizontal dotted green line indicates the corresponding property for the baseline reference model.

FIG. 4.

Analysis of constant NVT slab simulations for families of MDPD models. The top, middle, and bottom rows present the bulk liquid density, the bulk vapor density, and the interfacial tension, respectively. The columns present results for the A, B, c, L1.5, L2.5, Q1.2, and Q1.5 model families from left to right. The x-axis for each column indicates the value of the varying parameter. The error bars indicate estimates of the standard error. In each panel, the horizontal dotted green line indicates the corresponding property for the baseline reference model.

Close modal

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(ρI)=kρcρI2Θρcρ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.

FIG. 5.

Panels (a) and (b) present simulated pressure–density equations of state for the Q1.2 and Q1.5 model families, respectively. The black curve in each panel presents the equation of state for the baseline reference MDPD model. The colors indicate the value of the parameter k.

FIG. 5.

Panels (a) and (b) present simulated pressure–density equations of state for the Q1.2 and Q1.5 model families, respectively. The black curve in each panel presents the equation of state for the baseline reference MDPD model. The colors indicate the value of the parameter k.

Close modal

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

U(R)=ζλUζ(ψζλ(R)),
(26)

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 approach103–105 determines the potential that optimally approximates the many-body PMF by minimizing the force-matching functional,106,107

χ2[U]=13NIfI(r)FI(M(r))2,
(27)

where fI(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,,N, in an abstract vector space with FI(R) = −∇IU(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) = −dUζ(x)/dx, that is represented with simple basis functions, fζd(x), and constant coefficients, ϕζd,

Fζ(x)=dϕζdfζd(x).
(28)

The interaction potential Uζ determines a force on each site, I,

Fζ;I(R)=dϕζdGζd;I(R),
(29)

where

Gζd;I(R)=λfζd(ψζλ(R))Iψζλ(R)
(30)

and Gζd={Gζd;I(R)}I=1,,N defines a force field basis vector. The set of basis vectors, {Gζ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, ϕζd*, that minimize χ2.109 

Returning to Eq. (2), both Upair and ULD 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

w(R)=1R/Rc31+3R/RcΘRcR
(31)

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

USG(R)=IU(ρI)IρI2,

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 USG, we first define

F;I(R)IUSG(R)=J(I)F;IJ(R).
(32)

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

F;IJ(R)=F(ρI)GIJ;1(R)+F(ρJ)GIJ;2(R)+U(ρI)HIJ;1(R)+U(ρJ)HIJ;2(R),
(33)

where

GIJ;1(R)=w̄(RIJ)AI2R̂IJ,GIJ;2(R)=w̄(RIJ)AJ2R̂IJ,HIJ;1(R)=2g(RIJ)AIR̂IJR̂IJ2f(RIJ)AI,HIJ;2(R)=+2g(RIJ)AJR̂IJR̂IJ+2f(RIJ)AJ.
(34)

The function U(ρ) is represented by simple basis functions, fd(ρ), with constant coefficients, φd,

U(ρ)=dφdfd(ρ),
(35)

which determines a corresponding representation of F(ρ),

F(ρ)=dφdfd(ρ).
(36)

Given this representation for U and F, F∇;IJ may be expressed as

F;IJ(R)=dφdGd;IJ(R),
(37)

where

Gd;IJ(R)=fd(ρI)HIJ;1(R)fd(ρI)GIJ;1(R)+fd(ρJ)HIJ;2(R)fd(ρJ)GIJ;2(R).
(38)

Finally, this allows for a basis set expansion of the forces on each site I due to USG,

F;I(R)=dφdGd;I(R),
(39)

where

Gd;I(R)=J(I)Gd;IJ(R)
(40)

defines the force field basis vector, Gd={Gd;I(R)}I=1,,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.

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.4110 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 nm3 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, F2(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̄(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, Rc.

As noted above, we defined the local density with the Lucy function in Eq. (31), which depends upon a single cut-off parameter, Rc. At present, we lack a systematic methodology for automatically determining Rc. Consequently, we parameterized a series of CG potentials for varying cutoffs and determined the optimal Rc 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 LAMMPS96 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 thermostat116,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.

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.

FIG. 6.

Analysis of the mapped ensemble obtained by mapping AA constant NPT simulations of acetic acid to the 1-site mass center CG representation. Panel (a): The solid black curve presents the mapped rdf. The dashed curves present un-normalized Lucy weighting functions, w(r), with the cutoffs Rc = 0.55 nm (red), 0.70 nm (green), and 0.90 nm (blue). Panel (b): The solid black curve presents the distribution of global densities, ρ = N/V, sampled by the AA constant NPT simulation. The dashed curves present the sampled distribution of local densities, ρ1, defined by the various Lucy functions in panel (a). The magnitude of the global density distribution has been scaled for comparison with the local density distributions.

FIG. 6.

Analysis of the mapped ensemble obtained by mapping AA constant NPT simulations of acetic acid to the 1-site mass center CG representation. Panel (a): The solid black curve presents the mapped rdf. The dashed curves present un-normalized Lucy weighting functions, w(r), with the cutoffs Rc = 0.55 nm (red), 0.70 nm (green), and 0.90 nm (blue). Panel (b): The solid black curve presents the distribution of global densities, ρ = N/V, sampled by the AA constant NPT simulation. The dashed curves present the sampled distribution of local densities, ρ1, defined by the various Lucy functions in panel (a). The magnitude of the global density distribution has been scaled for comparison with the local density distributions.

Close modal

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, Rc, beyond which w(r) = 0. The dashed curves in Fig. 6(a) present Lucy functions with Rc = 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 Rc = 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 Rc = 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, Rc. 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 Rc, as the model that most accurately reproduced the density and compressibility of the AA model.

Figure 7 presents the results employed in optimizing Rc. 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 Rc = 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 Rc = 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.

FIG. 7.

Sensitivity of the MS-CG models to the cutoff, Rc, employed in the Lucy function that defines the local density. The x- and y-axes present the density and compressibility estimated from constant NPT simulations of each model. The single “x” presents results for the AA model of acetic acid. The single “o” presents results for the pair-additive MS-CG model. The diamonds connected by the maroon lines present results for various LD models, while the crosses connected by the dark green lines present results for the various SG models. The LD and SG models are labeled by the cutoff, Rc, employed in defining the local density.

FIG. 7.

Sensitivity of the MS-CG models to the cutoff, Rc, employed in the Lucy function that defines the local density. The x- and y-axes present the density and compressibility estimated from constant NPT simulations of each model. The single “x” presents results for the AA model of acetic acid. The single “o” presents results for the pair-additive MS-CG model. The diamonds connected by the maroon lines present results for various LD models, while the crosses connected by the dark green lines present results for the various SG models. The LD and SG models are labeled by the cutoff, Rc, employed in defining the local density.

Close modal

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.

FIG. 8.

Calculated potentials for select MS-CG models. Panels (a)–(c) present the calculated functions U2(R), Uρ(ρ), and U(ρ), respectively. The black curve presents the pair potential for the pair-additive model. The brown and red curves present the potentials for the LD models with Rc = 0.55 and 0.90 nm, respectively. The green and blue curves present the potentials for the SG models with Rc = 0.70 and 0.72 nm, respectively. The pair and LD potentials for the LD and SG models were transformed to maximize overlap with the pair force function of the pair-additive model.90 All potentials were calculated from the constant NVT AA slab simulation.

FIG. 8.

Calculated potentials for select MS-CG models. Panels (a)–(c) present the calculated functions U2(R), Uρ(ρ), and U(ρ), respectively. The black curve presents the pair potential for the pair-additive model. The brown and red curves present the potentials for the LD models with Rc = 0.55 and 0.90 nm, respectively. The green and blue curves present the potentials for the SG models with Rc = 0.70 and 0.72 nm, respectively. The pair and LD potentials for the LD and SG models were transformed to maximize overlap with the pair force function of the pair-additive model.90 All potentials were calculated from the constant NVT AA slab simulation.

Close modal

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 kBT 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 USG 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.

FIG. 9.

Analysis of constant NPT bulk simulations for various models. Panels (a) and (b) present the simulated rdf and pressure–density equation of state, respectively. The solid and dashed black curves correspond to the AA model and pair-additive MS-CG model, respectively. The brown and red curves correspond to the LD models with Rc = 0.55 and 0.90 nm, respectively. The green and blue curves correspond to the SG models with Rc = 0.70 and 0.72 nm, respectively. Error bars in the equation of state indicate the simulated standard error.

FIG. 9.

Analysis of constant NPT bulk simulations for various models. Panels (a) and (b) present the simulated rdf and pressure–density equation of state, respectively. The solid and dashed black curves correspond to the AA model and pair-additive MS-CG model, respectively. The brown and red curves correspond to the LD models with Rc = 0.55 and 0.90 nm, respectively. The green and blue curves correspond to the SG models with Rc = 0.70 and 0.72 nm, respectively. Error bars in the equation of state indicate the simulated standard error.

Close modal

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.

FIG. 10.

Analysis of constant NVT slab simulations for various models. Panel (a) presents the simulated density profiles, while panel (b) presents the error in the simulated density profile for each model, i.e., the difference with respect to the AA density profile. The solid and dashed black curves correspond to the AA model and pair-additive MS-CG model, respectively. The brown and red curves correspond to the LD models with Rc = 0.55 and 0.90 nm, respectively. The green and blue curves correspond to the SG models with Rc = 0.70 and 0.72 nm, respectively.

FIG. 10.

Analysis of constant NVT slab simulations for various models. Panel (a) presents the simulated density profiles, while panel (b) presents the error in the simulated density profile for each model, i.e., the difference with respect to the AA density profile. The solid and dashed black curves correspond to the AA model and pair-additive MS-CG model, respectively. The brown and red curves correspond to the LD models with Rc = 0.55 and 0.90 nm, respectively. The green and blue curves correspond to the SG models with Rc = 0.70 and 0.72 nm, respectively.

Close modal
TABLE II.

Width, δ, and surface tension, γ, of the simulated liquid–vapor interface for various acetic acid models. The 5–95 interfacial width, δ, is determined as the distance interval over which the simulated density profile transitions from ρvap + 0.05Δρ to ρvap + 0.95Δρ, where Δρ = ρliqρvap. The uncertainty in the simulated surface tension, γ, is estimated as the standard error obtained by dividing the simulation into five independent blocks.

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” ensemble119 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 potentials47–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.

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, Rc, 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 Rc.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 criterion101,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.

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.

The authors have no conflicts to disclose.

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

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

FIJ(R)=F2(RIJ)+Fρ;IJ(R)+S1;IJ(R)+S2;IJ(R)R̂IJ+A;IJ(R),
(A1)

where

S1;IJ(R)=w̄(RIJ)F(ρI)AI2+F(ρJ)AJ2,
(A2)
S2;IJ(R)=2g(RIJ)U(ρI)AIR̂IJ+U(ρJ)AJR̂JI,
(A3)
A;IJ(R)=2f(RIJ)U(ρI)AIU(ρJ)AJ
(A4)

correspond to the three rows contributing to F∇;IJ in Eq. (11). Note that the scalars F2(RIJ), Fρ;IJ(R), S1;IJ(R), and S2;IJ(R) are all symmetric with respect to IJ. Conversely, A∇;IJ is a vector that is antisymmetric with respect to IJ. Consequently, Eq. (A1) implies that FIJ(R) = −FJI(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

N=IRI×FI=IJ(>I)NIJ,
(A5)

where

NIJ=RIJ×FIJ=RIJ×A;IJ
(A6)

and RIJ = RIRJ. Because NIJ = NJI,

N=12IJ(I)NIJ.
(A7)

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

A;IJ=2f(RIJ)w̄(RIJ)U(ρI)+U(ρJ)R̂IJ2f(RIJ)K(I,J)U(ρI)w̄(RIK)R̂IKU(ρJ)w̄(RJK)R̂JK.
(A8)

The pair torque may then be expressed as

NIJ=2K(I,J)AI|JK+AJ|IK,
(A9)

where

AI|JK=U(ρI)w̄(RIJ)w̄(RIK)R̂IJ×R̂IK.
(A10)

Note that Eq. (A9) preserves the IJ symmetry, but AI|JK is antisymmetric with respect to JK, i.e., AI|JK = −AI|KJ. Substituting Eq. (A9) into Eq. (A7), one obtains

N=(I,J,K)AI|JK+AJ|IK=2(I,J,K)AI|JK,
(A11)

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 AI|JK, i.e., ∑JK(≠I)AI|JK = 0. Consequently, the total torque vanishes, N = 0, and the angular momentum of the system is preserved.

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, U2(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 U2 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.

FIG. 11.

Illustration of validation for the implemented methods in LAMMPS and BOCS. Panels (a)–(c) present results for u2(r), uρ(ρ), and u(ρ), respectively. The green curves present the known potentials simulated in LAMMPS. The black curves present the potential functions calculated by BOCS from the LAMMPS trajectory. The dashed green curves present the potential functions obtained after transforming the calculated potentials to align with the simulated pair force function.90 

FIG. 11.

Illustration of validation for the implemented methods in LAMMPS and BOCS. Panels (a)–(c) present results for u2(r), uρ(ρ), and u(ρ), respectively. The green curves present the known potentials simulated in LAMMPS. The black curves present the potential functions calculated by BOCS from the LAMMPS trajectory. The dashed green curves present the potential functions obtained after transforming the calculated potentials to align with the simulated pair force function.90 

Close modal
1.
M. L.
Klein
and
W.
Shinoda
, “
Large-scale molecular dynamics simulations of self-assembling systems
,”
Science
321
,
798
800
(
2008
).
2.
M.
Deserno
, “
Mesoscopic membrane physics: Concepts, simulations, and selected applications
,”
Macromol. Rapid Commun.
30
,
752
771
(
2009
).
3.
C.
Peter
and
K.
Kremer
, “
Multiscale simulation of soft matter systems
,”
Faraday Discuss.
144
,
9
24
(
2010
).
4.
M. G.
Saunders
and
G. A.
Voth
, “
Coarse-graining methods for computational biology
,”
Annu. Rev. Biophys.
42
,
73
93
(
2013
).
5.
W. G.
Noid
, “
Perspective: Coarse-grained models for biomolecular systems
,”
J. Chem. Phys.
139
,
090901
(
2013
).
6.
S.
Kmiecik
,
D.
Gront
,
M.
Kolinski
,
L.
Wieteska
,
A. E.
Dawid
, and
A.
Kolinski
, “
Coarse-grained protein models and their applications
,”
Chem. Rev.
116
,
7898
7936
(
2016
).
7.
S. J.
Marrink
,
V.
Corradi
,
P. C. T.
Souza
,
H. I.
Ingólfsson
,
D. P.
Tieleman
, and
M. S. P.
Sansom
, “
Computational modeling of realistic cell membranes
,”
Chem. Rev.
119
,
6184
6226
(
2019
).
8.
M.
Giulini
,
M.
Rigoli
,
G.
Mattiotti
,
R.
Menichetti
,
T.
Tarenzi
,
R.
Fiorentini
, and
R.
Potestio
, “
From system modeling to system analysis: The impact of resolution level and resolution distribution in the computer-aided investigation of biomolecules
,”
Front. Mol. Biosci.
8
,
676976
(
2021
).
9.
T.
Sun
,
V.
Minhas
,
N.
Korolev
,
A.
Mirzoev
,
A. P.
Lyubartsev
, and
L.
Nordenskiöld
, “
Bottom-up coarse-grained modeling of DNA
,”
Front. Mol. Biosci.
8
,
645527
(
2021
).
10.
J. G.
Kirkwood
, “
Statistical mechanics of fluid mixtures
,”
J. Chem. Phys.
3
,
300
313
(
1935
).
11.
A.
Liwo
,
C.
Czaplewski
,
J.
Pillardy
, and
H. A.
Scheraga
, “
Cumulant-based expressions for the multibody terms for the correlation between local and electrostatic interactions in the united-residue force field
,”
J. Chem. Phys.
115
,
2323
2347
(
2001
).
12.
R. L. C.
Akkermans
and
W. J.
Briels
, “
A structure-based coarse-grained model for polymer melts
,”
J. Chem. Phys.
114
,
1020
1031
(
2001
).
13.
C. N.
Likos
, “
Effective interactions in soft condensed matter physics
,”
Phys. Rep.
348
,
267
439
(
2001
).
14.
N. J. H.
Dunn
,
T. T.
Foley
, and
W. G.
Noid
, “
Van der Waals perspective on coarse-graining: Progress toward solving representability and transferability problems
,”
Acc. Chem. Res.
49
,
2832
2840
(
2016
).
15.
W. G.
Noid
, “
Systematic methods for structurally consistent coarse-grained models
,”
Methods Mol. Biol.
924
,
487
531
(
2013
).
16.
E.
Brini
,
E. A.
Algaer
,
P.
Ganguly
,
C.
Li
,
F.
Rodríguez-Ropero
, and
N. F. A.
van der Vegt
, “
Systematic coarse-graining methods for soft matter simulations—A review
,”
Soft Matter
9
,
2108
2119
(
2013
).
17.
J. F.
Rudzinski
and
W. G.
Noid
, “
The role of many-body correlations in determining potentials for coarse-grained models of equilibrium structure
,”
J. Phys. Chem. B
116
,
8621
8635
(
2012
).
18.
J. F.
Rudzinski
and
W. G.
Noid
, “
Bottom-up coarse-graining of peptide ensembles and helix-coil transitions
,”
J. Chem. Theory Comput.
11
,
1278
1291
(
2015
).
19.
S. J.
Wörner
,
T.
Bereau
,
K.
Kremer
, and
J. F.
Rudzinski
, “
Direct route to reproducing pair distribution functions with coarse-grained models via transformed atomistic cross correlations
,”
J. Chem. Phys.
151
,
244110
(
2019
).
20.
S. T.
John
and
G.
Csányi
, “
Many-body coarse-grained interactions using Gaussian approximation potentials
,”
J. Phys. Chem. B
121
,
10934
10949
(
2017
).
21.
J.
Wang
,
N.
Charron
,
B.
Husic
,
S.
Olsson
,
F.
Noé
, and
C.
Clementi
, “
Multi-body effects in a coarse-grained protein force field
,”
J. Chem. Phys.
154
,
164113
(
2021
).
22.
J.
Ghosh
and
R.
Faller
, “
State point dependence of systematically coarse-grained potentials
,”
Mol. Simul.
33
,
759
767
(
2007
).
23.
P.
Carbone
,
H. A. K.
Varzaneh
,
X.
Chen
, and
F.
Müller-Plathe
, “
Transferability of coarse-grained force fields: The polymer case
,”
J. Chem. Phys.
128
,
064904
(
2008
).
24.
A.
Chaimovich
and
M. S.
Shell
, “
Anomalous waterlike behavior in spherically-symmetric water models optimized with the relative entropy
,”
Phys. Chem. Chem. Phys.
11
,
1901
1915
(
2009
).
25.
K.
Farah
,
A. C.
Fogarty
,
M. C.
Böhm
, and
F.
Müller-Plathe
, “
Temperature dependence of coarse-grained potentials for liquid hexane
,”
Phys. Chem. Chem. Phys.
13
,
2894
2902
(
2011
).
26.
L.
Lu
and
G. A.
Voth
, “
The multiscale coarse-graining method. VII. Free energy decomposition of coarse-grained effective potentials
,”
J. Chem. Phys.
134
,
224107
(
2011
).
27.
S.
Riniker
,
J. R.
Allison
, and
W. F.
van Gunsteren
, “
On developing coarse-grained models for biomolecular simulation: A review
,”
Phys. Chem. Chem. Phys.
14
,
12423
12430
(
2012
).
28.
B.
Mukherjee
,
L.
Delle Site
,
K.
Kremer
, and
C.
Peter
, “
Derivation of coarse grained models for multiscale simulation of liquid crystalline phase transitions
,”
J. Phys. Chem. B
116
,
8474
8484
(
2012
).
29.
J.
Zhang
and
H.
Guo
, “
Transferability of coarse-grained force field for nCB liquid crystal systems
,”
J. Phys. Chem. B
118
,
4647
4660
(
2014
).
30.
F.
Cao
and
H.
Sun
, “
Transferability and nonbond functional form of coarse grained force field—Tested on linear alkanes
,”
J. Chem. Theory Comput.
11
,
4760
4769
(
2015
).
31.
T. D.
Potter
,
J.
Tasche
, and
M. R.
Wilson
, “
Assessing the transferability of common top-down and bottom-up coarse-grained molecular models for molecular mixtures
,”
Phys. Chem. Chem. Phys.
21
,
1912
1917
(
2019
).
32.
K. M.
Lebold
and
W. G.
Noid
, “
Systematic study of temperature and density variations in effective potentials for coarse-grained models of molecular liquids
,”
J. Chem. Phys.
150
,
014104
(
2019
).
33.
R. J.
Szukalo
and
W.
Noid
, “
Investigation of coarse-grained models across a glass transition
,”
Soft Mater.
18
,
185
(
2020
).
34.
J.
Jin
and
G. A.
Voth
, “
Ultra-coarse-grained models allow for an accurate and transferable treatment of interfacial systems
,”
J. Chem. Theory Comput.
14
,
2180
2197
(
2018
).
35.
J.
Jin
,
A.
Yu
, and
G. A.
Voth
, “
Temperature and phase transferable bottom-up coarse-grained models
,”
J. Chem. Theory Comput.
16
,
6823
6842
(
2020
).
36.
J.
Jin
,
Y.
Han
,
A. J.
Pak
, and
G. A.
Voth
, “
A new one-site coarse-grained model for water: Bottom-up many-body projected water (BUMPer). I. General theory and model
,”
J. Chem. Phys.
154
,
044104
(
2021
).
37.
F. H.
Stillinger
,
H.
Sakai
, and
S.
Torquato
, “
Statistical mechanical models with effective potentials: Definitions, applications, and thermodynamic consequences
,”
J. Chem. Phys.
117
,
288
296
(
2002
).
38.
A. A.
Louis
, “
Beware of density dependent pair potentials
,”
J. Phys.: Condens. Matter
14
,
9187
9206
(
2002
).
39.
M. E.
Johnson
,
T.
Head-Gordon
, and
A. A.
Louis
, “
Representability problems for coarse-grained water potentials
,”
J. Chem. Phys.
126
,
144509
(
2007
).
40.
H.
Wang
,
C.
Junghans
, and
K.
Kremer
, “
Comparative atomistic and coarse-grained study of water: What do we lose by coarse-graining?
,”
Eur. Phys. J. E: Soft Matter Biol. Phys.
28
,
221
229
(
2009
).
41.
A. J.
Clark
,
J.
McCarty
,
I. Y.
Lyubimov
, and
M. G.
Guenza
, “
Thermodynamic consistency in variable-level coarse graining of polymeric liquids
,”
Phys. Rev. Lett.
109
,
168301
(
2012
).
42.
J.
McCarty
,
A. J.
Clark
,
J.
Copperman
, and
M. G.
Guenza
, “
An analytical coarse-graining method which preserves the free energy, structural correlations, and thermodynamic state of polymer melts from the atomistic to the mesoscale
,”
J. Chem. Phys.
140
,
204913
(
2014
).
43.
M.
Guenza
, “
Thermodynamic consistency and other challenges in coarse-graining models
,”
Eur. Phys. J.: Spec. Top.
224
,
2177
2191
(
2015
).
44.
G.
D’Adamo
,
A.
Pelissetto
, and
C.
Pierleoni
, “
Predicting the thermodynamics by using state-dependent interactions
,”
J. Chem. Phys.
138
,
234107
(
2013
).
45.
J.
Lu
,
Y.
Qiu
,
R.
Baron
, and
V.
Molinero
, “
Coarse-graining of TIP4P/2005, TIP4P-Ew, SPC/E, and TIP3P to monatomic anisotropic water models using relative entropy minimization
,”
J. Chem. Theory Comput.
10
,
4104
4120
(
2014
).
46.
A.
Khot
,
S. B.
Shiring
, and
B. M.
Savoie
, “
Evidence of information limitations in coarse-grained models
,”
J. Chem. Phys.
151
,
244105
(
2019
).
47.
A.
Das
and
H. C.
Andersen
, “
The multiscale coarse-graining method. V. Isothermal-isobaric ensemble
,”
J. Chem. Phys.
132
,
164106
(
2010
).
48.
N. J. H.
Dunn
and
W. G.
Noid
, “
Bottom-up coarse-grained models that accurately describe the structure, pressure, and compressibility of molecular liquids
,”
J. Chem. Phys.
143
,
243148
(
2015
).
49.
N. J. H.
Dunn
and
W. G.
Noid
, “
Bottom-up coarse-grained models with predictive accuracy and transferability for both structural and thermodynamic properties of heptane-toluene mixtures
,”
J. Chem. Phys.
144
,
204124
(
2016
).
50.
A.
Lyubartsev
,
A.
Mirzoev
,
L. J.
Chen
, and
A.
Laaksonen
, “
Systematic coarse-graining of molecular models by the Newton inversion method
,”
Faraday Discuss.
144
,
43
56
(
2010
).
51.
M.
Jochum
,
D.
Andrienko
,
K.
Kremer
, and
C.
Peter
, “
Structure-based coarse-graining in liquid slabs
,”
J. Chem. Phys.
137
,
064102
(
2012
).
52.
C.
Dalgicdir
,
O.
Sensoy
,
C.
Peter
, and
M.
Sayar
, “
A transferable coarse-grained model for diphenylalanine: How to represent an environment driven conformational transition
,”
J. Chem. Phys.
139
,
234115
(
2013
).
53.
A.
Ghoufi
,
P.
Malfreyt
, and
D. J.
Tildesley
, “
Computer modelling of the surface tension of the gas–liquid and liquid–liquid interface
,”
Chem. Soc. Rev.
45
,
1387
1409
(
2016
).
54.
L.
Larini
,
L.
Lu
, and
G. A.
Voth
, “
The multiscale coarse-graining method. VI. Implementation of three-body coarse-grained potentials
,”
J. Chem. Phys.
132
,
164107
(
2010
).
55.
A.
Das
and
H. C.
Andersen
, “
The multiscale coarse-graining method. IX. A general method for construction of three body coarse-grained force fields
,”
J. Chem. Phys.
136
,
194114
(
2012
).
56.
C.
Scherer
and
D.
Andrienko
, “
Understanding three-body contributions to coarse-grained force fields
,”
Phys. Chem. Chem. Phys.
20
,
22387
22394
(
2018
).
57.
J. F.
Dama
,
A. V.
Sinitskiy
,
M.
McCullagh
,
J.
Weare
,
B.
Roux
,
A. R.
Dinner
, and
G. A.
Voth
, “
The theory of ultra-coarse-graining. 1. General principles
,”
J. Chem. Theory Comput.
9
,
2466
2480
(
2013
).
58.
J.
Jin
,
Y.
Han
, and
G. A.
Voth
, “
Ultra-coarse-grained liquid state models with implicit hydrogen bonding
,”
J. Chem. Theory Comput.
14
,
6159
6174
(
2018
).
59.
T.
Bereau
and
J. F.
Rudzinski
, “
Accurate structure-based coarse graining leads to consistent barrier-crossing dynamics
,”
Phys. Rev. Lett.
121
,
256002
(
2018
).
60.
J. F.
Rudzinski
and
T.
Bereau
, “
Coarse-grained conformational surface hopping: Methodology and transferability
,”
J. Chem. Phys.
153
,
214110
(
2020
).
61.
T.
Lemke
and
C.
Peter
, “
Neural network based prediction of conformational free energies—A new route toward coarse-grained simulation models
,”
J. Chem. Theory Comput.
13
,
6213
6221
(
2017
).
62.
L.
Zhang
,
J.
Han
,
H.
Wang
,
R.
Car
, and
W.
E
, “
DeePCG: Constructing coarse-grained models via deep neural networks
,”
J. Chem. Phys.
149
,
034101
(
2018
).
63.
J.
Wang
,
S.
Olsson
,
C.
Wehmeyer
,
A.
Pérez
,
N. E.
Charron
,
G.
De Fabritiis
,
F.
Noé
, and
C.
Clementi
, “
Machine learning of coarse-grained molecular dynamics force fields
,”
ACS Cent. Sci.
5
,
755
767
(
2019
).
64.
B. E.
Husic
,
N. E.
Charron
,
D.
Lemm
,
J.
Wang
,
A.
Pérez
,
M.
Majewski
,
A.
Krämer
,
Y.
Chen
,
S.
Olsson
,
G.
de Fabritiis
,
F.
Noé
, and
C.
Clementi
, “
Coarse graining molecular dynamics with graph neural networks
,”
J. Chem. Phys.
153
,
194101
(
2020
); arXiv:2007.11412.
65.
I.
Pagonabarraga
and
D.
Frenkel
, “
Dissipative particle dynamics for interacting systems
,”
J. Chem. Phys.
115
,
5015
5026
(
2001
).
66.
P.
Español
and
P. B.
Warren
, “
Perspective: Dissipative particle dynamics
,”
J. Chem. Phys.
146
,
150901
(
2017
).
67.
I.
Pagonabarraga
and
D.
Frenkel
, “
Non-ideal DPD fluids
,”
Mol. Simul.
25
,
167
175
(
2000
).
68.
S.
Merabia
and
I.
Pagonabarraga
, “
Density dependent potentials: Structure and thermodynamics
,”
J. Chem. Phys.
127
,
054903
(
2007
).
69.
R.
Evans
, “
Nature of the liquid-vapor interface and other topics in the statistical-mechanics of nonuniform, classical fluids
,”
Adv. Phys.
28
,
143
200
(
1979
).
70.
R. D.
Groot
and
P. B.
Warren
, “
Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation
,”
J. Chem. Phys.
107
,
4423
4435
(
1997
).
71.
P. B.
Warren
, “
Hydrodynamic bubble coarsening in off-critical vapor-liquid phase separation
,”
Phys. Rev. Lett.
87
,
225702
(
2001
).
72.
S. Y.
Trofimov
,
E. L. F.
Nies
, and
M. A. J.
Michels
, “
Thermodynamic consistency in dissipative particle dynamics simulations of strongly nonideal liquids and liquid mixtures
,”
J. Chem. Phys.
117
,
9383
9394
(
2002
).
73.
P. B.
Warren
, “
Vapor-liquid coexistence in many-body dissipative particle dynamics
,”
Phys. Rev. E
68
,
066702
(
2003
).
74.
A.
Ghoufi
and
P.
Malfreyt
, “
Calculation of the surface tension from multibody dissipative particle dynamics and Monte Carlo methods
,”
Phys. Rev. E
82
,
016706
(
2010
).
75.
A.
Ghoufi
and
P.
Malfreyt
, “
Mesoscale modeling of the water liquid-vapor interface: A surface tension calculation
,”
Phys. Rev. E
83
,
051601
(
2011
).
76.
M.
Hömberg
and
M.
Müller
, “
Main phase transition in lipid bilayers: Phase coexistence and line tension in a soft, solvent-free, coarse-grained model
,”
J. Chem. Phys.
132
,
155104
(
2010
).
77.
G. A.
Papoian
,
J.
Ulander
,
M. P.
Eastwood
,
Z.
Luthey-Schulten
, and
P. G.
Wolynes
, “
Water in protein structure prediction
,”
Proc. Natl. Acad. Sci. U. S. A.
101
,
3352
3357
(
2004
).
78.
W.
Zheng
,
N. P.
Schafer
,
A.
Davtyan
,
G. A.
Papoian
, and
P. G.
Wolynes
, “
Predictive energy landscapes for protein–protein association
,”
Proc. Natl. Acad. Sci. U. S. A.
109
,
19244
19249
(
2012
).
79.
E. C.
Allen
and
G. C.
Rutledge
, “
A novel algorithm for creating coarse-grained, density dependent implicit solvent models
,”
J. Chem. Phys.
128
,
154115
(
2008
).
80.
E. C.
Allen
and
G. C.
Rutledge
, “
Evaluating the transferability of coarse-grained, density-dependent implicit solvent models to mixtures and chains
,”
J. Chem. Phys.
130
,
034904
(
2009
).
81.
T.
Sanyal
and
M. S.
Shell
, “
Transferable coarse-grained models of liquid–liquid equilibrium using local density potentials optimized with the relative entropy
,”
J. Phys. Chem. B
122
,
5678
5693
(
2018
).
82.
D.
Rosenberger
,
T.
Sanyal
,
M. S.
Shell
, and
N. F. A.
van der Vegt
, “
Transferability of local density-assisted implicit solvation models for homogeneous fluid mixtures
,”
J. Chem. Theory Comput.
15
,
2881
2895
(
2019
).
83.
T.
Sanyal
and
M. S.
Shell
, “
Coarse-grained models using local-density potentials optimized with the relative entropy: Application to implicit solvation
,”
J. Chem. Phys.
145
,
034109
(
2016
).
84.
N.
Shahidi
,
A.
Chazirakis
,
V.
Harmandaris
, and
M.
Doxastakis
, “
Coarse-graining of polyisoprene melts using inverse Monte Carlo and local density potentials
,”
J. Chem. Phys.
152
,
124902
(
2020
).
85.
F.
Berressem
,
C.
Scherer
,
D.
Andrienko
, and
A.
Nikoubashman
, “
Ultra-coarse-graining of homopolymers in inhomogeneous systems
,”
J. Phys.: Condens. Matter
33
,
254002
(
2021
).
86.
J. D.
Moore
,
B. C.
Barnes
,
S.
Izvekov
,
M.
Lísal
,
M. S.
Sellers
,
D. E.
Taylor
, and
J. K.
Brennan
, “
A coarse-grain force field for RDX: Density dependent and energy conserving
,”
J. Chem. Phys.
144
,
104501
(
2016
).
87.
V.
Agrawal
,
P.
Peralta
,
Y.
Li
, and
J.
Oswald
, “
A pressure-transferable coarse-grained potential for modeling the shock Hugoniot of polyethylene
,”
J. Chem. Phys.
145
,
104903
(
2016
).
88.
J. W.
Wagner
,
T.
Dannenhoffer-Lafage
,
J.
Jin
, and
G. A.
Voth
, “
Extending the range and physical accuracy of coarse-grained models: Order parameter dependent interactions
,”
J. Chem. Phys.
147
,
044113
(
2017
).
89.
M. R.
DeLyser
and
W. G.
Noid
, “
Extending pressure-matching to inhomogeneous systems via local-density potentials
,”
J. Chem. Phys.
147
,
134111
(
2017
).
90.
M. R.
DeLyser
and
W. G.
Noid
, “
Analysis of local density potentials
,”
J. Chem. Phys.
151
,
224106
(
2019
).
91.
M.
DeLyser
and
W. G.
Noid
, “
Bottom-up coarse-grained models for external fields and interfaces
,”
J. Chem. Phys.
153
,
224103
(
2020
).
92.
P. M.
Chaikin
and
T. C.
Lubensky
,
Principles of Condensed Matter Physics
(
Cambridge University Press
,
2012
).
93.
A. P.
Thompson
,
S. J.
Plimpton
, and
W.
Mattson
, “
General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions
,”
J. Chem. Phys.
131
,
154107
(
2009
).
94.
M. E.
Tuckerman
,
Statistical Mechanics: Theory and Molecular Simulation
(
Oxford University Press
,
Oxford, UK
,
2013
).
95.
P.
Español
and
P.
Warren
, “
Statistical mechanics of dissipative particle dynamics
,”
Europhys. Lett.
30
,
191
(
1995
).
96.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
97.
G. J.
Martyna
,
D. J.
Tobias
, and
M. L.
Klein
, “
Constant pressure molecular dynamics algorithms
,”
J. Chem. Phys.
101
,
4177
4189
(
1994
).
98.
G. J.
Martyna
,
M. E.
Tuckerman
,
D. J.
Tobias
, and
M. L.
Klein
, “
Explicit reversible integrators for extended systems dynamics
,”
Mol. Phys.
87
,
1117
1157
(
1996
).
99.
S. Y.
Trofimov
,
E. L. F.
Nies
, and
M. A. J.
Michels
, “
Constant-pressure simulations with dissipative particle dynamics
,”
J. Chem. Phys.
123
,
144102
(
2005
).
100.
H. C.
Andersen
, “
Molecular dynamics simulations at constant pressure and/or temperature
,”
J. Chem. Phys.
72
,
2384
2393
(
1980
).
101.
M. E.
Fisher
and
D.
Ruelle
, “
The stability of many-particle systems
,”
J. Math. Phys.
7
,
260
270
(
1966
).
102.
D.
Ruelle
,
Statistical Mechanics: Rigorous Results
(
W. A. Benjamin, Inc.
,
1969
).
103.
S.
Izvekov
and
G. A.
Voth
, “
A multiscale coarse-graining method for biomolecular systems
,”
J. Phys. Chem. B
109
,
2469
2473
(
2005
).
104.
S.
Izvekov
and
G. A.
Voth
, “
Multiscale coarse graining of liquid-state systems
,”
J. Chem. Phys.
123
,
134105
(
2005
).
105.
L.
Lu
and
G. A.
Voth
, “
The multiscale coarse-graining method
,”
Adv. Chem. Phys.
149
,
47
81
(
2012
).
106.
W. G.
Noid
,
J.-W.
Chu
,
G. S.
Ayton
,
V.
Krishna
,
S.
Izvekov
,
G. A.
Voth
,
A.
Das
, and
H. C.
Andersen
, “
The multiscale coarse-graining method. I. A rigorous bridge between atomistic and coarse-grained models
,”
J. Chem. Phys.
128
,
244114
(
2008
).
107.
W. G.
Noid
,
P.
Liu
,
Y.
Wang
,
J.-W.
Chu
,
G. S.
Ayton
,
S.
Izvekov
,
H. C.
Andersen
, and
G. A.
Voth
, “
The multiscale coarse-graining method. II. Numerical implementation for molecular coarse-grained models
,”
J. Chem. Phys.
128
,
244115
(
2008
).
108.
J. F.
Rudzinski
and
W. G.
Noid
, “
A generalized-Yvon-Born-Green method for coarse-grained modeling
,”
Eur. Phys. J.: Spec. Top.
224
,
2193
2216
(
2015
).
109.
L.
Lu
,
S.
Izvekov
,
A.
Das
,
H. C.
Andersen
, and
G. A.
Voth
, “
Efficient, regularized, and scalable algorithms for multiscale coarse-graining
,”
J. Chem. Theory Comput.
6
,
954
965
(
2010
).
110.
M. J.
Abraham
,
T.
Murtola
,
R.
Schulz
,
S.
Páll
,
J. C.
Smith
,
B.
Hess
, and
E.
Lindahl
, “
GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers
,”
SoftwareX
1-2
,
19
25
(
2015
).
111.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
, “
Canonical sampling through velocity rescaling
,”
J. Chem. Phys.
126
,
014101
(
2007
).
112.
W. L.
Jorgensen
,
D. S.
Maxwell
, and
J.
Tirado-Rives
, “
Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids
,”
J. Am. Chem. Soc.
118
,
11225
11236
(
1996
).
113.
T.
Darden
,
D.
York
, and
L.
Pedersen
, “
Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems
,”
J. Chem. Phys.
98
,
10089
10092
(
1993
).
114.
M.
Parrinello
and
A.
Rahman
, “
Crystal structure and pair potentials: A molecular-dynamics study
,”
Phys. Rev. Lett.
45
,
1196
1199
(
1980
).
115.
N. J. H.
Dunn
,
K. M.
Lebold
,
M. R.
DeLyser
,
J. F.
Rudzinski
, and
W. G.
Noid
, “
BOCS: Bottom-up open-source coarse-graining software
,”
J. Phys. Chem. B
122
,
3363
3377
(
2018
).
116.
S.
Nosé
, “
A molecular dynamics method for simulations in the canonical ensemble
,”
Mol. Phys.
52
,
255
268
(
1984
).
117.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
,
1695
1697
(
1985
).
118.
G. J.
Martyna
,
M. L.
Klein
, and
M.
Tuckerman
, “
Nosé–Hoover chains: The canonical ensemble via continuous dynamics
,”
J. Chem. Phys.
97
,
2635
2643
(
1992
).
119.
K.
Shen
,
N.
Sherck
,
M.
Nguyen
,
B.
Yoo
,
S.
Köhler
,
J.
Speros
,
K. T.
Delaney
,
G. H.
Fredrickson
, and
M. S.
Shell
, “
Learning composition-transferable coarse-grained models: Designing external potential ensembles to maximize thermodynamic information
,”
J. Chem. Phys.
153
,
154116
(
2020
).
120.
M. C.
Villet
and
G. H.
Fredrickson
, “
Numerical coarse-graining of fluid field theories
,”
J. Chem. Phys.
132
,
034109
(
2010
).
121.
J. W.
Wagner
,
J. F.
Dama
,
A. E. P.
Durumeric
, and
G. A.
Voth
, “
On the representability problem and the physical meaning of coarse-grained models
,”
J. Chem. Phys.
145
,
044108
(
2016
).
122.
T.
Dannenhoffer-Lafage
,
J. W.
Wagner
,
A. E. P.
Durumeric
, and
G. A.
Voth
, “
Compatible observable decompositions for coarse-grained representations of real molecular systems
,”
J. Chem. Phys.
151
,
134115
(
2019
).
123.
K. M.
Lebold
and
W. G.
Noid
, “
Dual approach for effective potentials that accurately model structure and energetics
,”
J. Chem. Phys.
150
,
234107
(
2019
).
124.
K. M.
Lebold
and
W. G.
Noid
, “
Dual-potential approach for coarse-grained implicit solvent models with accurate, internally consistent energetics and predictive transferability
,”
J. Chem. Phys.
151
,
164113
(
2019
).
125.
J.
Towns
,
T.
Cockerill
,
M.
Dahan
,
I.
Foster
,
K.
Gaither
,
A.
Grimshaw
,
V.
Hazlewood
,
S.
Lathrop
,
D.
Lifka
,
G. D.
Peterson
,
R.
Roskies
,
J. R.
Scott
, and
N.
Wilkins-Diehr
, “
XSEDE: Accelerating scientific discovery
,”
Comput. Sci. Eng.
16
,
62
74
(
2014
).