Coarse-grained (CG) simulation models of condensed-phase systems can be derived with well-established methods that perform coarse-graining in space and provide an effective Hamiltonian with which some of the structural and thermodynamic properties of the underlying fine-grained (FG) reference system can be represented. Coarse-graining in time potentially provides CG models that furthermore represent dynamic properties. However, systematic efforts in this direction have so far been limited, especially for moderately coarse-grained, chemistry-specific systems with complicated conservative interactions. With the aim of representing structural, thermodynamic, and dynamic properties in CG simulations of multi-component molecular systems, we investigated a recently introduced method in which the force on a CG particle originates from conservative interactions with surrounding particles and non-Markovian dissipative interactions, the latter introduced by means of a colored-noise thermostat. We examined two different methods to derive isotropic memory kernels required for integrating the corresponding generalized Langevin equation (GLE) of motion, based on the orthogonal dynamics of the FG forces and on an iterative optimization scheme. As a proof of concept, we coarse-grain single-component molecular liquids (cyclohexane, tetrachloromethane) and ideal and non-ideal binary mixtures of cyclohexane/tetrachloromethane and ethanol/tetrachloromethane, respectively. We find that for all systems, the FG single particle velocity auto-correlation functions and, consequently, both the short time and long time diffusion coefficients can be quantitatively reproduced with the CG-GLE models. We furthermore demonstrate that the present GLE-approach leads to an improved description of the rate with which the spatial correlations decay, which is artificially accelerated in the absence of dissipation.

Soft matter systems, often encountered in chemistry and biology, exhibit various intriguing properties at distinct length and time scales. To investigate the underlying microscopic origins with computer simulations, systematic coarse-graining strategies are often employed in which several atoms are merged to define effective interaction sites whose interactions preserve information on the chemical specificity of the parent atomistic, or fine-grained (FG), system. This approach to coarse-graining in configuration space, also referred to as structure-based coarse-graining, has been an active field in soft matter research over the last two decades.1–6 Various methods have been proposed to systematically parameterize the molecular interactions in CG models,7–13 all with the aim of representing the structural and/or thermodynamic properties of the FG system at the CG level. Applications include, among others, polymeric1,14–17 and biomolecular4,5,18 systems. “Dynamic coarse-graining,” i.e., coarse-graining in time, by contrast, is still in its infancy and additionally targets representing (equilibrium) dynamics, which is ill-represented in molecular dynamics (MD) simulations of CG models derived with the above systematic coarse-graining methods.19,20 The underlying theory follows the Mori–Zwanzig (MZ) formalism,21–23 where a projection operator is employed to coarse-grain time reversible Hamiltonian dynamics of the reference FG system to time irreversible stochastic dynamics at the CG level. The resulting CG equation of motion (EoM) resembles a generalized Langevin equation (GLE), where the dynamics of the CG particles depend on their dynamical history through a non-Markovian memory kernel that distinctly couples to the various time scales in the system.

Despite a rigorous background in statistical mechanics, the application of dynamic coarse-graining as a computational strategy is so far limited.20 While there have been some efforts toward deriving exact CG EoMs,24 their non-Markovian nature makes it difficult to utilize in computer simulations. One way to evade this problem is to assume a complete time scale separation between the degrees of freedom present in the CG model and the faster (FG) ones that have been removed. This assumption renders the CG EoM Markovian. The dissipative interactions may then be modeled using bottom-up parameterized Langevin or dissipative particle dynamics (DPD)25 thermostats, which, unlike the GLE thermostat, use time-independent friction that equally couples to the various time scales in the system. Such bottom-up parameterized Markovian models have been developed for both model and molecular systems with varying success.24,26–29 Technical difficulties related to incomplete time scale separation may, however, complicate the derivation of time-independent (Markovian) frictions whose application may, moreover, lead to an over-prediction of the diffusion coefficient for the systems studied here, as will be further discussed below. Alternatively, Markovian thermostats can always be parameterized top-down such that the long time diffusion coefficients of the CG and FG models match by construction. However, in such a case the connection to the underlying FG system is lost and, thus, the CG model is not systematic.

Recent works30,31 have shown that relaxation mechanisms at short time scales can non-trivially affect the long time diffusive dynamics of a system. Given that the overall dynamics of the system is governed by a complex interplay of various relaxation time scales, Markovian friction may not always be adequate, which necessitates the use of non-Markovian models in CG molecular simulations. However, most of the studies along these lines have so far focused on simple model systems.32–37 A few studies that did focus on molecular (chemically-specific) systems reported very encouraging results but did not succeed in quantitatively reproducing the reference dynamical quantities38,39 (except for Ref. 40). Accurate computational modeling of these types of systems requires further development and testing of methodologies aiming to provide guidelines for parametrizing time-dependent dissipative interactions (non-Markovian memory kernels) that are consistent with the concurrent conservative interactions present in the various specific systems encountered in complex soft matter.

Herein, we coarse-grain the dynamics of molecular liquids and their ideal/non-ideal mixtures by extending the standard Newtonian EoM in CG MD simulation with a simple isotropic GLE thermostat. We use two recently proposed bottom-up approaches to parameterize such thermostats: a direct parameterization method30,41 and an iterative optimization method.40,41 While the systems considered are molecular liquids that can be well studied in FG simulations without the need for structural/dynamic coarse-graining, they do serve as proof of concept. The ultimate goal is for the GLE thermostat to “fix” the dynamical properties of computationally tractable CG models of complex multicomponent molecular systems. Dynamic coarse-graining on molecular mixtures, in which the dynamics of the participating species are coupled and composition dependent, provides a reasonable testing ground for methodological development. Moreover, non-ideal molecular mixtures bring us one step closer to investigating complex molecular systems such as ionic liquids,42 where local inhomogeneity gives rise to structural correlations that are beyond the capabilities of FG simulations to probe. Our results suggest that bottom-up parameterized GLE thermostats can efficiently correct single particle dynamics not only in pure molecular liquids but also in ideal and non-ideal binary mixtures, while quantitative agreement between FG and CG dynamical observables can be achieved by using the iterative scheme. Moreover, while the GLE thermostat exerts friction at the single bead level, it can reasonably reproduce the temporal decay of structural correlations in molecular mixtures that are otherwise drastically accelerated in CG MD simulations. These encouraging observations suggest the possibility of using simple and tractable GLE-based CG models toward dynamically consistent CG modeling of complex soft matter systems.

