Dynamic force matching: Construction of dynamic coarse-grained models with realistic short time dynamics and accurate long time dynamics

As computational power increases, it becomes possible to perform all-atom molecular dynamics simulations of increasingly larger systems of larger molecules for longer times. These simulations describe the dynamics of every atom in the system and the structure of the system at an atomic level of detail on a sub-picosecond time scale out to time scales currently as long as milliseconds. Even so, at the present time, there are problems of interest that are too demanding to be addressed in this way by available computational resources.

A strategy frequently used to deal with this is to develop coarse-grained (CG) models of the molecular systems of interest. These models have fewer degrees of freedom than fine grained (FG) all-atom models and generally simpler potential energy functions than FG models, thus facilitating the computations. Much effort has been directed towards choosing appropriate degrees of freedom for CG models as well as constructing CG potential energy surfaces,^1–21 and there also has been some effort directed toward developing dynamic CG models.^22–38 

The Multi-Scale Coarse-Graining (MS-CG) method developed by Voth and coworkers can be used to construct CG models and calculate the potential of mean force of the CG degrees of freedom.^39–44 In the MS-CG method, the CG degrees of freedom for an FG system are defined by (1) choosing which molecules are to be included in the CG model,⁴⁵ (2) partitioning the atoms on each of these molecules into sets, and (3) defining a CG “site” for each set. The CG site for a set of atoms is a point located at the position of the center of mass of the atoms in the set. The CG model is constructed by regarding each site as a CG “particle” whose mass is the total mass of the atoms in the corresponding set and whose position is the position of the corresponding site. For example, if a methanol molecule is partitioned into two sets of atoms, one containing the CH₃ group and one containing the OH group, and if a CG model is defined for this partition, then in the CG model each methanol molecule consists of one methyl “CG particle” and one hydroxyl “CG particle.”

The MS-CG method provides a variational procedure, sometimes called force matching, to calculate a many-dimensional potential of mean force of the CG model using FG simulation data. That potential of mean force is the potential that governs the canonical ensemble equilibrium distribution in phase space for the CG model. The MS-CG method can in principle provide the exact CG potential of mean force as a function of the CG coordinates, but in practice the variational result will be approximate because the set of basis functions used to represent the CG potential will likely not be complete enough.

The CG particles of an MS-CG model, moving according to ordinary Hamiltonian mechanics under the influence of the potential of mean force, will move much more rapidly through phase space than the corresponding CG sites in the FG system. This is because the other degrees of freedom of the FG system, which are absent from the MS-CG system, exert fluctuating forces on the CG sites in the FG system. These fluctuating forces act as a source of friction for the motion of the CG sites in the FG system.

In a recent paper,³² we have presented a method, which we call “dynamic force matching,” for constructing a dynamic MS-CG model that has some (but not all) of the dynamical properties of the FG system. To construct a dynamic CG model, we add degrees of freedom to the MS-CG model that will exert fluctuating forces on the CG particles but that will not change the potential of mean force of the CG particles. These additional degrees of freedom are represented as particles that we refer to as “fictitious particles.” The interactions are adjusted so that the forces that the fictitious particles generate have dynamical effects and statistical properties that are consistent with the forces experienced by the CG sites in the FG system. Their interactions with the CG particles are very simple and can be computed efficiently.

The fictitious particle models that we use are constructed so that they are special cases of a class of models discussed by Zwanzig.^46,47 These models have a number of important properties.

1. The canonical equilibrium distribution function of the CG variables in the dynamic CG model is the same as the equilibrium distribution function determined by the CG potential. Thus the fictitious particles affect the dynamics of the CG particles without changing their equilibrium distribution.

2. The fictitious particles undergo Langevin dynamics subject to the forces exerted on them by the CG particles, but the CG particles undergo Hamiltonian and Newtonian dynamics subject to the forces exerted on them by the fictitious particles.

3. Despite the fact that the Langevin forces on the fictitious particles are unbounded and discontinuous functions of time, the forces on the CG particles are finite and continuous functions of time and thus the velocities of the CG particles have finite continuous time derivatives. This leads to time correlation functions of various CG variables that are physically reasonable.

4. The time evolution of the CG degrees of freedom is not a Markovian process. Zwanzig’s work showed that the CG degrees of freedom by themselves satisfy a generalized Langevin equation. (The situation is very different in models in which the CG particles are directly connected to a Langevin heat bath. In that case the time evolution of the CG degrees of freedom is a Markovian process and satisfies a Langevin equation.²³)

5. The time evolution of the combined CG and fictitious particle degrees of freedom is Markovian and is easily calculated using standard programs since it involves Hamiltonian mechanics with only a subset of the degrees of freedom (the fictitious particles) subject to Langevin forces.

In the present paper, we present several developments of the dynamic force matching method.

1. We describe and use a new type of fictitious particle, called type-3, that exerts fluctuating forces on the center of mass of two CG particles. We expect that this type of particle, and simple generalizations of it, will be very important for all future use of these kinds of models.

2. We present improved algorithms for some of the procedures used in implementing the method.

3. We show how to construct fictitious particle models that are coarse-grained in time, i.e., we show that it is possible to systematically remove from a dynamic CG model the high frequency motions that complicate the dynamical calculations without having much net effect on the dynamics of the CG degrees of freedom at intermediate and long times.

4. We construct and parametrize several dynamic CG models for methanol liquid that have type-3 particles, including one whose two basic transport and relaxation properties, the self diffusion coefficient and the rotational relaxation time, have the same values as those of the FG model for methanol liquid.

In this section, we review and generalize some of the ideas in the previous paper³² regarding dynamic CG models. For further details, see that paper.

In this section, we discuss the use of fictitious particles to convert an MS-CG model into one that has a realistic dynamical behavior. We also discuss the properties of three types of fictitious particles.

An MS-CG model is constructed by starting with an FG system (usually an all-atom model), defining the CG sites that are chosen for each molecule in the system, and performing a variational calculation of the potential of mean force of those sites. The MS-CG model is defined as a system that has one CG point particle for each CG site in the FG system, and the potential energy of interaction of these CG particles is the potential of mean force of the sites in the FG system.

Let R^N denote the positions of the CG particles. Let P^N denote their momenta. Then the Hamiltonian of the MS-CG system is

Each CG particle has a label I, taking on values $1,2,…,N$⁠, where N is the number of CG particles in the system. M_I is the mass of CG particle I. V(R^N) is the MS-CG potential of mean force.

A dynamic CG model is created by the following:

1. adding fictitious particles to the MS-CG model;

2. choosing a Hamiltonian for the combined system of CG particles and fictitious particles;

