Comparison of Modern Langevin Integrators for Simulations of Coarse-Grained Polymer Melts
CCOMPARISON OF MODERN LANGEVIN INTEGRATORS FORSIMULATIONS OF COARSE-GRAINED POLYMER MELTS
JOSHUA FINKELSTEIN, GIACOMO FIORIN, AND BENJAMIN SEIBOLD
Abstract.
For a wide range of phenomena, current computational ability does not alwaysallow for fully atomistic simulations of high-dimensional molecular systems to reach timescales of interest. Coarse-graining (CG) is an established approach to alleviate the impact ofcomputational limits while retaining the same algorithms used in atomistic simulations. It isof importance to understand how algorithms such as Langevin integrators perform on non-trivial CG molecular systems, and in particular how large of an integration time step can beused without introducing unacceptable amounts of error into averaged quantities of interest.To investigate this, we examined three different Langevin integrators on a CG polymer melt:the recently developed BAOAB method by Leimkuhler and Matthews [17], the Grønbech-Jensen and Farago method [12], or G–JF, and the frequently used Br¨unger-Brooks-Karplusintegrator [6], also known as BBK. We compute and analyze key statistical properties foreach. Our results indicate that the three integrators perform similarly when using a smallfriction parameter; however, outside of this regime the use of large integration steps producessignificant deviations from the predicted diffusivity and steady-state distributions for allintegration methods examined with the exception of G–JF. Introduction
A central obstacle in using molecular dynamics (MD) simulations for quantitative predic-tions in material science and molecular biology is the presence of a wide range of time scalesthat are not well-separated. To investigate phenomena in such systems that occur over largetime scales, while accurately resolving smaller ones, a large number of integration steps isrequired. For example, in fully atomistic simulations, the size of the integration time step isconstrained by the fastest physical time scales and is typically on the order of femtoseconds.It becomes computationally inefficient, yet necessary, to use these relatively small time stepsfor integrating the medium to long-range time scale portion of the force field.There exist several methods for working around this bottleneck. One such method isNose-Hoover chains (with or without RESPA) [36, 37, 22]. This allows one to match thetime step size for each degree of freedom to its corresponding time scale and compose theresulting operators via a symmetric Trotter splitting. Unfortunately, the Nose-Hoover chainequations result in the introduction of many undetermined parameters and, depending onthe complexity of the system being simulated, may require different atoms to be coupledto different thermostats, thus complicating implementation. Other deterministic approacheswhich attempt to constrain the fast degrees of freedom are also used [30].A Langevin thermostat is a simple and efficient way of simulating constant temperatureconditions and has the intuitive physical interpretation of describing a molecular system in
Mathematics Subject Classification.
Key words and phrases.
Langevin integrator, coarse-grained, polymer melt, molecular dynamics.PACS numbers: 31.15.xv, 31.15.at, 36.20.Ey. a r X i v : . [ phy s i c s . c o m p - ph ] A p r J. FINKELSTEIN, G. FIORIN, AND B. SEIBOLD the presence of an implicit solvent or heat bath. The interaction between the heat bath andthe system is collapsed into the friction parameter γ , thus avoiding the need to representthe heat bath as a set of particles altogether. In the context of coarse-grained (CG) modelswith explicit solvent, the missing microscopic degrees of freedom can be described by suchstochastic forces, at least to a first approximation. For example, the dissipative-particle-dynamics (DPD) methodology [28] relies on this description to systematically build coarse-grained models with low computational cost.Because improvements in computer hardware are being introduced more slowly in recentyears, the usefulness of MD models that run at speeds higher than atomistic simulations willalso increase. As CG models become more commonplace, there is a need to systematicallyunderstand the numerical performance and the limitations of current temporal integrationschemes on CG systems, avoiding the reliance on rules of thumb that were derived primarilyfor atomistic simulations.In this work we examine methods used to numerically solve the Langevin equation, such asthe ones by Leimkuhler and Matthews [17] and Grønbech-Jensen and Farago [12], known asG–JF, with particular concern toward their respective sampling properties in CG simulations.Among the family of integrators described by Leimkuhler and Matthews [17], the BAOABmethod is the one characterised by the smallest configurational sampling error, and the onlyone here considered. Both the G-JF and BAOAB methods are weakly second-order accurate[3, 17] and produce the exact configurational mean, variance and co-variance of the harmonicoscillator. This important property is not produced by many other Langevin schemes [38, 26].To compare the two schemes with a representative of more traditional integration methods,our analysis includes the well-established Br¨unger-Brooks-Karplus method [6], or BBK. Outof the many formulations of BBK [19] we chose one (indicated here as BBK ∗ ) that performswell for the system under study: comparisons to the classical formulation are also made.Several numerical studies have been conducted on the performance of these integrators inatomistic simulations (e.g. see [13, 33]). However, CG models tend to use smoother potentialenergy functions than their fully atomistic counterparts, and thus allow for much larger inte-gration time steps. A key question that arises then is: how large can the time step h be madewithout introducing unacceptable levels of error into averaged static and dynamic quantities?In [4], G–JF is used to simulate a CG lipid bilayer in implicit solvent and averaged energyterms (both potential and kinetic) are examined: however, the system was simulated forrelatively short MD trajectories ( < γh values are explicitly required in its derivation. Our present work considersa wide range of values for γh . To date, we are not aware of any CG studies for BAOAB.There are several choices of CG-model resolutions to choose from: a recent survey ofseveral CG models [31] suggested that an adequate representation of the phase behavior seenin atomistic simulations is given by models of polyethylene chains with three or four methylenegroups per CG particle. The model by Klein and coworkers [34, 35], the MARTINI [21] andSalerno-Grest [31] models are significant examples of this level of resolution. Many othersuch models also exist [8, 11, 16, 5, 25, 24, 14]. It is also recognised that technologicalrequirements will motivate further effort to develop accurate CG models at lower levels ofresolution (mapping into fewer CG particles for the same system). Transitioning into suchmodels affects significantly the balance between Hamiltonian, stochastic and inertial terms OMPARISON OF LANGEVIN INTEGRATORS FOR COARSE-GRAINED POLYMER SIMULATIONS 3 in the equations of motion, and may require the use of existing Langevin methods outside oftheir typical range of parameter values.In fact, Langevin parameters near the high-friction limit are of particular concern. Here,the acceleration may be neglected and the second-order integrator can safely be substitutedwith a first-order one for improved numerical stability. However, due to large differences inhow existing CG models are formulated, as well as differences between potential energy termsin the same model, it is far from unusual to see applications of a Langevin thermostat thatapproach this region at least transiently. Unfortunately, any resulting biases in sampling areoften difficult to detect due to the heterogeneous nature of the physical system examined, orits proximity to a phase transition.Here we used as benchmark system a polyethylene melt (C ) modeled with three methylenegroups per CG particle. All three schemes, BAOAB, BBK and G–JF, were compared for awide-range of friction parameter values and time step sizes by examining relevant statisticalquantities from the simulations. The key finding of our study is that in the high-friction( γ ≈ . − ) regime, the G–JF method performs measurably better than BAOAB and BBKin reproducing molecular diffusivity and configurational distributions. The results obtainedprovide indications that Langevin integrators with similar properties to G–JF should beconsidered for use in CG simulations that aim to preserve dynamic properties and stationarydistributions equally accurately. Though, diffusivity notwithstanding, BAOAB and G–JFsample equally well the configurational distributions of the system considered here.2. Background and Theory
We consider an N -particle system with potential energy U , immersed in a heat bath withthe constant temperature T , modeled by the Langevin equation: d Q = M − P dt ,d P = − M − ∇ U ( Q ) dt − γ M − P dt + σ M − / d W . (1)Here σ = √ k b T γ is the noise coefficient, k b is Boltzmann’s constant, γ is the (spatiallyindependent) collision rate parameter (measured in units of fs − ), M a diagonal mass matrix, W is 3 N -dimensional Brownian motion, and H ( Q , P ) = P T M − P + U ( Q ) is the Hamiltonian.The Langevin equation is a stochastic differential equation (SDE), so we use capital lettersto remind us that position and velocity are stochastic processes.Usually in MD, it is not the exact dynamics generated by (1) that are of interest, but ratheran accurate sampling of the equilibrium distributions in phase space. Assuming that H ( q, p )is such that e −H ( q,p ) is integrable, one expects (1) to be ergodic and have the stationarydistribution µ ( dqdp ) = Q − e −H ( q,p ) /k b T dqdp , where Q is a normalisation constant and µ is the Boltzmann-Gibbs, or canonical, distribution.For realistic MD potentials, such as Lennard-Jones and/or Coulombic interaction forces, thesolution to (1) needs to be approximated numerically. Moreover, in many applications, such asthe one here, U is non-globally Lipschitz and singular. Consequently, many standard resultsin SDE theory do not apply, thus limiting the possibilities of a complete formal analysis ofnumerical schemes for (1). One must therefore study these schemes computationally.The fidelity to which numerical schemes reproduce the Boltzmann-Gibbs distribution µ is our principal interest. It is the hope that infinite time-averaged observables obtained J. FINKELSTEIN, G. FIORIN, AND B. SEIBOLD from these numerical schemes would reproduce correct statistical averages with respect to µ . However, even with ergodicity typically assumed, the steady states produced by eachnumerical scheme, µ G–JF , µ BAOAB , and µ BBK , respectively, will in general differ from µ , anddepend on the friction parameter γ and time step h . Therefore we expect two distinct sourcesof error in the calculation of distributions and observables: use of a finite trajectory instead ofan infinite one and the error due to the numerical scheme’s steady state distribution differingfrom the canonical distribution.Although, in general, statistical averages generated by numerical approximations of (1)cannot be derived analytically, exact formulas for mean, variance and correlation can be cal-culated in the flat and harmonic potential case, which is enough to characterise any stationarydistribution when starting with Gaussian initial conditions.2.1. Numerical Methods Studied.
In this work, we examined three different Langevinintegration schemes: G–JF, BAOAB, and the Br¨unger-Brooks-Karplus method [6], known asBBK, on a CG polymer melt. Both G–JF and BAOAB are included as options in LAMMPSand NAMD, respectively [1, 2]. The G–JF thermostat is a stochastic two-stage partitionedRunge-Kutta method [7, 3] and was shown to have highly desirable configurational properties[12]; particularly, Einstein’s diffusion relation holds exactly and the configurational averagesfor the harmonic oscillator are independent of both the time step h and the friction parameter γ . As described in [17], BAOAB is but one of many splitting schemes obtained by composingsolution operators in various orderings which evolve the A, B and O portions of the Langevinvector field (1): (cid:18) d Q d P (cid:19) = (cid:18) M − P (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) A dt + (cid:18) − M − ∇ U ( Q ) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) B dt + (cid:18) − γ M − P dt + σ M − / d W (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) O . BAOAB is the result of taking half steps for B and then A, a full step for O and then halfsteps again for A and then B.Similar to G–JF, BAOAB also reproduces exact sampling for the harmonic oscillator. More-over, it possesses an additional favorable configurational sampling property, termed “super-convergence” [17]: when γ is sufficiently large, the leading order error terms of averaged phasespace quantities will exhibit a 4th order error scaling in h for typical time step values, therebyyielding more accurate averages with the same computational effort. This property providessome motivation as to why BAOAB is of interest to practitioners.BBK has been a well-known Langevin discretisation method for the last three decades andis the default Langevin integrator in the popular MD suite, NAMD [2]. Similar to BAOAB,BBK is also a splitting method. It is weakly first order accurate [26] and in the free particlecase reproduces the Einstein relation. However, exact statistics are not recovered in the caseof the harmonic oscillator, in sharp contrast to G–JF and BAOAB. Equation (1) is oftenre-formulated in terms of position and velocity, instead of momentum: d Q = V dt ,d V = − M − ∇ U ( Q ) dt − γ V dt + σ M − / d W . (2)This form will serve as the governing equation for our forthcoming analysis and discussion.Set a := (1 − γh/ γh/ − and b := (1 + γh/ − , with h being the time step used for OMPARISON OF LANGEVIN INTEGRATORS FOR COARSE-GRAINED POLYMER SIMULATIONS 5 discretisation. Note that, for γh sufficiently small, one has a := 1 − γh/
21 + γh/ e − γh + O (( γh ) ) , which leads to a n = e − γt + O (( γh ) ) when nh = t . The quantities a and a n appear severaltimes in the subsequent paragraphs and sections.We start by displaying the recursion formulas for the considered numerical schemes dis-cretising (2). For a fixed time step h > Q , V ), the position andvelocity at time t = nh for each scheme are displayed below. Define (cid:101) σ = (cid:112) k b T (1 − e − γh ).Then, the G–JF update rule is: Q n = Q n − + bh V n − − bh M − ∇ U ( Q n − ) + bσh / M − / ξ n − , V n = a V n − − h M − ( a ∇ U ( Q n − ) + ∇ U ( Q n )) + bσ √ h M − / ξ n − , (3)the BAOAB method is: Q n = Q n − + h e − γh ) V n − − h e − γh ) M − ∇ U ( Q n − ) + ˜ σh M − / ξ n − , V n = e − γh V n − − h M − ( e − γh ∇ U ( Q n − ) + ∇ U ( Q n )) + ˜ σ M − / ξ n − , (4)and finally the BBK method is: Q n = Q n − + h (1 − γh/ V n − − h M − ∇ U ( Q n − ) + σh / M − / ξ n − , V n = a V n − − bh M − ( ∇ U ( Q n − ) + ∇ U ( Q n )) + bσ √ h M − / ( ξ n − + ξ n ) . (5)In its original formulation [6], the BBK numerical scheme is given for only the position, leavingsome ambiguity as to how the velocities are defined. Using the second order approximation V n ≈ ( Q n +1 − Q n − ) / h , the splitting formulation of BBK for both position and velocity isobtained (e.g. as seen in [15]). This is the same substitution one employs in transforming theposition-only Verlet integrator to velocity Verlet: V n − / = V n − + h M − (cid:0) − ∇ U ( Q n − ) − γ MV n − + σ √ h M / ξ n − (cid:1) , Q n = Q n − + h V n − / , V n = V n − / + h M − (cid:0) − ∇ U ( Q n ) − γ MV n + σ √ h M / ξ n (cid:1) . (6)It is then a simple exercise to derive (5) from (6). In particular, this method requires twoindependent random variables ξ n − and ξ n where ξ n is then re-used in the next step. However,[19] suggests the use of several variations of BBK which vary in how these random variables areselected. One such variation, which we denote by BBK ∗ , is obtained by taking ξ n − = ξ n in(6), and not conducting any re-use in the next step. This BBK ∗ variant is not equivalent to theoriginal version of BBK; and in fact, numerical tests indicated better all-around performancewith our particular CG molecular system of interest in the commonly used regime of γh ≤ . J. FINKELSTEIN, G. FIORIN, AND B. SEIBOLD
Analytical Properties of the Methods in One Dimension.
Before moving to thecomputational results, we first summarise the key statistical properties for each of the threeschemes in one dimension. This section attempts to highlight some key structural differences(and similarities) of the three schemes. Some schemes reproduce certain statistical quantitiesexactly whereas others reproduce such quantities only approximately.We consider the standard examples of a single particle diffusing in a heat bath, and astandard harmonic oscillator, modeling for instance a covalent bond between two particles.In these cases, we can directly compare statistical quantities generated by the numericalschemes with those generated by the true analytical solution of (2). Though simplistic, theseexamples illustrate some important properties. Most of the calculations for these exampleshave, in parts, been previously exposited [7, 12, 26, 17, 38, 18].2.2.1.
Harmonic Potential.
We consider a potential function of the form U ( Q ) = ωQ / ω >
0, so that H ( Q, P ) = P / m + ωQ /
2, and (2) becomes dQ = V dt ,dV = − ωQdt − γV dt + σm − / dW . (7)The following calculations show that the variance of position for the BAOAB and G–JFintegrators is independent of γ and h . This is not true for BBK. For the linear system(7), stationary distributions can be analytically derived for the three methods, denoted by µ BAOAB , µ BBK ∗ , and µ G–JF , respectively. We re-write (3) applied to (7) into matrix form: (cid:20) Q n +1 V n +1 (cid:21) = (cid:34) − bωh m bh − hωbm (1 − h ω m ) a − h ωb m (cid:35) (cid:20) Q n V n (cid:21) + (cid:34) bσh / √ mbσ √ h √ m (1 − h ω m ) (cid:35) ξ n . This is a two-dimensional ergodic Markov chain with a unique stationary measure, µ G–JF . Wecan then calculate a corresponding matrix equation for Q n +1 , V n +1 and Q n +1 V n +1 . Takingexpectations on both sides of this equation, and taking n → ∞ , yields a subsequent 3 × E ( Q ∞ ), E ( V ∞ ) and E (( QV ) ∞ ), the vector of steady-state averages. Solvingthe resultant linear system yields the G–JF stationary distribution: µ G–JF ( dqdv ) ∝ exp (cid:18) − β (cid:18) mv − h ω m ) + ωq (cid:19)(cid:19) dqdv . Using (4) and (5) one can derive analogous linear equations and expressions for BBK ∗ andBAOAB, as follows. • BBK ∗ : (cid:20) Q n V n (cid:21) = (cid:34) − h ω m h (1 − γh ) − hbωm (cid:0) − h ω m (cid:1) a (1 − h ω m ) (cid:35) (cid:20) Q n − V n − (cid:21) + (cid:34) σh / √ m ξ n − σb √ h √ m (1 − h ω m ) ξ n − (cid:35) (8)and µ BBK ∗ ( dqdv ) ∝ exp (cid:18) − β (cid:18) mv − h ω m )(1 − γh/ ωq (cid:19)(cid:19) dqdv . OMPARISON OF LANGEVIN INTEGRATORS FOR COARSE-GRAINED POLYMER SIMULATIONS 7 • BAOAB: (cid:20) Q n +1 V n +1 (cid:21) = (cid:34) − (1 + e − γh ) h ω m (1 + e − γh ) h − (1 − h ω m )(1 + e − γh ) hω m e − γh (1 − h ω m ) (cid:35) (cid:20) Q n V n (cid:21) + h (cid:113) k b Tm (1 − e − γh )(1 − h ω m ) (cid:113) k b Tm (1 − e − γh ) ξ n , (9)and µ BAOAB ( dqdv ) ∝ exp (cid:18) − β (cid:18) mv − h ω/ m ) + ωq (cid:19)(cid:19) dqdv . The variances and covariances for position and velocity are listed in Table 1. These expressionsshow the dependence on the dimensionless quantity γh and the time step h . In particular,both BAOAB and G–JF produce the exactly correct configurational variance and co-variance.Method (cid:104) Q n (cid:105) h,γ (cid:104) V n (cid:105) h,γ (cid:104) Q n V n (cid:105) h,γ Exact k b Tω k b Tm k b Tω k b Tm (1 − h ω/ m ) 0BBK ∗ k b Tω (1 − h ω m )(1 − γh/ − k b Tm k b Tω k b Tm (1 − h ω/ m ) 0 Table 1.
Numerical stationary averages, as functions of h and γ , for the three numerical methodsapplied to the one-dimensional harmonic oscillator. Both G–JF and BAOAB are exact for positionvariance, while BBK ∗ is not. Thermal Diffusion.
The simplest possible case for (2) is when F ≡
0, i.e., free diffusion.Albeit simple, it is insightful to understand the behavior of the integrators in this case. Inthermal diffusion, (2) reduces to dQ = V dt ,dV = − γV dt + σm − / dW , which can be solved analytically. The velocity V t is an Ornstein-Uhlenbeck process and hassolution: V t = e − γt V + σ √ m (cid:90) t e − γ ( t − s ) dW s . (10)The particle position, Q t , is then Q t = Q + 1 m (cid:90) t V u du = Q + γ (1 − e − γt ) V + σ √ m (cid:90) t (cid:18) (cid:90) u e − γ ( u − s ) dW s (cid:19) du . (11) J. FINKELSTEIN, G. FIORIN, AND B. SEIBOLD
Equation (11) is used to find the mean position and mean squared position. Assuming E ( Q t ) = 0, E ( Q t ) = E ( Q ) + 2 =0 (cid:122) (cid:125)(cid:124) (cid:123) E ( Q ) E (cid:18) (cid:90) t V u du (cid:19) + E (cid:18) (cid:90) t V u du (cid:19) = E ( Q ) + 1 γ (1 − e − γt ) E ( V ) + 2 D (cid:18) t − γ (1 − e − γt ) + 12 γ (1 − e − γt ) (cid:19) , where D := k b T /mγ is the diffusion coefficient. The large time asymptotic behavior of themean squared position for mean zero initial position is then V ( Q t ) ∼ Dt . (12)The velocity auto-correlation function and the covariance can also be computed from (10)and (11). These quantities are displayed in Table 2.2.2.3.
Diffusive Behavior of the Numerical Schemes.
We set ω = 0 in (8)–(9) to obtain theupdate rules for each numerical scheme in the zero potential case. A key quantity of interestis the mean square displacement for the particle position. For BBK ∗ , Q n = Q n − + h (1 − γh/ V n − + σh √ m ξ n − ,V n = aV n − + bσ √ h √ m ξ n − . (13)Given a fixed time step size h >
0, iterate the above recursive formula backwards to writethe position as a finite sum of independent Gaussians and the initial conditions: Q n = Q + (1 − γh ) h (1 − a n +1 )1 − a V + σh √ m n − (cid:88) k =0 ξ k (cid:18) a (1 − a k +1 )1 − a + 12 (cid:19) . Therefore, V ( Q n ) = V ( Q )+(1 − γh/ h (cid:18) − a ( n +1) − a (cid:19) V ( V )+ 2 D (cid:18) t − γ (1 − a n )(1 − γh/ + 12 γ (1 − a n )(1 − γh/ (cid:19) . where t = nh . Sending n → ∞ shows that the scheme preserves the Einstein diffusion relationin the limit. As with BBK ∗ , G–JF preserves the Einstein diffusion relation in the long timelimit as well. The position and its mean square displacement at time t = nh are: Q n = Q + bh (1 − a n +1 )1 − a V + σbh / √ m n − (cid:88) k =0 ξ k (cid:18) (1 − a k +1 ) b − a + 12 (cid:19) , V ( Q n ) = V ( Q ) + b h (1 − a ( n +1) ) (1 − a ) V ( V ) + 2 D (cid:18) t − γ (1 − a n ) a + 12 γ (1 − a n ) a (cid:19) , which has the same limiting behavior as in [12]. Additionally, both BBK ∗ and G–JF generatethe correct steady-state behavior for velocity. For a fixed time step h , the BAOAB scheme OMPARISON OF LANGEVIN INTEGRATORS FOR COARSE-GRAINED POLYMER SIMULATIONS 9 approximates thermal diffusion via: Q n = Q n − + h (1 + e − γh ) V n − + h (cid:113) k b T (1 − e − γh ) m ξ n − V n = e − γh V n − + (cid:113) k b T (1 − e − γh ) m ξ n − . An important distinction here is that V n +1 is given by the exact Ornstein-Uhlenbeck flow(in law), whereas the velocity updates are only approximate for BBK ∗ and G–JF. As donewith BBK ∗ and G–JF, by iterating backwards, we can write the position as a finite sum ofindependent identically distributed random variables: Q n = Q + h (cid:18) e − γh − e − γh (cid:19) (1 − e − γt ) V + h k b Tm (1 − e − γh ) n − (cid:88) k =0 ξ k (cid:18) e − γh − e − γh (1 − e − γh ( k +1) ) + 1 (cid:19) . (14)Then V ( Q n ) is V ( Q ) + γh (cid:18) e − γh − e − γh (cid:19) (1 − e − γt ) V ( V )+ 2 D (cid:18) γh e − γh ) − e − γh t − γh − e − γh )(1 + e − γh )(1 − e − γh ) (1 − e − γt ) e − γh + 12 γ (1 − e − γt ) e − γh (cid:19) . This appears to be vastly different than the other schemes, but one can check that (12) isrecovered when sending γh →
0. So whenever γ and h are fixed such that γh is sufficientlysmall, lim n →∞ V ( Q n ) /nh ≈ D , and BAOAB produces an acceptable approximation to thecorrect diffusive behavior. More particularly, the time evolution of the variance for BAOABat large times evolves according to 2 (cid:101) Dt where (cid:101) D = D (cid:18) γh e − γh ) − e − γh (cid:19) , is the effective BAOAB diffusion coefficient, showing how the calculated diffusion for BAOABdeviates from theory when γh is sufficiently large. The quantity inside the parentheses tendsto 1 as γh → γh is increased. Key statistical quantities foreach numerical scheme are tabulated in Table 2 for simple initial conditions. In fact, among allmethods within the aforementioned A,B,O family of splitting schemes, when a single randomsample per time step is desired, the incorrect diffusive behavior elucidated above is universal.In the free particle case, B induces the identity operator, so that the number of possiblelettered combinations reduce to AO, OA, and AOA. A quick calculation reveals that none ofthese methods reproduces the Einstein relation.3. Computational Methodology
In our computational study we considered a collection of 128 poly-ethylene chains (C H ),simulated at a fixed temperature of 450 K (i.e., well above the melting point) in a box withside lengths of 58 .
065 ˚A and periodic boundary conditions. This resulted in a density of0 . / ˚A = 0 . Method (cid:104) Q n (cid:105) h,γ (cid:104) V n (cid:105) h,γ (cid:104) Q n V n (cid:105) h,γ C h,γ ( nh )Exact ∼ Dnh k b T /m (1 − e − γhn ) D e − γhn k b Tm G–JF ∼ Dnh k b T /m (1 − a n ) D a n k b Tm BBK ∗ ∼ Dnh k b T /m (1 − a n ) D (cid:0) γh (cid:1) a n k b Tm BAOAB ∼ ( γh ) (cid:0) e − γh − e − γh (cid:1) Dnh k b T /m (1 − e − γhn ) D (cid:0) γh e − γh − e − γh (cid:1) e − γhn k b Tm Table 2.
The analytic and numerical schemes’ averages in the case of Brownian motion (Einsteindiffusion) with δ –distributed Q and Maxwell-Boltzmann-distributed V . The averages for thenumerical schemes are given as functions of h and γ . Both G–JF and BBK ∗ are exact for positionvariance, while BAOAB is not. C h,γ ( nh ) is the velocity autocorrelation function of the numericalscheme with time step h and friction parameter γ at time nh . directly from liquid-phase physical properties [35]. The hydrocarbon chains are modeled incoarse-grained resolution, where the − CH2CH2CH2 − and CH3CH2CH2 − groups are mappedto spherical “beads”, called CM and CT, respectively [35]. Compared to atomistic resolution,this reduces the system from a total of 18,432 atoms to 2,048 CG beads. Initial positionsof the CM and CT particles are taken from their respective centers-of-mass, and are evolvedaccording to the following interaction energy: U ( q ) = (cid:88) a ∈ angles k a ( θ − θ ) + (cid:88) b ∈ bonds k b ( r − r ) + U ( α,β ) LJ ( q ) , where ( α, β ) ∈ { (CT , CT) , (CT , CM) , (CM , CM) } and U ( α,β ) LJ is the 9–6 Lennard-Jones poten-tial [35] U ( α,β ) LJ ( q ) = N (cid:88) i (cid:54) = j (cid:15) α,β (cid:18)(cid:18) σ α,β | q i − q j | (cid:19) − (cid:18) σ α,β | q i − q j | (cid:19) (cid:19) · {| q i − q j | <δ } , with (cid:15) CT,CT = 0 .
42 kcal/mol and σ CT,CT = 4 .
506 ˚A for the CT − CT interaction, (cid:15)
CT,CM = 0 .
444 kcal/mol and σ CT,CM = 4 . − CM interaction, and (cid:15)
CM,CM =0 .
469 kcal/mol and σ CM,CM = 4 .
585 ˚A for the CM − CM interaction. Here δ represents theLennard-Jones cut-off distance of 15 ˚A. The function 1 {| q i − q j | <δ } is defined to be 1 for allpairs i, j satisfying | q i − q j | < δ and 0 otherwise. The CM − CM bonds have force constant k b = 6 .
16 kcal/mol and equilibrium length r = 3 .
64 ˚A, and the CM − CT bonds have forceconstant k b = 6 .
16 kcal/mol and equilibrium length r = 3 .
65 ˚A. The force constants forthe CM − CM − CM and CM − CM − CT angles were the same value: k a = 1 .
19 kcal/mol · rad ,equilibrium angles were θ = 173 ◦ and θ = 175 ◦ , respectively.To construct a reference ensemble, we first identified a value for the integration time stepthat was guaranteed not to introduce artifacts. The fastest CG bond oscillation is of thetype CM − CM. The mass of a CM particle is 42.7097 amu, giving a frequency of oscillationbetween two CM particles of ν bond = 12 π (cid:114) . µ fs − = 0 . − with reduced mass µ = 21 . /ν bond ≈ .
44 fs. Thus a time step of 5 fs resolves
OMPARISON OF LANGEVIN INTEGRATORS FOR COARSE-GRAINED POLYMER SIMULATIONS 11
Figure 1.
Graphical representation of the polymer melt inside the unit cell (blue): 128 CG C polymer chains. This rendering was made using VMD after an initial thermalisation run. well the time evolution of the CG bond vibration forces, almost 50 time steps per oscillation.Meanwhile, the Lennard-Jones forces have a characteristic frequency of about ν LJ = τ − = (cid:114) (cid:15)mσ ≈ − , showing that the LJ forces are extremely well-resolved for all choices of time step. Both ν bond and ν LJ suggest that our choice of h = 5 fs is sufficiently small.We used the molecular dynamics engine LAMMPS [29], and implemented the integratorsconsidered here using its Python interface ( fix python/move ) to the underlying data struc-tures. The BAOAB, BBK and BBK ∗ integrators are not available as packages in LAMMPSand were instead implemented using this Python interface. Although a G–JF option is avail-able for the langevin fix command, the version of LAMMPS at the time of writing doesnot implement G–JF in same way as given in [12]. Instead, it uses uniform random variablesto approximate the Gaussian noise, in an effort to increase computational speed. However,such an approximation is only valid for small enough time steps [9]. As we are interested inthe large time step regime, we implemented the original version of G–JF with Gaussian noiseusing the Python wrapper based on equations (20) and (21) in [12].For each choice of γ and h , 100 independent simulations were conducted to reduce the errorassociated to finite length simulations in approximating phase space averages. Each simulationwas performed from an identical spatial configuration for approximately 250 ns. This startingconfiguration was obtained by an initialisation run using LAMMPS’ fix npt command, whichimplements an MTK thermostat/barostat [23], for 100 ns with a temperature of 450 K andpressure set to 1 bar. All simulations used to benchmark the Langevin integrators were runin the NVT ensemble. 4. Simulation Results
The goal of our study is to understand how faithfully the different integrators reproducerelevant statistical averages of the coarse-grained model, particularly in the regime of largetime steps and γ values. To that end, numerical experiments were performed using a range of friction parameters and time steps. We examined the cases of γ ∈ { . , . , . , . } and h ∈ { , , , , , , } with units of fs − and fs, respectively.In the simple case of a Brownian particle in a fluid, γ represents the rate of collision of theBrownian particle with bath particles. So the dimensionless quantity γh gives a measure asto how many collisions occur over the length of the time step h . Similarly here, the number γh determines the strength of the interaction between the system and the heat bath and is afundamental quantity of the dynamics.In applications of Langevin dynamics with atomistic models, friction rate parameters ofthe order of .
01 fs − or less are typically used, as they give a reasonable approximation ofthe experimental diffusion coefficients of small molecules. Indeed in [6], the authors usedBBK in an atomistic water simulation as their benchmark test, and γ was chosen quite small,0 . − . In [17], the BAOAB method was tested and compared to other integratorsbased on its performance on an atomistic alanine dipeptide molecule in water with γ =0 .
001 fs − .However, in CG simulation models, much of the magnitude of the inter-atomic forces shiftsfrom the conservative to the stochastic terms, and a higher friction rate γ may be needed toretain the same diffusivity. Alternatively, high friction rates are also used simply to improvenumerical stability of MD simulations near particular conditions (for example, near phasetransitions). Therefore, we consider here a relatively broad range of γ values, 0 . − to0 . − : given the choice of integration time step h used in the following, the upper end ofthis interval may result in values of γh larger than 1.For smaller γ values, i.e., γ ≤ .
001 fs − , the three integrators become numerically unstablearound h ≈
38 fs. Hence, 35 fs was chosen as the upper bound for our range of time stepvalues. This stability limit for G–JF and BAOAB increases significantly for the largest choiceof γ = 0 . − to slightly more than 50 fs. In contrast, this larger value for γ seemed to nothave as much of an effect on BBK ∗ — simulations still exhibited instability at 40 fs.4.1. Diffusive Behavior.
To characterise diffusion, the mean squared displacement (MSD)of individual molecules was computed as a function of the simulation time. After each timestep in the simulation, we computed the MSD of the center-of-mass for each polymer chainover that time step, using the LAMMPS command compute msd , and then averaged thisresult over all chains and added this to the same calculation from the previous step, i.e., wecomputed: D h ( t k ) := 1 N N (cid:88) i =1 | Q CM i ( t k ) − Q CM i ( t k − ) | + D h ( t k − ) , where t k = kh , 1 ≤ k ≤ N , and Q CM i represents the center of mass position for the i -thpolymer chain ( i = 1 , .., γ and h , yielding5000 independent samples (each block was taken long enough to allow for correlations to dieoff). This was done for each of the three integrators. The MSD was observed to be linear withrespect to time within statistical error. These linear plots were fitted with regression lines andthe means of the slopes of these lines are plotted as a function of the time step h in Figure 2. Asdiscussed earlier, the diffusion coefficient for BAOAB in the case of the one-dimensional singleparticle with zero net external potential does not adhere to the Einstein diffusion relation (seeTable 1). In that same spirit, the diffusion coefficient increases linearly in the upper left plot OMPARISON OF LANGEVIN INTEGRATORS FOR COARSE-GRAINED POLYMER SIMULATIONS 13 time step [fs]12345 s l o p e [ ˚A / n s ] γ = .1 G–JFBAOABBBK* time step [fs]1719212325 s l o p e [ ˚A / n s ] γ = .01 G–JFBAOABBBK* time step [fs]90105120135150 s l o p e [ ˚A / n s ] γ = .001 G–JFBAOABBBK* time step [fs]150175200225250 s l o p e [ ˚A / n s ] γ = .0001 G–JFBAOABBBK*
Figure 2.
Diffusion coefficients (with units of ˚A /ns) were calculated from numerical simulationfor different γ values as a function of the time step h . The plot in the upper left displays the samekind of linear behavior for BAOAB as the method exhibits in the one-dimensional free particle.In the remaining cases, all three methods produce the same calculated slopes within statisticalerror as a function of h . of Fig. 2 as h (and thus γh ) increases. The computed diffusion coefficients for the G–JF andBBK ∗ integrators are relatively unchanged as a function of the time step, again, in line withthe behavior in the simple one-dimensional case. This is a desirable property, as it providesevidence that using larger time steps with G–JF and BBK does not corrupt the system’sdiffusive behavior for any choice of γ . For γ ≤ .
01 fs − , all integrators exhibit statisticallysimilar diffusive behavior, giving confidence that the choice of integrator should not influencediffusion in this regime.We would like to stress the fact that in this study, we observe only classical diffusion,unlike some previous studies. In [27] and [31], polymer melts with CG particles composedof three CH2 monomers, as considered here, were studied. Sub-diffusion was observed forthe quantity D h for simulation times up to and exceeding our simulation time of 250 ns(although the polymer lengths were at least double ours). Therefore we initially consideredthe possibility of anomalous diffusion during our simulation, however no such power law wasobserved. Further examination of error residuals with MSD data and regression lines did notprovide evidence of non-linear relationships between MSD and simulation time.4.2. Configurational Averages.
An important quantity typically used in statistical ther-modynamics is the radial distribution function (RDF). To discriminate the distributions of r [˚A] . . . . . . . g ( r ) Figure 3.
For each γ value and each integrator, the CM—CM radial distribution function for h = 5 fs was calculated. These radial distribution functions were then used as a proxy for thetrue CG radial distribution function. For a given γ value, we denote these reference RDF’s as g γ G–JF ( r ), g γ BAOAB ( r ) and g γ BBK ∗ ( r ). These were visually indistinguishable and so we only show g . ( r ) here. inter-molecular contacts from intra-molecular ones, we restricted the computation of the RDFto pairs of particles from distinct chains. The CM—CM distributions are here examined: be-cause the system is mainly made up of CM particles, this RDF is the most fully sampled.For each friction parameter γ and time step h , the RDF for the intermolecular CM − CMparticle pairing was calculated every 1000 time steps, and these calculated distributions werethen averaged over time. This was done for each of the 100 independent simulations, followedby a final averaging over these 100 simulations in order to reduce sampling error.Lacking analytical expressions for the true RDF, computations conducted with time step h = 5 fs are used as a reference solution. Given how small h = 5 fs is relatively to the timescale of the overall processes, it is reasonable to assume that the true intermolecular CM − CMRDF for the CG system will be well-approximated by the one calculated for h = 5 fs. For each γ value, we calculate the reference RDF’s: g γ G–JF ( r ), g γ BAOAB ( r ) and g γ BBK ∗ ( r ). Then the L relative differences between the reference RDFs and the RDFs obtained for the other choicesof h are computed and plotted as functions of h . These results are displayed in Figure 4.In the large γ regime, BBK ∗ exhibits more and more deviation from its baseline RDF as h is increased, while the other two methods remain unchanged. In fact, although not displayedin the plot, for h = 35 fs, the BBK ∗ error is larger than that of BAOAB and G–JF by an orderof magnitude. If one considered here the original BBK formulation, the resulting error wouldbe slightly lower than BBK* (Fig. 7), but still much higher than BAOAB and G–JF. When γ ≤ .
01 fs − , all three integrators (BAOAB, BBK ∗ and G–JF) perform almost identically.This fact should be useful to the practitioner working in the small γ regime when tryingto strike a balance between diffusivity and adherence to the canonical distribution. Notealso that for BBK*, this is a distinguishing property of the formulation used in this paper:the original BBK formulation has a larger error in this regime (Fig. 7). Some intuition isavailable. If we consider the harmonic potential, the mean square position for BBK is given by k b T ω − (1 − h ω m ) − ; therefore in the simple linear case, the configurational statistics dependonly on the size of the time step, much like what is seen in Fig 7. OMPARISON OF LANGEVIN INTEGRATORS FOR COARSE-GRAINED POLYMER SIMULATIONS 15
10 15 20 25 30 35 time step [fs] R e l a t i v e L e rr o r γ = .1 G–JFBAOABBBK* 10 15 20 25 30 35 time step [fs] R e l a t i v e L e rr o r γ = .01 G–JFBAOABBBK*10 15 20 25 30 35 time step [fs] R e l a t i v e L e rr o r γ = .001 G–JFBAOABBBK* 10 15 20 25 30 35 time step [fs] R e l a t i v e L e rr o r γ = .0001 G–JFBAOABBBK*
Figure 4.
Relative error of radial distribution functions for the three numerical schemes for fourdifferent values of γ . For the largest γ value, we notice minimal change in error as the time stepis increased for G–JF and BAOAB; but a more significant error for BBK ∗ . Although outside therange of the panel (a), for the largest integration step, we recorded an error of approximately0.75% in the RDF, which can make an impact on the quality of results. For smaller choices of γ , we observe behavior consistent with Hamiltonian dynamics: an increasing time step leads toincreased error. The divergence of the BBK ∗ RDF from the other RDFs when γ is largest, led us to questionwhat happens when an intermediate γ value is used. Figure 5 displays these results. A smoothtransition of error occurs between γ = .
01 fs − and γ = . − . This suggests that the BBK ∗ RDF exhibits a γh dependence that is not present in the other RDFs—which is not surprising,given that a similar behavior occurs in the simple 1-d harmonic oscillator case. Again, westress that in most atomistic applications, γ is simply not large enough for the γh dependenceto be noticeable. However, this may no longer be the case for CG dynamics where the larger γ regime becomes more relevant.Another configurational quantity of interest for polymer chains is the bond angle distribu-tion. Many properties of polymer melts, as well as their transition between liquid and solidphases, are determined by the propensity of the polymer chains to align. Failing to accuratelyreproduce the angle distribution can dramatically influence the accuracy of the simulation’sthermodynamic properties.Of the two possible triplets CM—CM—CT and CM—CM—CM, the latter was the bestsampled and so was more amenable to statistical inference. For larger γ , G–JF and BAOABproduce minimal to no variation in error as h is increased. Figure 6 displays the results.
10 15 20 25 30 35 time step [fs] R e l a t i v e L e rr o r Radial distribution function γ = . fs − γ = . fs − γ = . fs − γ = . fs − γ = . fs −
10 15 20 25 30 35 time step [fs]
Bond angle distribution γ = . fs − γ = . fs − γ = . fs − γ = . fs − γ = . fs − Figure 5.
Relative error of the radial and bond angle distribution for the BBK ∗ method as γ varies between 10 − fs − and 10 − fs − , taking on the values 10 − , − . , − . , − . , − fs − . We see a smooth transition as γ becomes smaller, not an abrupt phase transition, furthersuggesting dependence of the steady state distribution on the dimensionless quantity γh . BBK ∗ again exhibits systematic errors for larger γ values. The top row in Fig. 6 indicates alarge relative distortion of the bond angle distribution for BBK ∗ . This evidence may lead oneto use caution in applying BBK ∗ to liquid phase CG systems with larger γ values, especiallysince this behavior for BBK ∗ is also observed with the CM—CM intermolecular RDF. Forsmaller γ , that is as the system gets closer to pure Hamiltonian dynamics, again all threeintegrators perform similarly, indicating no significant preference of integrator for this regime.Table 1 showed that in a harmonic potential, the configurational statistics were dependenton γh for BBK ∗ . So it is not unreasonable to expect that a similar γh dependence for otherconfigurational quantities may occur in more complicated situations when using BBK ∗ , asseen in our simulations.We observe very good agreement of BBK ∗ with G–JF and BAOAB for γ ≤ .
01 fs − inFigures 4 and 6. However, for all choices of γ , the original BBK exhibits a clear trend: therelative error increases as the time step h is increased.5. Conclusions
In this paper we systematically studied three different Langevin integrators on a coarse-grained (CG) polymer melt. For the ideal cases of the Brownian motion and harmonicoscillator, key statistical properties were calculated analytically for each integrator, whichprovided guiding insights into diffusive and statistical behavior of realistic molecular systems.In particular, for pure Brownian motion, both BBK ∗ and G–JF capture the true diffusivebehavior exactly for all choices of γ and h ; but BAOAB is only approximate, with the diffusioncoefficient depending on the dimensionless parameter γh . This carried over to the CG-polymersimulation results. In Section 4.1, the calculated diffusion coefficient as a function of timestep, D h , was found to be statistically independent of the time step h for BBK ∗ and G–JFfor all γ values, whereas BAOAB displayed the same type of linear behavior for γh values of O (1) as in the free particle case. OMPARISON OF LANGEVIN INTEGRATORS FOR COARSE-GRAINED POLYMER SIMULATIONS 17
10 15 20 25 30 35 time step [fs] R e l a t i v e L e rr o r γ = .1 G–JFBAOABBBK* 10 15 20 25 30 35 time step [fs] R e l a t i v e L e rr o r γ = .01 G–JFBAOABBBK*10 15 20 25 30 35 time step [fs] R e l a t i v e L e rr o r γ = .001 G–JFBAOABBBK* 10 15 20 25 30 35 time step [fs] R e l a t i v e L e rr o r γ = .0001 G–JFBAOABBBK*
Figure 6.
Relative error of the bond distribution for the three numerical schemes when γ =0 . − , .
01 fs − , .
001 fs − , . − . The reference is taken to be h = 5 fs. As with theRDF, minimal change in error occurs when the time step is increased for G–JF and BAOAB, butfor BBK ∗ it is much larger. At 35 fs, BBK ∗ experiences an almost 1% error in its bond angledistribution. Convergence to Hamiltonian dynamics is seen in the bottom row. The computational results indicate that G–JF is the only integrator among those consideredhere that describes equally well configurational distributions and diffusive behavior over all γ choices. As expected, BBK ∗ and BBK perform poorly for the largest choice of γ . BAOABsamples equally well the configurational distributions, but exhibits a spurious dependence ofthe diffusivity on the integration time step near the high-friction regime. These conclusionshave implications for CG MD simulations, where large ratios between friction and Hamiltonianforces are more frequently encountered than in atomistic ones. The evidence presented in thispaper should be useful to the practitioner in supporting the use of the G–JF thermostat forCG simulations. 6. Acknowledgments
We would like to thank Richard Berger for technical assistance. Part of this research wasfunded by the US Army Research Laboratory under contract number W911NF-16-2-0189.B. Seibold wishes to acknowledge support by NSF grant DMS–1719640. Calculations werecarried out on Temple University’s HPC resources and thus were supported in part by theNational Science Foundation through major research instrumentation grant number 1625061. γ = .1 BBK*BBK γ = .01
10 15 20 25 30 35 time step [fs] γ = .001
10 15 20 25 30 35 time step [fs] γ = .0001 R e l a t i v e L e rr o r Intermolecular RDF γ = .1 BBK*BBK γ = .01
10 15 20 25 30 35 time step [fs] γ = .001
10 15 20 25 30 35 time step [fs] γ = .0001 R e l a t i v e L e rr o r Bond angle distribution
Figure 7.
A comparison of radial and angle distribution errors for BBK and BBK ∗ . OMPARISON OF LANGEVIN INTEGRATORS FOR COARSE-GRAINED POLYMER SIMULATIONS 19
References [1] LAMMPS source code, 2017. http://lammps.sandia.gov/ .[2] Sequencer.c, NAMD source code, 2018. .[3] S. Anmarkrud, K. Debrabant, and A. Kværnø. General order conditions for stochastic partitioned Runge–Kutta methods.
BIT Numerical Mathematics , 58(2):257–280, 2018.[4] E. Arad, O. Farago, and N. Grønbech-Jensen. The G-JF thermostat for accurate configurational samplingin soft-matter simulations.
Isr. J. Chem. , 56(8):629–635, 2016.[5] G. Brannigan, L. Lin, and F. Brown. Implicit solvent simulation models for biomembranes.
Eur. Biophys.J. , 35(2):104–124, 2006.[6] A. Br¨unger, C. Brooks III, and M. Karplus. Stochastic boundary conditions for molecular dynamicssimulations of ST2 water.
Chem. Phys. Lett. , 105(5):495–500, 1984.[7] K. Burrage and G. Lythe. Accurate stationary densities with partitioned numerical methods for stochasticdifferential equations.
SIAM J. Numer. Anal. , 47(3):1601–1618, 2009.[8] I. Cooke, K. Kremer, and M. Deserno. Tunable generic model for fluid bilayer membranes.
Phys. Rev. E ,72(1):011506, 2005.[9] B. D¨unweg and W. Paul. Brownian dynamics simulations without Gaussian random numbers.
Int. J. Mod.Phys. C , 2(03):817–827, 1991.[10] H. Flyvbjerg and H. G. Petersen. Error estimates on averages of correlated data.
J. of Chem. Phys. ,91(1):461–466, 1989.[11] R. Goetz and R. Lipowsky. Computer simulations of bilayer membranes: self-assembly and interfacialtension.
J. Chem. Phys. , 108(17):7397–7409, 1998.[12] N. Grønbech-Jensen and O. Farago. A simple and effective Verlet-type algorithm for simulating langevindynamics.
Mol. Phys. , 111(8):983–991, 2013.[13] N. Grønbech-Jensen, N. Hayre, and O. Farago. Application of the G-JF discrete-time thermostat for fastand accurate molecular simulations.
Comput. Phys. Commun. , 185(2):524–527, 2014.[14] K. Hadley and C. McCabe. A structurally relevant coarse-grained model for cholesterol.
Biophys. J. ,99(9):2896–2905, 2010.[15] J.A. Izaguirre, D.P. Catarello, J.M. Wozniak, and R. Skeel. Langevin stabilization of molecular dynamics.
J. Chem. Phys. , 114(5):2090–2098, 2001.[16] S. Izvekov and G. Voth. A multiscale coarse-graining method for biomolecular systems.
J. Phys. Chem.B , 109(7):2469–2473, 2005.[17] B. Leimkuhler and C. Matthews. Robust and efficient configurational molecular sampling via Langevindynamics.
J. Chem. Phys. , 138(17):05B601 1, 2013.[18] B. Leimkuhler and C. Matthews. Robust and efficient configurational molecular sampling via Langevindynamics.
J. Chem. Phys. , 138(17):05B601 1, 2013.[19] B. Leimkuhler and C. Matthews.
Molecular Dynamics . Springer, 2015.[20] Hans-J¨org Limbach, A. Arnold, B. Mann, and C. Holm. ESPResSo — an extensible simulation packagefor research on soft matter systems.
Comput. Phys. Commun. , 174(9):704–727, 2006.[21] S. Marrink, H. Risselada, S. Yefimov, D. Tieleman, and A. De Vries. The MARTINI force field: coarsegrained model for biomolecular simulations.
J. Phys. Chem. B , 111(27):7812–7824, 2007.[22] G. J. Martyna, M. Klein, and M. Tuckerman. Nos´e–hoover chains: The canonical ensemble via continuousdynamics.
J. Chem. Phys. , 97(4):2635–2643, 1992.[23] G. L. Martyna, D. Tobias, and M. Klein. Constant pressure molecular dynamics algorithms.
J. Chem.Phys. , 101(5):4177–4189, 1994.[24] T. Murtola, M. Karttunen, and I. Vattulainen. Systematic coarse graining from structure using internalstates: Application to phospholipid/cholesterol bilayer.
J. Chem. Phys. , 131(5):08B601, 2009.[25] M. Orsi, D. Haubertin, W. Sanderson, and J. Essex. A quantitative coarse-grain model for lipid bilayers.
J. Phys. Chem. B , 112(3):802–815, 2008.[26] R. Pastor, B. Brooks, and A. Szabo. An analysis of the accuracy of Langevin and molecular dynamicsalgorithms.
Mol. Phys. , 65(6):1409–1419, 1988.[27] B. Peters, K. M. Salerno, A. Agrawal, D. Perahia, and G. Grest. Coarse-grained modeling of polyethylenemelts: effect on dynamics.
J. Chem. Theory Comput. , 13(6):2890–2896, 2017.[28] I. Pivkin, B. Caswell, and G. Karniadakis.
Reviews in Computational Chemistry, Chapter 2 , volume 27.John Wiley & Sons, 2011.[29] S. Plimpton. Fast parallel algorithms for short-range molecular dynamics.
J. Comput. Phys. , 117(1):1–19,1995. [30] Jean-Paul Ryckaert, G. Ciccotti, and H. Berendsen. Numerical integration of the cartesian equations ofmotion of a system with constraints: molecular dynamics of n-alkanes.
J. Comput. Phys. , 23(3):327–341,1977.[31] K. M. Salerno, A. Agrawal, D. Perahia, and G. S. Grest. Resolving dynamic properties of polymers throughcoarse-grained computational studies.
Phys. Rev. Lett. , 116(5):058302, 2016.[32] T. Schneider and E. Stoll. Molecular-dynamics study of a three-dimensional one-component model fordistortive phase transitions.
Phys. Rev. B , 17(3):1302, 1978.[33] X. Shang, M. Kr¨oger, and B. Leimkuhler. Assessing numerical methods for molecular and particle simu-lation.
Soft Matter , 13(45):8565–8578, 2017.[34] J. Shelley, M. Shelley, R. Reeder, S. Bandyopadhyay, and M. Klein. A coarse grain model for phospholipidsimulations.
J. Phys. Chem. B , 105(19):4464–4470, 2001.[35] W. Shinoda, R. DeVane, and M. Klein. Multi-property fitting and parameterization of a coarse grainedmodel for aqueous surfactants.
Mol. Simulat. , 33(1-2):27–36, 2007.[36] M. Tuckerman.
Statistical Mechanics and Molecular Simulations . Oxford University Press, 2008.[37] M. Tuckerman, B. Berne, and G. Martyna. Reversible multiple time scale molecular dynamics.
J. Chem.Phys. , 97(3):1990–2001, 1992.[38] W. Wang and R. Skeel. Analysis of a few numerical integration methods for the Langevin equation.
Mol.Phys. , 101(14):2149–2156, 2003.(Joshua Finkelstein)
Department of Mathematics, Temple University, 1805 North Broad Street,Philadelphia, PA 19122
E-mail address : [email protected] (Giacomo Fiorin) Institute for Computational Molecular Science, Temple University, 1805North Broad Street, Philadelphia, PA 19122 and National Heart, Lung and BloodInstitute, Bethesda, MD
E-mail address : [email protected] (Benjamin Seibold) Department of Mathematics, Temple University, 1805 North Broad Street,Philadelphia, PA 19122
E-mail address : [email protected] URL ::