We performed atomistic simulations of single-component molecular liquids and ideal and non-ideal binary mixtures. The two single-component systems consisted of 1000 molecules of tetrachloromethane (CCl4; henceforth denoted as TCM) and cyclohexane (C6H12; henceforth denoted as CHX). The three ideal (homogeneous) CHX-TCM mixtures consisted of 2000 molecules in the number ratios 25:75, 50:50, and 75:25, which are known to show minor deviations from ideality.43 The non-ideal system was a mixture of ethanol (C2H5OH; henceforth denoted as ETH) and TCM in a ratio of 75:25. Binary ETH-TCM mixtures are known to be non-ideal,43 and the composition here was chosen such that the CG conservative interaction can be extracted using systematic coarse-graining. All the systems were modeled using the OPLS/AA force field.44 The simulations were performed using the Gromacs package45 using an integration time step of 0.5 fs. The LINCS algorithm46 was used to constrain bonds involving hydrogen atoms. A cut-off distance of 1.4 nm and geometric combination rules were used for the van der Waals interactions. The long range electrostatic interactions were computed using the PME algorithm47 with a real space cut-off of 1.4 nm, a grid spacing of 0.12 nm, and an interpolation order of 4. Long range dispersion corrections were applied to calculate the energy and pressure of the system. After the initial energy minimization using the steepest descent method, a 0.5 ns MD run was performed under constant NPT conditions, followed by a 0.5 ns MD run under constant NVT conditions. Five cycles of such NPT-NVT runs were performed to properly equilibrate the system. During these runs, the systems were maintained at a temperature of 300 K and a pressure of 1 bar using the velocity-rescale thermostat48 and the Berendsen barostat49 with time constants of 1 and 2 ps, respectively. This procedure was followed by a 3 ns production NPT simulation using the Nośe–Hoover thermostat50 and the Parrinello–Rahman barostat51 with time constants of 2 and 4 ps, respectively. Subsequently, a 10 ns production run was performed under NVT conditions using the Nośe–Hoover thermostat, which was used to calculate the radial distribution functions (RDFs) using the molecular center of mass (CoM) as reference sites and a spatial resolution of 0.01 nm. Thereafter, a 100 ps NVT run was performed in which the forces and velocities were stored at an interval of 1 fs. This short trajectory was used to parameterize the GLE thermostats, as discussed in Sec. II B.

We employ a previously studied EoM in which each CG particle i is coupled to an isotropic GLE-based thermostat according to30,40,41
Fi(t)=FiC(t)mi0tK̃(ts)Vi(s)ds+F̃iR(t),
(1)
where Fi(t), FiC(t), F̃iR(t), mi, and Vi(s) are the total force, the force due to the conservative interactions and the random force acting on the CG particle i at time t, and its mass and velocity, respectively. K̃(t) denotes the single particle memory kernel, which is related to the random force through the relation F̃iR(t)F̃iR(0)=3mikBTK̃(t), where kB is the Boltzmann constant and T is the temperature, ensuring canonical sampling. Parameterizing the above CG EoM, therefore, boils down to parameterizing two quantities: the conservative interactions and the memory kernels. It should be noted that, in general, the memory kernel can be position dependent.52,53 However, here we employ a computationally tractable model using memory kernels with no explicit position dependence.

The conservative interactions were optimized so as to reproduce the FG RDFs using the iterative Boltzmann inversion (IBI) method,8 as implemented in the VOTCA package.54,55 For pure liquids, the RDFs converged around 50 iterations (Fig. S1). For the CHX-TCM mixtures, the three intra- and inter-species RDFs were found to converge within 150 iterations (Figs. S2 and S3). Figure S2 shows snapshots of the two components in the 75:25 mixture. Their homogeneous spatial distribution indicates the ideal nature of this mixture. The three RDFs for the non-ideal ETH-TCM mixture were found to converge around 230 iterations (Fig. S4). Although the RDFs were visually close, the agreement was farther from that observed in the case of ideal mixtures (compare against Figs. S2 and S3). Figure S4 also shows representative snapshots of the ETH and TCM components, which clearly indicate non-ideal mixing behavior. In all cases, the interaction potentials at the corresponding converged iteration were taken as the CG interaction potentials (Fig. S1 through Fig. S4).

The memory kernels K̃(t) of the GLE thermostats were parameterized using: (1) a direct parameterization method30 and (2) an iterative optimization method.40 The direct parameterization method provides insight regarding the amplitudes and time scales of fluctuating forces that contribute to K̃(t) and the effective particle friction at the CG level. The iterative method optimizes K̃(t) such that its application in CG GLE simulations leads to a quantitative representation of the reference FG velocity auto-correlation function (VACF). These two methods have been successfully applied in studies of a model star polymer melt,30 of model polymer-colloid41 systems, and more recently, of liquid water.40 The exact methodology behind the two schemes has been extensively discussed in the original references and is only briefly summarized below.