3. defining the dynamics of the combined system to be obtained from the Hamiltonian equations of motion for the CG and fictitious particles, supplemented by Langevin forces acting on each of the fictitious particles;

4. choosing values for the parameters that describe the essential properties of the fictitious particles.

We shall use the terms “dynamic CG model” and “fictitious particle model” interchangeably in the following to refer to a model that contains both CG particles and fictitious particles.

Each fictitious particle has a label α, taking on values $1,2,…,M$⁠, where M is the number of fictitious particles in the system. Each fictitious particle is of a specific type. The type of particle α is denoted as $sα$⁠, where $sα$ takes on values that are small integers. Three types have been defined, and the corresponding values of s are 1, 2, and 3. (In the future, additional types of fictitious particles may be defined and used.) A particle whose type equals s will be called a “type-s particle.” Type-1 and type-2 particles were defined and used in the previous paper. The type-3 particle is discussed for the first time in the present paper.

Each fictitious particle interacts with only a small number of CG particles.

1. A type-1 particle interacts with only one CG particle. If α is the label of a type-1 fictitious particle, we define $pα$ to be the label of the CG particle that fictitious particle α interacts with.

2. A type-2 or a type-3 particle interacts with only two CG particles on the same CG molecule. If α is the label of a type-2 or type-3 fictitious particle, we define $pα$ and $qα$ to be the labels of the two CG particles that fictitious particle α interacts with, where $pα<qα$⁠.

For this first step in constructing a fictitious particle model, the fundamental choices that must be made are which types of fictitious particles should be used and what CG particles should they be coupled to.

The total Hamiltonian of the combined system is of the following form:

S^M and Q^M denote the positions and momenta of the M fictitious particles and

where

Each fictitious particle α is a harmonic oscillator with mass $μα$⁠, angular frequency $ωα(>0)$⁠, and Langevin damping constant $γα$ (see Appendix B 3 for the definition of the damping constant). The location of the minimum of its potential energy of interaction is determined by the positions of the CG particle or particles it interacts with. The value of this potential energy of interaction is zero at the minimum.

1. For a type-1 particle, the minimum is located at the position of the one CG particle it interacts with (see Eq. (3)).

2. For a type-2 particle, the minimum is located at the vector distance between the two CG particles that it interacts with (see Eq. (4)).

3. For a type-3 particle, the minimum is located at the center of mass of the two CG particles that it interacts with (see Eq. (5)).

In a computer simulation using classical mechanics and periodic boundary conditions, the fictitious particles will behave in the following ways.

1. A type-1 fictitious particle will stay near the CG particle it interacts with and hence near the molecule that contains that CG particle.

2. A type-3 particle will stay near the center of mass of the two CG particles it interacts with. Since these two CG particles are on the same molecule, the fictitious particle will stay near the molecule that contains the two CG particles.

3. A type-2 particle will have an unusual behavior. The two CG particles with which it interacts are on the same molecule and hence are constrained by the MS-CG potential to stay near one another; the distance between them fluctuates in a range that is consistent with the intramolecular interactions of the original FG model. $|𝐑pα−𝐑qα|$ is in effect bounded by a value that is independent of the size of the system. $𝐒α$⁠, the set of coordinates of fictitious particle α, is similarly constrained by its harmonic potential to stay near the value of $𝐑pα−𝐑qα$⁠. Thus, the vector coordinate of a type-2 particle will, in effect, be bounded in magnitude, and the value of the bound is independent of the size of the system.

Particles of type-1 and type-3 can be visualized as part of the CG molecule to which they are attached. They move along with that molecule and satisfy periodic boundary conditions of the same type that the CG particles satisfy.

Particles of type-2, however, cannot be visualized in this way. They are in a three dimensional space of their own. Each type-2 particle tends to stay near the origin of that space, regardless of the motions of the two CG particles with which it interacts. As the molecule that contains these two CG particles moves from one cell to another under periodic boundary conditions, the coordinates of the type-2 particle are not affected.

With this choice of Hamiltonian, the dynamic CG model is in the class of models defined by Zwanzig.^46,47 The features of the Hamiltonian that are sufficient to guarantee this are that the fictitious particles are harmonic oscillators with constant frequency, the fictitious particles do not interact with each other, and the location of the minimum in the potential energy of each fictitious particle is a linear function of the coordinates of the CG particles. Moreover, with this choice of Hamiltonian, it is straightforward to show that the fictitious particles have no effect on the equilibrium probability distribution of the positions of the CG particles.

The dynamics of the fictitious particle model are those of classical mechanics of the CG and fictitious particles generated by the Hamiltonian above, with each fictitious particle being also subject to Langevin forces of the usual type. (It is straightforward to show that even with Langevin forces that act on the fictitious particles only, the model is still in the class of models defined by Zwanzig and still has the property that the fictitious particles have no effect on the equilibrium probability distribution of the positions of the CG particles.)

There are several important features to this dynamics.

1. Since there are no Langevin forces acting on the CG particles, the velocity of each CG particle has a continuous first derivative with regard to time, and hence the dynamics and the correlation functions of the CG sites do not have some of the unphysical features that would appear if they were directly coupled to a Langevin heat bath.

2. A type-1 particle exerts fluctuating forces only on the one CG particle it interacts with.

3. A type-2 particle exerts fluctuating forces only on the two CG particles it interacts with. The two forces are equal in magnitude and opposite in directions. They can exert torques on the two particles and can act to increase or decrease the distance between those particles, but the net momentum transferred to the two particles is exactly zero at all times.

4. A type-3 particle exerts fluctuating forces only on the two CG particles it interacts with. These forces are in the same direction for both particles and the ratio of their magnitudes is equal to the ratio of the masses of the two CG particles. As a result, the two particles receive exactly the same acceleration from the fictitious particle. This acts to accelerate the center of mass of the two CG particles, but it exerts no torques on them and it does not change the distance between the two particles.

5. The Hamiltonian in Eq. (2) is unchanged if the system is translated in space such that as follows:

   • each CG particle coordinate and each type-1 and type-3 fictitious particle are incremented by the same vector displacement;

   • the coordinates of each type-2 fictitious particle are unchanged;

   • the momentum of each particle is unchanged.

6. From this it can be shown that, in Hamiltonian dynamics, there is a conserved momentum that is equal to the sum of the momenta of the CG particles and the type-1 and type-3 particles. When Langevin forces act on the fictitious particles, this momentum is not conserved but instead fluctuates about its initial value, which in practice should be set to zero in a simulation.

7. Similarly, since there is no explicit time dependence on the Hamiltonian, the Hamiltonian is a conserved energy in Hamiltonian dynamics, but Langevin forces on the fictitious particles will make this energy fluctuate about its initial value and eventually achieve a canonical equilibrium distribution.