The direct parameterization method is based on the MZ theory21–23 and evaluates the single particle memory kernel K(t) from the Q-projected time correlation functions of forces, estimated at the mapping centers in the CG representation of a FG trajectory, following the fluctuation–dissipation theorem. For a single-component liquid, K(t) is related to its self-diffusion coefficient D according to 0K(t)dt=kBT/miD. The direct parameterization method splits K(t) into contributions originating from fluctuating FG DoFs that were removed in the CG model, the remaining CG DoFs, and the cross-correlation between the two.30 This splitting allows for a reasonable a priori estimate of the contribution from the conservative interactions to the total single-particle memory kernel K(t) in CG simulations, which in turn allows us to estimate an appropriate memory kernel K̃(t) for the GLE thermostat parameterization.30,41 A representative set of memory kernels with all relevant contributions is provided in Fig. S5 in the supplementary material. All the time correlation functions were evaluated using the backward orthogonal dynamics (BOD) method,56 which provides a generic framework to compute time correlation functions of arbitrary Q-projected dynamical observables.

The iterative approach follows the iterative optimization of memory kernels (IOMK) method40 to quantitatively reproduce the FG VACF. The iterative scheme is given as
G̃i+1(t)=Gtarget(t)a(t)GiCG(t)a(t)G̃i(t),
(2)
where G(t)=0tK(s)ds denotes the integrated single particle memory kernel, which is related to the VACF through the Volterra-type equation,40,57
CVV(t)CVV(0)=0tG(ts)CVV(s)ds
(3)
and, therefore, can be numerically evaluated from the VACF (CVV). Here, Gtarget indicates the integrated effective single-particle memory kernel evaluated from the VACF (calculated using the mapped FG trajectory) using the numerical procedure discussed in Ref. 57, which is based on discretization in the time domain. GiCG indicates the integrated single particle memory calculated in the CG GLE simulation, and G̃i indicates the integrated single particle memory exerted by the GLE thermostat at the i-th iteration. The term a(t) denotes the effectively integrated memory in the absence of a GLE thermostat, i.e., a priori calculated in a CG MD simulation. Equation (2) can be interpreted as a quasi-Newton optimization scheme that converges within a few iterations for single-component systems.40 

Herein, GLE thermostats were implemented based on the auxiliary variables approach,58,59 which employs a set of auxiliary momenta to render the non-Markovian EoM Markovian in extended phase space. The integrated memory kernels obtained with the two approaches were fitted to a combination of six damped harmonic oscillators. The resulting coupling matrix with 12 auxiliary momenta was used in the GLE thermostat to fix the temperature at 300 K in CG GLE simulations. The non-bonded interactions, optimized using the IBI method,8 were used as tabulated potentials with a cutoff of 2 nm in both CG GLE and CG MD simulations. We note that an alternative, recently reported method for structural coarse-graining based on a Hypernetted-Chain-Gauss-Newton (HNCGN) approach allows for significantly shorter interaction cutoffs.13 The CG simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package60,61 under NVT conditions using a time step of 1 fs and the velocity-Verlet time integration algorithm. The number of CG beads and the corresponding box volumes in CG simulations were kept the same as their FG NVT MD counterparts. In the case of binary mixtures, two memory kernels were parameterized, one for each component, and two thermostats were used in the CG GLE simulation.

For the systems studied here, the VACFs of the CG CoMs and the corresponding single particle memory kernels K̃(t) decay over similar time scales [shown in Fig. 1(a), Figs. S6(a), S5(a), and S5(c) for the two single-component liquids], indicating non-Markovian behavior. Despite the lack of time scale separation, we derived bottom-up parameterized Markovian frictions to check if the long time diffusive dynamics could be represented by using a Langevin thermostat in CG simulations. We extracted the time independent friction from K̃(t) (evaluated using the BOD method) from the long time asymptotic value of the fitting function to the integrated single particle memory kernel G̃(t)=0tK̃(s)ds, i.e., γ̃=limtG̃(t). Integrated kernels G̃(t) for pure liquids are reported in Figs. S5(b) and S5(d) in the supplementary material. CG simulations using such bottom-up parameterized Langevin thermostats (hereafter denoted as CG LE simulations) were performed for the two molecular liquids and their three ideal mixtures using the LAMMPS package.60,61

FIG. 1.

(a) VACFs and (b) their running integrals, D(t), computed with FG MD, CG MD, CG LE (BOD), and CG GLE (BOD) simulations of the single-component CHX liquid. The inset shows a zoomed in view of the VACF minima.

FIG. 1.

(a) VACFs and (b) their running integrals, D(t), computed with FG MD, CG MD, CG LE (BOD), and CG GLE (BOD) simulations of the single-component CHX liquid. The inset shows a zoomed in view of the VACF minima.

Close modal

1. Dynamics in CG simulations: The direct parameterization approach

Figure 1(a) compares the VACFs calculated with FG MD, CG MD, CG LE (BOD), and CG GLE (BOD) simulations of the single-component CHX system. The corresponding running integral 0tdsV(s)V(0) [i.e., 1/30tdsV(s)V(0)] indicates the time dependent diffusion coefficient D(t) and is shown in Fig. 1(b). The corresponding plots for the pure TCM system are provided in Fig. S6. The absence of friction in the CG MD simulation leads to a large deviation between the VACFs obtained with the FG and CG MD models, and consequently in D(t). As compared to CG MD simulations, CG LE simulations do not significantly improve the VACFs but match the D(t) of the FG systems at large t reasonably well. The long time limiting values of D(t) predicted by the CG LE (BOD) model slightly overestimate the target FG MD value. This is caused by the negative and slowly varying long time tail of K̃ [Figs. S5(a)–S5(c)], which reduces the friction in the LE model but couples weakly to velocity fluctuations in the GLE model. By contrast, CG GLE (BOD) simulations not only reproduce the initial decay of the FG VACF but also its overall profile, including the first minimum. Subsequently, D(t) is very close to the FG MD values in the CG GLE (BOD) simulation.