The last step in constructing a dynamic CG model is that of choosing the numerical values of the parameters for fictitious particles. This involves the following choices:

1. for each CG site, the number of type-1 fictitious particles that interact with it;

2. for each pair of CG sites on the same molecule, the numbers of type-2 and of type-3 fictitious particles that interact with it;

3. for each fictitious particle α, the numerical values of the oscillator frequency $ωα$⁠, the mass $μα$⁠, and the Langevin friction constant $γα$ (see Appendix B 3).

The parametrization is based on making some of the dynamical properties of the fictitious particle model be the same as those of the FG model. The previous paper discusses two ways of parametrizing such models. They are called CG[1] and CG[2]. These methods are discussed further in Sec. II C.

The parametrization of dynamic CG models makes use of information about the “constrained dynamics” of both the FG and the dynamic CG systems. Here we review some of the properties and characteristics of constrained dynamics. (For a more precise mathematical discussion of the dynamics and statistical mechanics associated with constrained dynamics, see the previous paper.)

Consider an FG system for which we have chosen the various CG sites that will be used to construct a coarse-grained model. Each such site is the center of mass of a certain set of atoms on the same molecule.

In constrained FG dynamics, the motion is determined by the usual FG Hamiltonian subject to the constraint that each CG site on each molecule does not move. Thus, if a CG site is defined in terms of one atom, that atom is not allowed to move. If a CG site is defined as the center of mass of two or more atoms on the same molecule, the individual atoms associated with that site can move, but only in such a way that their center of mass remains fixed.

Information about constrained FG dynamics can be obtained only by performing constrained simulations of the FG system. The purpose of performing such constrained simulations is to evaluate the correlation functions of the fluctuating forces that act on the CG sites when the sites are constrained not to move. Such correlation functions are the result of the interaction of the CG degrees of freedom with those degrees of freedom that are eliminated in the coarse-graining process.

Constrained dynamics of a dynamic CG model is very similar. The particles in the dynamic CG model move according to the Hamiltonian in Eq. (2) and the fictitious particles are subject to Langevin forces; however, the CG sites are constrained not to move. There is no need to simulate constrained dynamics for a dynamic CG model, because the information of interest, namely, the correlation functions of the fluctuating forces acting on the CG particles, can be expressed analytically in terms of the parameters of the particles. (This is discussed in the previous paper.)

The CG[1] method of parametrizing a fictitious particle model uses the time dependence of the correlation functions of the fluctuating forces on CG sites that are obtained from constrained dynamics simulation of the FG model. The method finds the parameters of the fictitious particles that make the fictitious particle model have fluctuation force correlation functions for constrained dynamics whose time dependence is the same as that of the FG model.

For times short enough that the CG sites have not moved very far, it is reasonable to expect that the fluctuating forces in a constrained FG simulation would be similar to the fluctuating forces in normal (unconstrained) dynamics. Moreover, in the case that the CG sites are massive and move very slowly, the similarity between constrained and unconstrained dynamics may extend to rather long times. This is part of the motivation for use of the CG[1] method.

CG[1] is also a very practical method of obtaining all the parameters of the model in a straightforward way. It requires only a constrained dynamical simulation of the FG system to obtain the correlation functions and an optimization method for choosing the parameters of the fictitious particles. It requires no simulations of the CG system.⁴⁸ 

The GC[2] method starts with a CG[1] parametrization and modifies a small number of the CG[1] parameters so that certain properties of the fictitious particle model obtained from normal (unconstrained) dynamical simulations are the same as those obtained from normal (unconstrained) simulations of the FG model. This makes it possible to insure that the dynamic model has values of some important dynamic properties, such as the self-diffusion coefficient, that are the same as those of the FG model. This method requires that some simulations of the fictitious particle model be performed.

The two-site model of liquid methanol that we will discuss is described in the Introduction. Each molecule consists of two CG particles: a methyl CG particle and a hydroxyl CG particle. The fine grained (FG) all-atom system on which the model is based is a system of 1000 methanol molecules in a cubic box of length 40.9 Å, whose molecules interact with each other with the Optimized Potentials for Liquid Simulations all-atom (OPLS-AA)⁴⁹ force field, and which is at equilibrium at 300 K. The density corresponds approximately to the density of real methanol at this temperature and 1 atmosphere pressure.

Unless otherwise indicated, all quantities related to the simulation of FG and CG systems are expressed using the Angstrom, the picosecond, and the kcal/mol as the units of length, time, and energy. Thus, for example, force has units of (kcal/mol) (Å)⁻¹, and frequency has units of ps⁻¹.

The construction of an MS-CG force field for the CG system is discussed in Sec. VI of the previous paper.³² See the supplementary material for that paper for the MS-CG potentials and radial distribution functions of this model. The MS-CG potential consists of a sum of intermolecular potentials for interactions of CH₃ with CH₃, CH₃ with OH, and OH with OH, and an intramolecular CH₃–OH potential.

A dynamic CG model for liquid methanol was constructed by adding fictitious particles to each CG molecule. In the present work, we have used type-2 and type-3 fictitious particles, which are discussed above in Sec. II A. We have parametrized the model using the CG[1] and CG[2] parametrization methods discussed in Sec. II C.

In this section, we discuss this model and its two parametrizations. We also discuss less elaborate models with fewer fictitious particles that are obtained by coarse-graining in time, which is achieved by removing the higher frequency fictitious particles.

CG[1] parametrization of a dynamic CG model requires choosing certain correlation functions of the CG model and requiring them to be equal to the corresponding functions in the FG model. The principles of the CG[1] parametrization are summarized in Sec. II C and discussed in detail in the previous paper.³² The generalization of these principles to include type-3 particles is summarized in  Appendix A of the present paper. The following is a brief discussion of the way we have applied these principles in the present work. For more details see  Appendix A.

Consider the FG model for methanol. For each molecule, we define a CG site labeled K, located at the center of mass of the atoms in the methyl group, and a CG site labeled L, located at the center of mass of the atoms in the hydroxyl group.

In constrained dynamics of the FG model, there are fluctuating forces exerted on the two sites. Let $δFK(t)$ and $δFL(t)$ refer to the fluctuating forces on the two sites. Consider the following correlation functions:

(Here the angular brackets $⟨…⟩cFG$ refer to a canonical ensemble average for the constrained FG system. See the previous paper for a discussion of the equilibrium ensemble for the constrained system.) $Ξ+FG(K,L,t)$ is the autocorrelation function of the sum of the fluctuating forces acting on the two distinct sites of the same molecule when the dynamics is constrained. $Ξ−FG(K,L,t)$ is a similar autocorrelation function for the difference of the fluctuating forces on the two sites. These functions can be obtained straightforwardly from constrained simulations of the FG system.

Analogous correlation functions can also be defined for the dynamic coarse-grained system that has both CG particles and fictitious particles. Let $δFK(t)$ and $δFL(t)$ refer to the fluctuating forces of two CG particles on the same CG molecule in a constrained system,

The CG[1] parametrization that we have used is based on the requirement that these two correlation functions for the dynamic CG model, under constrained dynamics, should have the same time dependence as the analogous functions of the FG model under constrained dynamics. That is, we require that

These two requirements, on Ξ ₊ and Ξ ₋, are sufficient to determine the properties of the two-site fictitious particle model for methanol when only type-2 and type-3 fictitious particles are allowed in the model.

Thus, in the CG[1] parametrization, the number of fictitious particles and the values of their parameters are chosen so that the functions $Ξ±CG(K,L,t)$⁠, calculated from analytic formulas, fit the data for the $Ξ±FG(K,L,t)$ obtained from constrained FG simulations. In practice, we work in the Fourier domain and we require that the Fourier transforms of the $Ξ±CG(K,L,t)$ functions fit the Fourier transforms of the $Ξ±FG(K,L,t)$ data for all frequencies below a cutoff frequency.

The analytic formulas for $Ξ±CG(K,L,t)$ and the procedures and algorithms for applying this criterion to perform the CG[1] parametrization are discussed in Appendices A 3 and  B.

Based on the method described above, we performed a CG[1] parametrization that resulted in a model that contains 11 fictitious particles, which we shall refer to as “CG[1] 11fp.” Fig. 1 shows the Fourier transforms of $Ξ±FG$ functions obtained from constrained FG simulations and the best fits to these functions obtained by using the functional forms of the $Ξ±CG$ functions and varying the number of fictitious particles of each type and their parameters. (In the fitting process, no attempt was made to fit the functions for frequencies greater than 150 for the graph on the left and 100 for the graph on the right.) Fig. 2 shows the same functions in the time domain. It is evident that the fits obtained this way filter out the very fast oscillations.

We performed molecular dynamics simulation of the CG[1] 11fp model, and the results will be discussed below. In preparation for that discussion, we note a number of features of the parameters obtained for the fictitious particles. (See Table I for values of the parameters and their definitions. Also see Appendix B 3 for a more detailed discussion of these parameters.)

1. The zero frequency Fourier transform of an autocorrelation function of the fluctuation force exerted on an FG site in a constrained system can be regarded as a friction coefficient for the FG sites in the unconstrained system if the CG sites are massive enough so that they move very short distances on the time scale of the decay of the autocorrelation function. The ff0 of a fictitious particle is proportional to its contribution to fitting the zero frequency Fourier transform of the $Ξ±FG$ functions. For this model, the fictitious particles with the smallest frequencies have the largest values of ff0 and hence largest frictional effect.

2. The t = 0 value of the autocorrelation function of the force exerted on an FG site in a constrained system is the same as that for an unconstrained system. The $vt$0 of a fictitious particle is its contribution to fitting the zero time values of the $Ξ±FG$ functions. For the CG[1] 11fp model, the fictitious particles with the largest frequencies exert the largest forces. Despite this, they have the smallest frictional effect on the CG particles, according to the ff0 criterion discussed above.

3. The parameter aof of a fictitious particle is an approximation for the oscillation frequency of the force exerted by that fictitious particle. The contribution of a fictitious particle to fitting the Fourier transform of the functions has a peak approximately at aof. In general, if there are peaks in the Fourier transforms of the $Ξ±FG$ functions, the fictitious particle model obtained by CG[1] parametrization will have fictitious particles whose aof is approximately equal to the location of the observed peaks.

Coarse-grained models constructed using the MS-CG method are designed so that they do not describe the properties of the molecules on very short distance scales. (In the case of the two-site model of methanol, the precise locations of the H atoms relative to the C and O atoms are not described by the CG model.) This elimination of short distance information from the model is deliberate and is done in the hope of obtaining a more economical description of the system on longer distance scales. In this sense, we can say that the MS-CG method constructs models that are “coarse-grained in space.”

The dynamic CG model we have constructed contains fictitious particles that exert forces on the CG particles on a range of time scales, as can be seen in Figs. 1 and 2. This raisesthe question of whether we can systematically remove the high frequency behavior and thereby generate models that are “coarse-grained in time.”

In fact this can be done by systematically removing the high frequency fictitious particles from the CG[1] model and by eliminating intramolecular vibration in the CG molecule by constraining the CG bond length to be constant.

Table I shows that the CG[1] 11fp model has fictitious particles with apparent oscillation frequencies aof as large as about 120 ps⁻¹.

Eliminating all the fictitious particles with aof greater than 20 ps⁻¹ gives a model with four fictitious particles, which we call CG[1] 4fp. All four of these fictitious particles are overdamped. Note that removing the higher frequency particles eliminated only the particles with the smallest values of ff0 (see the discussion in Sec. III B 1).

In the CG[1] 4fp model, there are two type-2 particles and two type-3 particles. For both types, one of the two particles had a much larger value of ff0. In both cases that particle is also the one with the smaller frequency ω. This suggests that for each of the types only the particle with the larger value of ff0 is important.

Retaining only this one particle for each type of fictitious particle, we obtain what we call the CG[1] 2fp model.

The CG intramolecular potential is very stiff and the methyl-hydroxyl bond length fluctuates (at high frequency) in a narrow range of values. It is straightforward to perform dynamical calculations in which the bond length is constrained to be fixed. We have done this by using the molecular dynamics algorithm RATTLE, which is a generalization of the velocity Verlet algorithm used in this work so far.

We have performed simulations of such a rigid model based on each of the three CG[1] models discussed so far. For all the properties we calculated, the results for a vibrating model and the corresponding rigid model were indistinguishable from each other on the scale of the figures presented below (except for a small oscillation in the force correlation functions for the vibrating molecule). From this we conclude that for the two-site model of methanol, constraining the bond from vibrating has little effect on the translational and rotational motion of the molecule, and a rigid model is as accurate as a vibrating model.

The MS-CG force matching method as presented by Noid et al.⁴² shows how to start with a vibrating FG model and derive a MS-CG vibrating model. The method presented in this paper provides a way of starting with such an MS-CG vibrating model and obtain a vibrating dynamic CG model by adding fictitious particles. Once that fictitious particle model is constructed, it is straightforward to convert the fictitious particle model into one in which the CG molecules are rigid, and such a conversion may be useful. The only question is whether the final conversion has a significant effect on the important properties of the model. For the two-site methanol models discussed here, the conversion has little effect on the translational and rotational motion of the molecules. This indicates that the effect of the fictitious particles on the translational and rotational motions is not significantly changed when this two-site model of methanol is made rigid. In particular, the forces they exert on the center of mass of the molecule and the torques they exert are essentially unchanged.