Next, we study the dynamics of the two components in CHX-TCM ideal mixtures at number ratios of 25:75, 50:50, and 75:25. The VACFs and the corresponding D(t) for the three mixtures are provided in Fig. S7 through Fig. S9 in the supplementary material. Here, we use two directly parameterized isotropic GLE thermostats in the CG GLE (BOD) simulations and their time independent Markovian counterparts in the CG LE (BOD) simulations. As the dynamics of the participating species in a mixture are coupled, it remains to be verified if two thermostats that are independently coupled to each component of the mixture lead to a correct representation of dynamics in the CG simulation. Interestingly, we observe very good agreement between the VACFs [Figs. S7–S9(a) and S7–S9(c)] and D(t) [Figs. S7–S9(b) and S7–S9(d)] of the two components in the binary mixtures calculated with FG MD and CG GLE (BOD) simulations. As shown above, CG LE (BOD) simulations do not noticeably improve the VACFs, and D(t) deviates stronger from the FG results for the molecular mixtures as compared to what was observed for single-component liquids [see, for example, Fig. 1(b) and Fig. S6(b) vs Figs. S7–S9(b) and S7–S9(d)].

In Fig. 2(a), we compare the diffusion coefficients D = limtD(t), calculated as the averaged plateau value of D(t), obtained from FG MD, CG MD, and CG GLE (BOD) simulations of the two single-component liquids and their three ideal mixtures. In Fig. 2(b), we compare the results from CG LE (BOD) simulations against CG GLE (BOD) simulations. Data points along the x = y lines, shown as the black dashed lines in the plots, would indicate precise agreement between the FG and CG values. As expected, CG MD simulations highly overestimate D in all cases. By contrast, the values of D predicted with CG GLE (BOD) simulations are close to the FG MD values. Moreover, CG LE (BOD) simulations also lead to an overprediction of D [Fig. 2(b)]. It is interesting to note that Lei and co-workers,26 in their bottom-up parameterized CG LE simulations of clusters of LJ particles, reported the diffusion coefficient to be almost four times smaller than the FG counterpart. Based on our work, we believe that this discrepancy can be attributed to their scheme for parametrizing the LE friction, where they used the auto-correlation function of the forces due to the fast DoFs (i.e., the auto-correlation of instantaneous random forces) as the memory kernel. In our scheme, we include another contribution to our directly parameterized memory kernel, which is the cross-correlation between the forces due to the fast and slow DoFs (i.e., the cross-correlation between the instantaneous random force and the CG conservative force30). This contribution is found to exhibit a slowly decaying negative tail [see Figs. S5(a) and S5(c)]. While using the memory kernel K̃ in CG GLE simulations leads to a satisfactory representation of single-particle dynamics, CG LE simulations performed with the corresponding Markovian frictions overestimate the diffusion coefficients. This happens because the negative and slowly varying long time tail of K̃ [Figs. S5(a) and S5(c)] reduces the effective friction in the LE model. Similar observations were previously reported for a model star polymer system.30 

FIG. 2.

Comparison of diffusion coefficients calculated with (a) FG MD, CG MD, and CG GLE (BOD) simulations and (b) FG MD, CG LE (BOD), and CG GLE (BOD) simulations of single-component liquids (black markers) of CHX and TCM and their mixtures (colored markers) in 25:75, 50:50, and 75:25 number ratios. The x = y dashed black lines would indicate precise agreement between FG and CG values.

FIG. 2.

Comparison of diffusion coefficients calculated with (a) FG MD, CG MD, and CG GLE (BOD) simulations and (b) FG MD, CG LE (BOD), and CG GLE (BOD) simulations of single-component liquids (black markers) of CHX and TCM and their mixtures (colored markers) in 25:75, 50:50, and 75:25 number ratios. The x = y dashed black lines would indicate precise agreement between FG and CG values.

Close modal

On zooming in on the VACFs (Fig. 3), discrepancies are observed between the FG MD (in maroon) and CG GLE (BOD) (in blue) VACFs at small t. While the GLE model quantitatively captures the initial decay, the exact profile of the VACFs at small t is not captured. However, the resulting D(t) matches quite well with those from FG MD simulations due to error compensation (also see Fig. S10). As discussed in detail in the original publications concerning the GLE parametrization method,30,41 this remaining deviation at small t can be largely attributed to the inaccuracy in the CG conservative interactions.

FIG. 3.

VACF (zoomed in) calculated with FG MD and CG GLE (IOMK) simulations of (a) single-component TCM and (b) TCM in the 75:25 CHX-TCM mixture. The VACFs at the zeroth iteration are those calculated with the CG GLE (BOD) simulation (shown in Figs. S6 and S7).

FIG. 3.

VACF (zoomed in) calculated with FG MD and CG GLE (IOMK) simulations of (a) single-component TCM and (b) TCM in the 75:25 CHX-TCM mixture. The VACFs at the zeroth iteration are those calculated with the CG GLE (BOD) simulation (shown in Figs. S6 and S7).

Close modal

Here, the ability of the isotropic GLE model to reproduce single particle dynamics in many-body simulations is noteworthy as (1) so far, studies aiming at reproducing CG dynamics in non-Markovian simulations have mostly focused on model systems rather than on chemically specific systems in which conservative interactions can play a decisive role,20 and (2) existing efforts along these lines have only been moderately successful in terms of quantitatively reproducing FG dynamical observables.38,39 Moreover, efforts to represent the dynamics of molecular mixtures (in which the motions of the different components are coupled) with CG models provide new routes to dynamically consistent CG modeling of more complex molecular systems.

2. Reproducing FG dynamics quantitatively: the iterative optimization approach

While the IOMK method40 can closely reproduce the FG VACFs by construction, it has so far only been applied to single-component systems.40,41 Thus, its applicability to mixtures and complex molecular systems, in general, remains to be verified. As a starting point, we apply this method to coarse-grain molecular liquids and binary mixtures. Although the convergence of the IOMK method is essentially independent of the initial guess, which can be far from optimal, here we use the memory kernels from the direct parametrization (BOD) method as the initial guess. Evidently, these are expected to be very close to the optimal memory kernels as the VACFs are close to the target FG VACFs (Figs. S6–S9). Consequently, as shown in Fig. 3(a), the CG VACF for TCM in the single-component liquid matches quantitatively with its FG counterpart within ten iterations. The plots for the single-component CHX liquid are shown in Fig. S10(a). A similarly good agreement has recently been shown in the case of CG water using the IOMK method.40 Consequently, D(t) for both single-component systems calculated with FG MD and CG GLE (IOMK) simulations is found to be in very good agreement (not shown here).

Next, we applied the IOMK method to coarse-grain dynamics in the binary mixtures of CHX and TCM. It should be reiterated here that in a mixture the motions of the participating species are coupled and, therefore, updating the (integrated) memory kernel of one component at a particular iteration affects the dynamics of the other. Therefore, at any given iteration, we keep one of the two memory kernels fixed and only update the other. Figure 3(b) shows the CG VACF for TCM in the 75:25 mixture, which is in very good agreement with its FG counterparts. Fig. S10(b) in the supplementary material shows the corresponding plot for CHX in a 75:25 mixture, and Figs. S11 and S12 show the plots for the two components in 25:75 and 50:50 mixtures, respectively. We notice that the VACFs match well within ten iterations, during which both the memory kernels get updated five times. Beyond this, updating one kernel does not seem to affect the other, indicating convergence of the iterative scheme.

The observations here are encouraging, as a simple iterative method, which does not even require a good initial guess to begin with,40 is found to quantitatively reproduce the FG VACF in the CG GLE (IOMK) simulation of molecular mixtures. The implementation of the IOMK method in molecular simulation is rather straight-forward, as, in contrast to previously proposed iterative methods,37 it does not require any predetermined parameters. Moreover, it shows better stability and faster convergence.40 The relative simplicity, usability, and robustness have been a few reasons why an iterative scheme such as IBI8 has been a popular approach to bottom-up systematic coarse-graining. In the same way, the IOMK method may serve as an efficient and robust tool for dynamically consistent CG modeling.

3. Beyond single particle dynamics: Time evolution of spatial correlations

The dynamic evolution of the pair correlation between particles in a system can be accessed in terms of the distinct van Hove function,62,
Gd(R,t)=ρgd(R,t)=1NiNjiNδRRj(t)+Ri(0),
(4)
where ρ indicates the particle number density and N indicates the total particle number. The quantity ρgd(R, t)dR equals the number of particles at R at time t, given another particle was at the origin at time zero. Evidently, gd(R, 0) corresponds to the pair correlation function, g(R), and gd(R, t) indicates its temporal evolution, i.e., the decay of the spatial correlations over time. As the dynamics of any single pair of particles depends on the evolution of their surrounding particles in time and space, gd(R, t) encodes the collective structural relaxation in the system. Herein, we employ an isotropic GLE thermostat, which couples to each particle in the system independently. Thus, the model does not include any information on the dynamics of inter-particle correlations. Despite this, we recently found that the decay of pair correlations in liquid water, assessed in terms of gd(R, t) in FG MD, could be reasonably well reproduced using an isotropic GLE thermostat at the CG level.40 Here, we intend to investigate if a similar agreement can be achieved in a molecular mixture where two DoFs are independently coupled to two separate GLE thermostats.

Figures 4(a)4(c) show the distinct van Hove function corresponding to CHX-CHX correlations in the 75:25 binary CHX-TCM mixture calculated with the FG MD, CG MD, and CG GLE (BOD) models. At very small times (t = 0.01 ps), the profile of gd(R, t) is essentially similar to the CHX-CHX g(R) [Fig. S2(a)]. For larger times t, the spatial correlation between the CHX particle pairs decays, and consequently, gd(R, t) loses its structure. In the CG MD simulation, the absence of friction leads to a faster decorrelation between the particles, resulting in a faster decay of the van Hove function as compared to the FG MD case. At time t = 1 ps, where gd(R, t) obtained from the FG MD simulation still shows structural correlations, gd(R, t) obtained from the CG MD simulations is almost featureless. By contrast, gd(R, t) calculated from the CG GLE (BOD) simulation is found to be quite close to the FG MD counterpart. This shows that the time-dependent spatial correlations in the CG model can be represented in good agreement with the FG system when time-dependent friction is introduced through the GLE thermostat. Figure 4(d) shows a close-up of the first peak of gd(R, t) at t = 0.5 ps for FG MD, CG MD, CG LE (BOD), and CG GLE simulations using the direct parametrization (BOD) and iterative optimization (IOMK) schemes. While the peak corresponding to the CG MD simulation has decayed considerably as compared to the FG MD case, it is significantly better represented with the CG LE (BOD) model. The data obtained with the CG GLE (BOD) model indicate a slight over-structuring of the peak [Fig. 4(d)], whereas the CG GLE (IOMK) model shows excellent agreement with the results from the FG MD model.

FIG. 4.