The CG[1] parametrizations used information about the FG system obtained only from constrained simulations.

As discussed in Sec. II C and the previous paper,³² we also performed CG[2] parametrizations. A CG[2] parametrization starts from a CG[1] parametrization and uses information obtained from dynamic simulations of the FG system to make modifications of a small number of parameters in the CG[1] model. This information is long time information and transport coefficient information that we would want the most useful CG model to be in agreement with.

For methanol, there are two pieces of information of this type: the self-diffusion coefficient and the rotational relaxation time. The self-diffusion coefficient is obtained from the slope at long times of the mean squared displacement of the center of mass of a molecule as a function of time. The rotational relaxation time is the negative of the reciprocal of the slope at long times of a semilog plot of the orientational correlation function vs. time (see Eq. (17) of the previous paper³²).

For CG[2] parametrization, we start with a CG[1] model and vary only the Langevin damping constant γ of the type-2 and type-3 fictitious particles whose γ is the smallest for its type (namely, the first particle for each type in Table I). For various choices of these two parameters, we perform molecular dynamics simulations to calculate the self-diffusion coefficient and the rotational relaxation time. With just a small number of such simulations, it was possible to converge on the γ values that gave results for these two properties that agreed with the results of the FG simulation.

We did this parametrization starting with the CG[1] 4fp model and the CG[1] 2fp model. The parameters for the latter model, which plays an important role in the following discussion, are in Table II.

We have parametrized and simulated a variety of dynamic CG models. Some were vibrating and some were rigid. Some had 11 fictitious particles, some had 4, and some had 2, with fewer particles implying a larger extent of coarse-graining in time. Some were parametrized using only short time information obtained from FG simulations (the various CG[1] parametrizations), and others were parametrized using, in addition, the self-diffusion coefficient and the orientational relaxation time obtained from FG simulations (the CG[2] parametrizations).

The two most useful models to discuss in detail are the CG[1] 11fp model and the CG[2] 2fp model. The former uses a great deal of short time information obtained from the FG simulation of methanol, including some data representing high frequency motion. The latter is highly coarse-grained in time and was parametrized to fit the FG self-diffusion coefficient and the orientation relaxation time, which contain information about the long time behavior of the FG model. A model of this type is more likely to be useful in performing very long simulations of the CG model because it accurately describes the time scales for translational and rotational motion of the molecules. The other models give results that are similar to, and intermediate between, those of the two models just mentioned.

Figs. 3–6 and Table III present dynamic properties for the CG[1] 11fp model and the CG[2] 2fp model and compare them with the results for the FG (all atom) model and the MS-CG model, which is the CG model with no fictitious particles. For each CG model, the data for the rigid version are presented, whereas the data for the FG model are for the original vibrating model. Table III also contains some information about the other dynamic CG models investigated.

The four dynamic properties are as follows: the autocorrelation function of the total force on a molecule (Fig. 3), the autocorrelation function of the velocity of the center of mass of a molecule (Fig. 4), the self-diffusion coefficient (Fig. 5), and the rotational relaxation time (Fig. 6). The first two properties are correlation functions that decay to zero in times of about 0.4 and 4 ps, respectively. The last two are a transport coefficient and a relaxation time that are related to the long time behavior of translational and rotational motion, respectively.

The CG models have force autocorrelation functions that are qualitatively similar to the FG results (with the exception of the time region from about 0.15 to 0.31, which will be discussed further below).

For short times (⁠$t<0.14$⁠), the force correlation function of the CG[1] 11fp model is in good quantitative agreement with the FG results. For such short times, this model is a significant improvement over the MS-CG model, which has no fictitious particles. This is confirmation of a basic principle of the CG[1] parametrization, which is that for times so short that the CG sites can move only a short distance, the fluctuating forces that CG sites feel in an FG system are similar to those they would feel in constrained dynamics.

The CG[2] 2fp model is similar at these short times except that the coarse graining in time has led to less curvature (i.e., less high frequency behavior) of the function. In particular, coarse graining in time has reduced the value of the function at short times and decreased the depth of the minimum of the function.

For $0.14<t<0.3$⁠, the force autocorrelation functions for the MS-CG model and the fictitious particle models are very similar to each other but they differ significantly from the FG results. The possible reasons for this are discussed in Sec. V.

The three CG models have velocity autocorrelation functions that have the qualitative features of the FG results. The only exception is that the CG[2] 2fp velocity autocorrelation function is positive (but small) for times of about 2.5. The way that this model was forced to fit the self-diffusion coefficient of the FG model distorted the shape of the correlation function.

For short times (⁠$t<0.14$⁠) the CG[1] 11fp model has a velocity autocorrelation function that is in good quantitative agreement with the FG model. This is consistent with the fact that for such short times the CG[1] 11fp model has an accurate force autocorrelation function. (The force autocorrelation function is, except for a constant factor, equal to the second time derivative of the velocity autocorrelation function.) However, for longer times the CG[1] 11fp velocity autocorrelation function is much too negative.

For $t>0.5$⁠, the MS-CG velocity autocorrelation function is very similar to the FG function. It is much more accurate than either of the two fictitious particle models for these times. This is the case even though for shorter times the MS-CG model is less accurate than the fictitious particle models. We have no good explanation of this behavior.

The MS-CG model significantly overestimates the mobility of molecules. It has a self-diffusion coefficient that is larger than the FG value by a factor of 2.4 and a rotational relaxation time that is smaller than the FG value by a factor of 2.3.

The CG[1] 11fp model, which is not fit to long time behavior, gives better values for the self-diffusion coefficient and the rotational relaxation time (42% and 19% below the FG values, respectively) than the MS-CG model.

The CG[2] 2fp model has been parametrized to fit the FG results for the self-diffusion coefficient and the orientational relaxation time. As discussed above (Sec. III D), we have done this in a very simple way by adjusting only two parameters. The resulting CG[2] 2fp model has a velocity autocorrelation function whose shape is physically reasonable but not quantitatively accurate. However, for long time simulations, this model would be very useful because it has the correct self-diffusion coefficient and orientational relaxation time. It is conceivable that a more elaborate reparametrization might give a better shape to the velocity autocorrelation function, but we have not investigated this.

In this paper we have presented techniques for enhancing the usefulness of the class of dynamic CG models discussed in a previous paper³² and have applied them to a 2-site CG model of liquid methanol.

The new techniques presented here are as follows:

1. a new type of fictitious particle that exerts fluctuating forces on the center of mass of a molecule without exerting torques on the molecule or influencing vibrational motion;

2. a new algorithm for parametrizing fictitious particle models using Fourier techniques;