Distinct van Hove function calculated for CHX in the 75:25 CHX-TCM binary mixture obtained with (a) FG MD, (b) CG MD, (c) CG GLE (BOD) simulations, and (d) the first peak of gd(R, t) at t = 0.5 ps calculated with FG MD, CG MD, CG LE, and CG GLE simulations using the direct parametrization (BOD) and iterative optimization (IOMK) schemes.

FIG. 4.

Distinct van Hove function calculated for CHX in the 75:25 CHX-TCM binary mixture obtained with (a) FG MD, (b) CG MD, (c) CG GLE (BOD) simulations, and (d) the first peak of gd(R, t) at t = 0.5 ps calculated with FG MD, CG MD, CG LE, and CG GLE simulations using the direct parametrization (BOD) and iterative optimization (IOMK) schemes.

Close modal

The evolution of gd(R, t) can be comprehensively presented in terms of the temporal decay of the distinct peaks. Fig. S13 shows the decay of the first peak of gd(R, t) at R = 0.62 nm for CHX and at R = 0.57 nm for TCM in the 75:25 CHX-TCM mixture, calculated with FG MD, CG MD, CG LE (BOD), and CG GLE (BOD and IOMK) simulations. The plots for the other two mixtures are shown in Figs. S14 and S15. The CG MD simulation data show a rapid decay due to the absence of friction. This decay slows down significantly in the CG LE (BOD) simulation. The decay of the first peak height computed with the GLE (BOD) model is quite close to the FG MD counterpart. CG GLE (IOMK) simulations lead to a further improvement with excellent agreement at around t = 0.3 ps and beyond. A detailed discussion on the performance of all the models in capturing the decay of the first peak of gd(R, t) is provided in the supplementary material.

As a step closer to complex molecular systems, we also derived a coarse-grained model for a non-ideal ETH-TCM mixture at a number ratio of 75:25. It should be noted that this system deviates quite strongly from an ideal mixture,43 which is reflected in the inhomogeneous spatial distribution of the two solution components as shown in Fig. S4. The applicability of various systematic coarse-graining methods to derive CG conservative interactions for non-ideal mixtures is not straightforward.63,64 Here, we test our methodology to derive dynamical coarse-grained models for the 75:25 ETH-TCM system, for which we were able to converge the CG conservative IBI interactions (Fig. S4).

We derived the memory kernels for the GLE thermostats using the IOMK method. The difference between the target integrated effective single-particle memory kernel Gtarget(t) and the effective integrated memory in CG MD simulation a(t), i.e., Gtarget(t) − a(t) (Fig. S16), was fitted to a combination of six damped harmonic oscillators, and the resulting fit function was taken as the initial guess at the zeroth iteration, i.e., G̃0(t). As in the case of the CHX-TCM mixture, the two memory kernels were optimized alternately (i.e., one kernel was kept fixed at any iteration and the other was optimized). Figure S17 shows the VACFs and D(t) for the two components. Excellent agreement with their FG counterparts was obtained within 24/25 iterations. In this sequence of iterations, each memory kernel was updated twelve times. Figure 5 shows the decay of the first peak of gd(R, t) for the two components at the corresponding peak positions: R = 0.41 nm for ETH and R = 0.56 nm for TCM. For ETH, gd(R, t) calculated with the CG MD model decays significantly faster in comparison with TCM. This observation may be attributed to the inability of the conservative (IBI) interactions to capture the rate at which the hydrogen bonded structure around the polar ETH molecule breaks up. By contrast, the CG GLE (IOMK) model predicts a slower decay of the structural correlations between the ETH particles as compared to the FG model. The first peak of the ETH component has fully decayed after 1 ps with all models, which in the FG case would correspond to a complete alteration of the hydrogen bond network around the central molecule. For the TCM component, the decay of the first peak observed with the CG GLE (IOMK) model is in close agreement with the FG model on time scales larger than 0.5 ps.

FIG. 5.

Decay of the first peak of the distinct van Hove functions for (a) ETH and (b) TCM obtained from FG MD and CG GLE (IOMK) simulations of the 75:25 ETH-TCM mixture.

FIG. 5.

Decay of the first peak of the distinct van Hove functions for (a) ETH and (b) TCM obtained from FG MD and CG GLE (IOMK) simulations of the 75:25 ETH-TCM mixture.

Close modal

In a recent study on liquid water,40 it has been shown that the agreement between gd(R, t) computed with the CG GLE (IOMK) and FG models can be improved when the conservative CG interactions between water molecules are modeled with a multibody potential rather than with an IBI pair potential. We, therefore, expect that the discrepancy between the data obtained with the GLE IOMK model and the FG MD model in Fig. 5(a) is caused by the limited accuracy of the ETH-ETH IBI pair potential.

The satisfactory agreement achieved in representing the single-particle dynamics as well as gd(R, t) with an effective IBI pair potential and an isotropic GLE IOMK thermostat is encouraging. It suggests that, in addition to single-particle dynamics, this CG modeling approach provides a way to investigate the dynamic processes of complex solutions in which correlated particle motions are important.

In this work, we derived dynamical coarse-grained models for pure component molecular liquids and ideal and non-ideal binary mixtures. We found that FG dynamical observables such as the VACF, its running time integral D(t), and the long time diffusion coefficients D = limtD(t) can be quantitatively represented with these CG models by employing an isotropic GLE thermostat. Two bottom-up schemes were used to parameterize the GLE thermostat: the direct parametrization (BOD)30 method and the iterative parametrization (IOMK)40,41 method. Based on an analysis of the distinct van Hove function, we demonstrated that this approach also provides an accurate representation of the rate at which the structure of the liquid breaks up, in particular when the IOMK-based GLE is employed.