3. a systematic way of removing high frequency motion from a CG model, leading in effect to a way of constructing a CG model that is also coarse-grained in time.

The CG[2] 2fp Rigid model is a good candidate for CG simulations of liquid methanol because it has the following properties.

1. It is extremely simple because it has only two fictitious particles on each CG molecule.

2. It is a realistic dynamic model in that it has force and velocity autocorrelation functions that are similar to the functions of the FG model.

   • At very short times (⁠$t<0.14$ ps), both functions are qualitatively correct.

   • At long times both functions decay to zero (to within the statistical error of the simulations) on the correct time scale (approximately 0.4 ps for the force function and approximately 4 ps for the velocity function).

   • However, the model’s results are not in quantitative agreement with the FG data. (This is discussed in more detail below.)

3. The model has the same self-diffusion constant and the same rotational diffusion constant as the FG model and for long times its behavior should be very similar to the FG model. CG[2] fictitious models of the type discussed in this paper (and appropriate generalizations of them) are flexible enough so that when applied to multicomponent systems it will be possible to construct models for which the self-diffusion coefficient and rotational diffusion coefficient of each individual component in the system can be independently adjusted to fit the values of these coefficients as obtained from FG simulations.

The behavior of the correlation functions for the CG[2] 2fp model, with the exception of the very short time behavior of the force autocorrelation function, is qualitatively correct at best. The most striking discrepancy is in the force autocorrelation function for times between about 0.14 and 0.31, see Fig. 3.

Note that in this range of times, the CG[2] 2fp model, the CG[1] 11fp model, and the MS-CG model have very similar force autocorrelation functions that are very different from the FG results. The similarity of the three CG models suggests that the fluctuating forces generated by the fictitious particles are having a very small effect on the force autocorrelation function for this range of times. This is not a result of the coarse-graining in time that was used in constructing the CG[2] 2fp model, since the CG[1] 11fp model, which is much less coarse-grained in time, gives similar results.

It is likely that this discrepancy is due to the inaccuracy of the MS-CG model itself. The MS-CG model was constructed using a variational calculation^42,43 of the CG potential that assumes that the intermolecular part of that potential can be represented as a sum of pairwise additive potentials of interaction between CG particles on different molecules. This is a standard approximation that is reasonably accurate for various CG models of many liquids, and the equilibrium pair correlation functions calculated from simulations of such an MS-CG model are often in good agreement with the FG data for the system. This is the case, for example, for a one-site model of liquid methanol.⁵⁰ However, for the two-site model of methanol, the approximate MS-CG model calculated in this way has pair correlation functions that are significantly different from the FG data.^50,51

Another potential source of discrepancies between fictitious particle models and the FG data is the fact that the simple fictitious particles used here cannot describe time dependent correlations among the fluctuating forces acting on different molecules. It is not clear how important such correlations are for the properties calculated here.

A third potential source of discrepancies is that FG information that is used in a CG[1] parametrization (namely, correlation functions for constrained dynamics) may be informative about the unconstrained dynamics only for systems in which the CG sites are massive enough to be moving very slowly compared to the degrees of freedom that are eliminated in the coarse-graining process.

A fourth potential source of discrepancies is that CG[2] models are constructed by modifying only a limited number of the parameters of a GG[1] model. It might be possible to modify other parameters to improve the model, but that would be an optimization procedure that would not scale well with the complexity of the system, and we have not explored this for two-site methanol.

Nevertheless, as discussed above, by using CG[2] parametrization, an approximate dynamic MS-CG model can be constructed that has a physically reasonable behavior at short and intermediate times and that correctly and quantitatively describes the long time behavior of translational and rotational diffusion.

The ultimate goal of the present work is to develop dynamic CG models for large and complicated systems, including those containing a variety of biomolecules.

As a start toward this goal, we have focussed on simpler systems of smaller molecules to determine what types of features are useful in a dynamic CG model and what types of algorithms are useful.

In the previous paper,³² we discussed several simple systems in which each CG molecule was represented by a single CG particle. These included a mixture of Lennard-Jones particles, a one-site model of ethanol in water, and a one-site model of liquid methanol. The CG[1] and CG[2] models for these systems contained type-1 particles and were very good quantitative models of the corresponding FG systems.

In that paper we also discussed a two-site methanol CG model and a dynamic CG model that incorporated type-1 and type-2 particles. We found that we could not construct a CG[2] version of the model that could fit both the self-diffusion coefficient and the rotational diffusion coefficient. In the present work, using type-2 particles and type-3 particles, we were able to fit both diffusion constants for such a model.

As we move on to construct models with three or more CG sites on a molecule, a new type of particle that couples to the center of mass of the entire molecule would clearly be worthwhile. It will also be necessary to find ways of combining the various types of fictitious particles to get a model that is not only flexible enough to describe the dynamics but also easy to parametrize and easy to simulate. Finally, it is possible to devise fictitious particles that can, in principle, describe hydrodynamic interactions between different CG sites in an approximate way, but it is not yet clear whether a practical algorithm for parametrization of such particles can be developed for the kinds of systems to which the MS-CG method would be applied.

In this paper and the previous one, we have not focussed on the computational efficiency of the individual dynamic CG models. The simple cases we have studied are useful for testing ideas for constructing dynamic models, but they are not examples in which efficiency is an important consideration. However, for larger molecules and more aggressive coarse-graining, it might be possible to devise models that are much more computationally efficient than FG simulations. There are three possible origins for this efficiency.

For highly coarse-grained systems, one CG particle plus a relatively small number of fictitious particles might be sufficient to represent a very large number of atoms.⁵² Since each fictitious particle has simple harmonic interactions with CG particles, and the fictitious particles do not interact with one another, it is likely that the force computation for one time step will be much faster than in the FG simulations.

In some systems, the solvent plays an important role in the way that solutes interact and move, and the solvent degrees of freedom are a large fraction of the degrees of freedom of the system. If the solvent is polar, however, the interactions of the solvent are long ranged and require a significant amount of computation. When such a system is coarse-grained, it is possible to assign no CG sites to any solvent molecule.^53,54 The MS-CG potential then provides an approximate description of the solvent-mediated interaction between solutes, and the fictitious particles provide an approximate representation of the fluctuating forces generated by the solvent. The CG model will then have many fewer degrees of freedom than the FG model, and the force computation for one time step will be much faster than in FG simulations.

For some problems, efficiency of the simulations can be enhanced by taking advantage of the time scales of the various particles. Suppose a CG[1] parametrization of a dynamic CG model is constructed, and then a corresponding CG[2] parametrization is constructed. The fictitious particles could then be divided into three categories:

• those that have negligible effect on the CG degrees of freedom on longer time scales, because they are of such high frequency;

• those that have a significant effect on the CG particles but whose relaxation time is short compared with the time scale of the motion of the CG particles;

• those that have a significant effect on the CG particles but whose relaxation time is not short compared with the time scale of the motion of the CG particles.

Category I fictitious particles can be deleted from the fictitious particle model.

Each category II fictitious particle can be replaced by Langevin forces acting directly on the CG particles that the fictitious particle interacts with.

The result is a Markovian description of the combined set of CG particles and category III fictitious particles. The time step for the simulation of the model can be larger than what would have been necessary for the original CG[2] model since fictitious particles with very high frequency and/or very short relaxation times have been removed from the model.

See supplementary material for the derivation of Eqs. (B3)–(B5) and the details of the evaluation of $Ξ±FG(K,L,t)$ using Eq. (B6).

This research was supported by the National Science Foundation (NSF) through Grant No. CHE-1465248.

The general principles of CG[1] parametrization are discussed in the previous paper. In this section we extend that discussion to include type-3 fictitious particles. Appendix A 1 and A 2 apply to any fictitious particle model that uses only type-2 and type-3 particles. Appendix A 3 is specific for a two-site model of methanol that uses type-2 and type-3 fictitious particles.

Let A(rⁿ) and B(rⁿ) be dynamical variables that are functions of the positions rⁿ of the atoms in the FG system. The fundamental time correlation functions of the constrained FG system that we are concerned with are of the following form:⁵⁵ 

Here rⁿ(t) denotes the atomic positions as a function of time for a trajectory, calculated using constrained dynamics, whose positions and momenta at time 0 are rⁿ and pⁿ, respectively.

This average depends on the positions R^N of the CG sites, which are constrained not to move during the dynamics. If we average over the canonical equilibrium distribution of positions of the CG sites, we end up with another correlation function,

where

is the constrained FG canonical equilibrium distribution function.

Analogous correlation functions, denoted

and

can be defined for constrained CG systems. Note that in a constrained CG system, the CG particle positions R^N are independent of time, but the fictitious particle positions S^M are time dependent.

The only such correlation functions that appeared in the previous paper were the following:

See Eqs. (11), (12), (14), and (15) in the previous paper. Here $δFI(t)$ represents the fluctuating force on CG site (or particle) I, i.e., the difference between the actual force and the force calculated from the coarse-grained potential.

The basic idea of the CG[1] parametrization is to choose the number of fictitious particles and their parameters such that the $ΦIJCG(t)$ correlation functions are as similar as possible to the $ΦIJFG(t)$ correlations, as determined from FG simulations.

As a first step toward parametrizing fictitious particle models with only type-2 and type-3 fictitious particles, we need general expressions for the CG fluctuation functions $ΦIICG(t)$ for all CG particles I in the system and $ΦIJCG(t)$ for all pairs of distinct CG particles, I and J, that are on the same molecule.⁵⁶ (See Sec. III B and Appendix B 3 of the previous paper.³²)

Using the methods discussed in Appendix B 3 of the previous paper, we obtain the following results. Here I and J refer to CG particles in the CG model, α refers to fictitious particles in the CG model, and $sα$ denotes the type (2 or 3) of fictitious particle α,

An explicit expression for $Cα(t)$ is given in Eq. (B7) of the previous paper. These equations imply the following results that are analogous to Eqs. (C1) and (C2) of the previous paper:

Here

For the case of two-site methanol, all molecules have two distinct sites and have equivalent properties. Let K and L denote the methyl and hydroxyl CG particles, respectively, on the same molecule. All fictitious particles of type-3 have the same value of $λα$⁠, namely, M_K/(M_K + M_L). Then it is straightforward to show that

Using these results, Eqs. (A2) and (A3) for all molecules are the same and can be written in the following form:

where

The functions Σ(2, K, L, t) and Σ(3, K, L, t) describe the fluctuating forces generated by the fictitious particles of type-2 and type-3, respectively. The time dependence of these functions is determined by the number of fictitious particles of each type and the parameters of the fictitious particles. An explicit expression for $Cα(t)$ is given in Eq. (B7) of the previous paper. Thus we have explicit analytic expressions for these functions. (Explicit expressions for the Fourier transform of $Cα(t)$ can be derived, so we can obtain explicit analytic expressions for the Fourier transforms of the Σ functions.)

It is clear that in general there is no choice of the two Σ functions that can generate the three functions on the left of Eqs. (A4)–(A6) that are in agreement with the corresponding $ΦFG(t)$ functions. The three distinct $ΦCG(t)$ functions are clearly linearly related, but there is no reason to think that the same is the case for the three distinct $ΦFG(t)$ functions. So the most we can accomplish in parametrizing this fictitious particle model is to satisfy two distinct equalities relating $ΦCG(t)$ and $ΦFG(t)$ functions.

We define the following fluctuation force time correlation functions for the cFG and cCG system:

Using Eqs. (A4)–(A7), these functions for the cCG system can be expressed as

Criterion for parametrization of fictitious particles. We choose the number of fictitious particles and their parameters so that

In practice, we work in the Fourier domain and require that

Here a caret denotes the cosine transform of the function. Using Eqs. (A8) and (A9), as well as Eq. (B7) of the previous paper,³² it is straightforward to obtain explicit analytic expressions for the $Σ^$ functions in these equations in terms of the number of fictitious particles and their parameters.

The method for calculating the $Ξ±FG(t)$ functions from constrained FG simulations is discussed in Appendix B 1. The method for choosing the number of fictitious particles and the values of the parameters that lead to Σ functions that satisfy Eqs. (A13) and (A14) is discussed in Appendix B 2.

The $Ξ±FG$ correlation functions on the right of Eqs. (A13) and (A14) are defined in Eq. (A10) as time correlation functions for a constrained FG canonical ensemble discussed in the previous paper³² (see Sec. II C and Appendix A of that paper). The time dependence of such correlation functions can be calculated from constrained simulations of systems whose initial states are sampled from this ensemble (see Appendix A 2 d of the previous paper for the algorithm for performing such simulations).

As discussed in the previous paper, for FG systems, the total force on CG site I, $FI(t)$⁠, can be written as the sum of two terms: the MS-CG force F_I(R^N) and the fluctuation force $δFI(t)$⁠,

The fluctuation force on individual sites $δFI(t)$ appears in Eq. (A10). The calculation of such time correlation functions using simulation data is made difficult by the fact that although the total forces $FI(t)$ on CG sites are calculated during the simulations, very accurate values of the MS-CG forces F_I(R^N) for a constrained FG trajectory are not known during the simulation, and hence, the fluctuation forces $δFI(t)$ are not known.