The accuracy of direct parameterization methods is generally limited by approximations as, in general, exact CG EoMs with non-linear conservative interactions are known to take more complex forms than Eq. (1).52,53,65,66 While direct parameterization using the BOD method is shown to yield quite accurate results, the IOMK method is easy to implement and converges rapidly. It can be used as a convenient and useful tool for dynamic coarse-graining purposes in future studies. This is especially relevant for chemically specific systems, where the conservative interactions are, in general, non-linear. The successful application to single-component liquids and ideal and non-ideal molecular mixtures reported herein is encouraging and suggests that the dynamics of structural relaxations in multi-component systems can be addressed with CG models augmented with isotropic GLE thermostats that couple to the dynamics at the single-particle level.

See the supplementary material for additional data supporting the findings of this work that are referenced but not included in the main text. We also include discussions on the contributions to effective single particle memory kernels and the decay of van Hove functions.

This work was funded by the Deutsche Forschungsgemeinschaft (DFG) via SFB TRR 146 (Grant No. 233630050).

The authors have no conflicts to disclose.

M.T. and V.K. contributed equally to this work.

Madhusmita Tripathy: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (supporting); Software (supporting); Writing – original draft (lead); Writing – review & editing (equal). Viktor Klippenstein: Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (lead); Software (lead); Writing – review & editing (equal). Nico F. A. van der Vegt: Conceptualization (equal); Funding acquisition (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material and openly available in TU Datalib at https://doi.org/10.48328/tudatalib-1177.67 

2.
C.
Peter
and
K.
Kremer
,
Faraday Discuss.
144
,
9
(
2010
).
3.
E.
Brini
,
E. A.
Algaer
,
P.
Ganguly
,
C.
Li
,
F.
Rodríguez-Ropero
, and
N. F. A.
van der Vegt
,
Soft Matter
9
,
2108
(
2013
).
4.
W. G.
Noid
,
J. Chem. Phys.
139
,
090901
(
2013
).
5.
J.
Jin
,
A. J.
Pak
,
A. E. P.
Durumeric
,
T. D.
Loose
, and
G. A.
Voth
,
J. Chem. Theory Comput.
18
,
5759
(
2022
).
6.
W. G.
Noid
,
J. Phys. Chem. B
127
,
4174
(
2023
).
7.
A. P.
Lyubartsev
and
A.
Laaksonen
,
Phys. Rev. E
52
,
3730
(
1995
).
8.
D.
Reith
,
M.
Pütz
, and
F.
Müller-Plathe
,
J. Comput. Chem.
24
,
1624
(
2003
).
9.
S.
Izvekov
and
G. A.
Voth
,
J. Phys. Chem. B
109
,
2469
(
2005
).
10.
E.
Brini
,
V.
Marcon
, and
N. F. A.
van der Vegt
,
Phys. Chem. Chem. Phys.
13
,
10468
(
2011
).
11.
M. S.
Shell
,
J. Chem. Phys.
137
,
084503
(
2012
).
12.
J. F.
Rudzinski
and
W. G.
Noid
,
J. Phys. Chem. B
118
,
8295
(
2014
).
13.
M. P.
Bernhardt
,
M.
Hanke
, and
N. F. A.
van der Vegt
,
J. Chem. Theory Comput.
19
,
580
(
2023
).
14.
W.
Tschöp
,
K.
Kremer
,
J.
Batoulis
,
T.
Bürger
, and
O.
Hahn
,
Acta Polym.
49
,
61
(
1998
).
15.
D.
Fritz
,
V. A.
Harmandaris
,
K.
Kremer
, and
N. F. A.
van der Vegt
,
Macromolecules
42
,
7579
(
2009
).
16.
H. A.
Karimi-Varzaneh
,
N. F. A.
van der Vegt
,
F.
Müller-Plathe
, and
P.
Carbone
,
Chem. Phys. Chem.
13
,
3428
(
2012
).
18.
H. I.
Ingòlfsson
,
C. A.
Lopez
,
J. J.
Uusitalo
,
D. H.
de Jong
,
S. M.
Gopal
,
X.
Periole
, and
S. J.
Marrink
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
225
(
2014
).
20.
V.
Klippenstein
,
M.
Tripathy
,
G.
Jung
,
F.
Schmid
, and
N. F. A.
van der Vegt
,
J. Phys. Chem. B
125
,
4931
(
2021
).
21.
R.
Zwanzig
,
J. Chem. Phys.
33
,
1338
(
1960
).
22.
H.
Mori
,
Prog. Theor. Phys.
33
,
423
(
1965
).
23.
R.
Zwanzig
,
Nonequilibrium Statistical Mechanics
(
Oxford University Press
,
New York
,
2001
).
24.
C.
Hijón
,
P.
Español
,
E.
Vanden-Eijnden
, and
R.
Delgado-Buscalioni
,
Faraday Discuss.
144
,
301
(
2010
).
25.
R. D.
Groot
and
P. B.
Warren
,
J. Chem. Phys.
107
,
4423
(
1997
).
26.
H.
Lei
,
B.
Caswell
, and
G. E.
Karniadakis
,
Phys. Rev. E
81
,
026704
(
2010
).
27.
G.
Deichmann
,
V.
Marcon
, and
N. F. A.
van der Vegt
,
J. Chem. Phys.
141
,
224109
(
2014
).
28.
S.
Trèment
,
B.
Schnell
,
L.
Petitjean
,
M.
Couty
, and
B.
Rousseau
,
J. Chem. Phys.
140
,
134113
(
2014
).
29.
G.
Deichmann
and
N. F. A.
van der Vegt
,
J. Chem. Phys.
149
,
244114
(
2018
).
30.
V.
Klippenstein
and
N. F. A.
van der Vegt
,
J. Chem. Phys.
154
,
191102
(
2021
).
31.
J.
Kappler
,
J. O.
Daldrop
,
F. N.
Brünig
,
M. D.
Boehle
, and
R. R.
Netz
,
J. Chem. Phys.
148
,
014903
(
2018
).
32.
Z.
Li
,
X.
Bian
,
X.
Li
, and
G. E.
Karniadakis
,
J. Chem. Phys.
143
,
243128
(
2015
).
33.
Z.
Li
,
H. S.
Lee
,
E.
Darve
, and
G. E.
Karniadakis
,
J. Chem. Phys.
146
,
014104
(
2017
).
34.
Y.
Yoshimoto
,
Z.
Li
,
I.
Kinefuchi
, and
G. E.
Karniadakis
,
J. Chem. Phys.
147
,
244110
(
2017
).
35.
S.
Wang
,
Z.
Li
, and
W.
Pan
,
Soft Matter
15
,
7567
(
2019
).
36.
G.
Jung
,
M.
Hanke
, and
F.
Schmid
,
J. Chem. Theory Comput.
13
,
2481
(
2017
).
37.
G.
Jung
,
M.
Hanke
, and
F.
Schmid
,
Soft Matter
14
,
9368
(
2018
).
38.
A.
Davtyan
,
J. F.
Dama
,
G. A.
Voth
, and
H. C.
Andersen
,
J. Chem. Phys.
142
,
154104
(
2015
).
39.
A.
Davtyan
,
G. A.
Voth
, and
H. C.
Andersen
,
J. Chem. Phys.
145
,
224107
(
2016
).
40.
V.
Klippenstein
and
N. F. A.
van der Vegt
,
J. Chem. Theory Comput.
19
,
1099
(
2023
).
41.
V.
Klippenstein
and
N. F. A.
van der Vegt
,
J. Chem. Phys.
157
,
044103
(
2022
).
42.
R.
Hayes
,
G. G.
Warr
, and
R.
Atkin
,
Chem. Rev.
115
,
6357
(
2015
).
43.
G.
Guevara-Carrion
,
T.
Janzen
,
Y. M.
Muñoz-Muñoz
, and
J.
Vrabec
,
J. Chem. Phys.
144
,
000001
(
2016
).
44.
W. L.
Jorgensen
,
D. S.
Maxwell
, and
J.
Tirado-Rives
,
J. Am. Chem. Soc.
118
,
11225
(
1996
).
45.
M. J.
Abraham
,
T.
Murtola
,
R.
Schulz
,
S.
Páll
,
J. C.
Smith
,
B.
Hess
, and
E.
Lindahl
,
SoftwareX
1-2
,
19
(
2015
).
46.
B.
Hess
,
J. Chem. Theory Comput.
4
,
116
(
2008
).
47.
T.
Darden
,
D.
York
, and
L.
Pedersen
,
J. Chem. Phys.
98
,
10089
(
1993
).
48.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
,
J. Chem. Phys.
126
,
014101
(
2007
).
49.
H. J. C.
Berendsen
,
J. P. M.
Postma
,
W. F.
van Gunsteren
,
A.
DiNola
, and
J. R.
Haak
,
J. Chem. Phys.
81
,
3684
(
1984
).
51.
M.
Parrinello
and
A.
Rahman
,
J. Appl. Phys.
52
,
7182
(
1981
).
52.
H.
Vroylandt
and
P.
Monmarché
,
J. Chem. Phys.
156
,
244105
(
2022
).
53.
H.
Vroylandt
,
Europhys. Lett.
140
,
62003
(
2022
).
54.
V.
Rühle
,
C.
Junghans
,
A.
Lukyanov
,
K.
Kremer
, and
D.
Andrienko
,
J. Chem. Theory Comput.
5
,
3211
(
2009
).
55.
V.
Rühle
and
C.
Junghans
,
Macromol. Theory Simul.
20
,
472
(
2011
).
56.
A.
Carof
,
R.
Vuilleumier
, and
B.
Rotenberg
,
J. Chem. Phys.
140
,
124103
(
2014
).
57.
A. V.
Straube
,
B. G.
Kowalik
,
R. R.
Netz
, and
F.
Höfling
,
Commun. Phys.
3
,
126
(
2020
).
58.
M.
Ceriotti
,
G.
Bussi
, and
M.
Parrinello
,
Phys. Rev. Lett.
102
,
020601
(
2009
).
59.
M.
Ceriotti
,
G.
Bussi
, and
M.
Parrinello
,
J. Chem. Theory Comput.
6
,
1170
(
2010
).
60.
61.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in ’t Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
,
R.
Shan
,
M. J.
Stevens
,
J.
Tranchida
,
C.
Trott
, and
S. J.
Plimpton
,
Comput. Phys. Commun.
271
,
108171
(
2022
).
62.
J. P.
Hansen
and
I. R.
McDonald
,
Theory of Simple Liquids
(
Academic Press
,
Oxford
,
2013
).
63.
A.
Villa
,
C.
Peter
, and
N. F. A.
van der Vegt
,
J. Chem. Theory Comput.
6
,
2434
(
2010
).
64.
T.
Sanyal
and
M. S.
Shell
,
J. Phys. Chem. B
122
,
5678
(
2018
).
65.
F.
Glatzel
and
T.
Schilling
,
Europhys. Lett.
136
,
36001
(
2022
).
66.
C.
Ayaz
,
L.
Scalfi
,
B. A.
Dalton
, and
R. R.
Netz
,
Phys. Rev. E
105
,
054138
(
2022
).
67.
M.
Tripathy
,
V.
Klippenstein
, and
N. F. A.
van der Vegt
(
2023
), “
Data of the publication ‘dynamical coarse-grained models of molecular liquids and their ideal and non-ideal mixtures’
,” TU Datalib.

Supplementary Material