To address this problem, we consider the corresponding correlation functions of the total force,

A straightforward but detailed argument in the supplementary material leads to the following relationships:

Hence

Eq. (B6) expresses the $Ξ±FG(K,L,t)$ functions in terms of quantities that can be calculated from FG trajectories. The first term on the right can be evaluated straightforwardly. The second term on the right is more difficult to estimate and is very sensitive to statistical error in the FG simulation data. In the supplementary material we describe the method that we use to determine the second term and then construct useful representations of $Ξ±FG(K,L,t)$ using the simulation data.

The algorithm that we used for optimization of the number of fictitious particles and their parameters is the following.

1. Find the $Ξ±FG(t)$ functions using the numerical method described in Sec. III B 1.

2. For each function calculate the Fourier transforms $Ξ^+FG$ and $Ξ^−FG$ numerically.

3. For each function, choose a reasonable value of the largest frequency to be considered in the fitting process.

4. Choose a reasonable number of type-3 fictitious particles and choose a reasonable initial set of $μα$⁠, $ωα$⁠, and $γα$ parameters.

5. Vary the parameters of the type-3 particles so that $Σ^(3,K,L,ω)$ satisfies Eq. (A13) in a least squares sense.

6. Increase the number of type-3 fictitious particles and repeat step 5 until $Σ^(3,K,L,ω)$ satisfies Eq. (A13) with sufficient accuracy.

7. Use the optimized parameters for the type-3 fictitious particles to calculate $Σ^−(3,K,L,ω)$⁠.

8. Choose a reasonable number of type-2 fictitious particles and a reasonable set of initial parameters.

9. Vary the parameters of the type-2 particles so that $Σ^(2,K,L,ω)$ satisfies Eq. (A14) in a least squares sense.

10. Increase the number of type 2 fictitious particles until $Σ^(2,K,L,ω)$ satisfies Eq. (A14) with sufficient accuracy.

The results of this algorithm for the current problem of two-site methanol are given in Fig. 1.

The result of the process described just above in Appendix B 2 is the number of fictitious particles of type-2, the number of fictitious particles of type-3, and the numerical values of the fundamental parameters of the fictitious particles. All of this information is needed to perform simulations of the dynamic CG model.

The numerical values of the parameters can be useful for deciding which fictitious particles are very important for the accuracy of the model and which are of less importance and can perhaps be deleted. They also give information about the time scales associated with each fictitious particle, which can be used in deciding which fictitious particles to delete when constructing models that are coarse-grained in time, as discussed in Sec. III C.

In this section, we discuss various properties of fictitious parameters.

Each fictitious particle α has a mass $μα$ and a harmonic oscillator frequency $ωα$ that appear in the Hamiltonian in Eq. (2). The Langevin damping $γα$ is a parameter with units of (time)⁻¹. The form of the Langevin equation for one component of the coordinate of a fictitious particle is

$γα−1$ is thus the relaxation time of the velocity of the particle due to Langevin forces.

The three parameters $(μω,ωα,γα)$ appear in the expressions above for the functions Σ(2, K, L, t) and Σ(3, K, L, t), see Eqs. (A8) and (A9). Each of these two functions has a separate additive contribution from each fictitious particle that interacts with the CG particles K and L. These functions are linearly related to the functions $Ξ±CG(K,L,t)$ that are central to the CG[1] parametrization method in that they are forced to fit the $Ξ±FG(K,L,t)$ functions obtained from constrained FG simulations.

As a result, when the optimization process is complete, the following equation holds:

Thus, each fictitious particle makes a separate additive contribution to the function that is used to fit $Ξ+FG(K,L,t)$ and $Ξ−FG(K,L,t)$⁠, the constrained FG correlation functions of the fluctuating force.

It follows that the relative importance of a fictitious particle in fitting the FG data can be judged from the magnitude of its contribution to Σ(2, K, L, t) for type-2 particles or to Σ(3, K, L, t) for type-3 particles. There are three useful things to know about the contribution of an individual fictitious particle to the Σ function for particles of its type.

1. t0. The t0 of a fictitious particle is defined as its contribution to the value of the Σ function of its type at t = 0. This determines its contribution to the t = 0 values of the $Ξ$^FG functions. This in turn is a measure of its importance to fitting the mean squared fluctuating forces in the constrained FG system.

2. ff0. The ff0 of a fictitious particle is defined as its contribution to the value of the zero frequency Fourier transform of the Σ function of its type. This is a measure of its importance to fitting the integral of the fluctuating force correlation functions in the constrained FG system. This is a measure of its contribution to the total friction on the CG degrees of freedom exerted on them by the degrees of freedom that are not retained in the coarse-graining process.

3. aof. The aof of a fictitious particle is its approximate oscillation frequency. It is defined as $(ω2−γ2/4)1/2$ for underdamped oscillators and 0 for overdamped oscillators. It is an approximation to the frequency of the peak of the Fourier transform of its contribution to the Σ function of its type.

All-atom simulations of liquid methanol containing 1000 molecules were performed with the LAMMPS⁵⁷ package, using OPLS-AA force field⁴⁹ and the velocity Verlet algorithm. The Particle-Particle Particle-Mesh (PPPM) method⁵⁸ was used for long range electrostatic interactions.

The system was equilibrated at 300 K and constant volume of 68 417.9 ų.

We performed ten simulations of the constrained FG system in order to calculate the cFG canonical ensemble average fluctuation force functions required to parametrize the CG[1] model (see Appendix A 2 d of the previous paper³² for a discussion of the algorithm for performing the ensemble average). In these simulations the center of mass of each CG site was constrained not to move and no temperature control was applied. Each run was 100 ps long. The force fluctuation functions $Ξ±FG(t)$ were calculated based on these data.

We also performed a number of unconstrained simulations, up to 1 ns long, to calculate the unconstrained force and velocity autocorrelation functions as well as the self diffusion constant and the orientational relaxation time for the methanol molecules.

The time step used for the constrained and unconstrained FG simulations was 1 fs.

Simulations of the various dynamic CG models were also performed with the LAMMPS⁵⁷ package. The tabulated potentials derived from MS-CG method were used as the interactions between CG sites. The interactions between CG sites and the fictitious particles were computed using custom angle potentials. All simulations of dynamic CG models were performed using Hamiltonian dynamics for the CG and fictitious particles and additional Langevin forces on the fictitious particles. The time step for the CG simulations was 2 fs.

We did two types of simulations.

1. 100 ps long simulations were performed with data collected for every time step. Those simulations were used for calculations of molecular center of mass force and velocity autocorrelation functions and the rotational diffusion constant.

2. 1 ns long simulations were performed with data collected every 1000 time steps. Those simulations were used to calculate the translational self diffusion constant.