Group-based evaluation of temperature and pressure for molecular dynamics simulation with a large time step

In molecular dynamics (MD) simulations, under the isothermal and isobaric condition, one evaluates instantaneous temperature and pressure (T/P) and updates the thermostat and barostat by comparing the evaluated T/P values with the target ones.¹ The instantaneous T/P values should be accurately evaluated in MD simulations to obtain reliable results. If instantaneous temperature is overestimated, the simulation system will be overcooled than that required by the thermostat. Similarly, the overestimation of pressure can increase the system size than expected in barostat updates. The underestimations of temperature and pressure cause the opposite effects. These inaccurate T/P evaluations can introduce unexpected phase transitions of the simulation system. For instance, a solvated lipid bilayer in the liquid phase is changed to that in the gel phase around the phase transition point if temperature is largely overestimated. Inaccurate temperature evaluation may break the equipartition of temperature, which introduces the difference in temperatures between the solute and the solvent. This unexpected hot-solvent and cold-solute (or vice versa) phenomenon is widely known as the “hot-solvent/cold-solute” problem in MD simulations.² This unexpected result from the inaccurate T/P evaluations may reduce the validation of MD simulations.³ 

Instantaneous temperature should be evaluated as a product of momentum and time derivative of position based on Tolman’s equipartition theorem.^4–6 Pressure should be evaluated to satisfy the virial theorem when there is no external force.⁷ Unfortunately, commonly used temperature evaluation based on kinetic energy and pressure evaluations based on kinetic energy and virial are accurate up to the first order of a MD integration time step (δt). They cannot be used in isothermal–isobaric MD simulations with a large time step. By combining full- and half-time step kinetic energies, we derived instantaneous T/P that are accurate up to the third order of δt in our previous studies.^5,7 The evaluations of instantaneous T/P work well in isothermal–isobaric MD simulations with a short time step (δt = 1 fs or 2 fs) as well as a large time step (δt = 5 fs) in the velocity Verlet integrator.

In several MD software, the group approach is applied to evaluate instantaneous pressure.^8–10 Here, all the atoms in a simulation system are partitioned into disjoint groups, each of which consists of one heavy atom and bonded hydrogen atoms. The groups instead of atoms are used to evaluate instantaneous pressure. The advantage of this approach over the conventional one (hereafter, we refer to it as the atomic approach) is that better accuracy can be obtained in pressure evaluation because it does not include high-frequency vibrational motions within a group. To our knowledge, there is no implementation of the group-based evaluation of temperature in well-known MD software.

In this paper, we aim at combining the time integrator that includes the instantaneous T/P evaluations up to the third order of δt^5,7 with the group-based T/P evaluations (group T/P) for further improving the efficiency of the MD integrator. Based on our T/P evaluations, the group-based approach can further increase the efficiencies for several reasons. First of all, in the updates of the thermostat and barostat, self-consistent iterations are not required even if SHAKE¹¹ or RATTLE¹² constraints are applied to the bond lengths. Second, a multiple time step (MTS) integration, such as reversible reference system propagation algorithm (r-RESPA),¹³ can be applied by avoiding the ambiguity of T/P evaluations. Our approach is expected to outperform the other integrators in terms of performance and accuracy.

This paper is organized as follows: In Sec. II, we describe T/P evaluations, which are accurate up to the third order of δt, and the isothermal–isobaric MD integrator using the T/P evaluations. Then, we combine the integrator with the group-based T/P evaluations. In Sec. III, we show the extensive MD simulations of two systems: a pure water box and DPPC (1,2-dipalmitoyl-sn-phosphatidylcholine) lipid bilayer systems (80 DPPC and 160 DPPC). A pure water box is chosen as a test system for the isotropic pressure coupling. 80 DPPC and 160 DPPC are selected because the bilayer structures are sensitive to the system properties around the transition temperature. They are also suitable to understand structural properties related to semi-isotropic pressure couplings. Structural properties of 80 DPPC and 160 DPPC are examined in the MD simulations under different conditions {different T/P definitions, atomic T/P vs group T/P, different time steps (δt = 2.0 fs or 5.0 fs), different integrators [velocity Verlet or multiple time step (MTS)], different thermostat algorithms [Berendsen,¹⁴ Bussi,¹⁵ and Nose–Hoover chain (NHC)¹⁶], hydrogen mass repartitioning (HMR) schemes,^17,18 and so on}. Based on the theoretical derivations and numerical simulation results, we conclude in Sec. IV that the proposed T/P evaluations are robust for different simulation conditions and give reliable MD simulation trajectories in the isothermal–isobaric condition.

Here, we derive the group-based T/P evaluations in the framework of previously proposed T/P evaluations in the MD integrator, which are accurate up to the third order of δt. In Sec. II A, we derive the accurate T/P evaluations suggested in our previous works.^5,7 In this derivation, we newly discuss the order of accuracy in the T/P evaluations. In Sec. II B, MD integration with isothermal–isobaric conditions is described. From Sec. II C, we start theoretical derivations of the group T/P approach and explain how to introduce it to our previous T/P evaluation framework.

One can derive the accurate T/P evaluations in MD based on the Boltzmann distributions.⁶ For an N particle system with momentum p ≡ {p₁, p₂, …, p_N} and position vector r ≡ {r₁, r₂, …, r_N}, the canonical ensemble partition function Z(N,V,T) at temperature T is given by

where $β=1kBT$⁠, k_B is the Boltzmann constant, and H is the Hamiltonian. By replacing $∂H∂pie−βH$ with $−1β∂e−βH∂pi=−kBT∂e−βH∂pi$ and applying the integration by parts, $pi⋅ṙi$ and $ri⋅ṗi$ result in

where d is the number of dimensions of each particle. For a Hamiltonian system, average kinetic energy or virial are related to temperature with the following relationship:

where N_f is the number of degrees of freedom. In MD simulations, this relationship cannot be directly applied because of truncation error with a finite time step.

To understand the relationship between temperature and kinetic energy or viral in MD, we first consider the velocity Verlet integration scheme using a time step δt. In this integration, one cycle of position and momentum vector updates is

where $p̃it+12δt$ is the intermediate momentum of the ith particle propagated from p_i(t) by $12δt$⁠. We do not write it as $pit+12δt$ because $p̃it+12δt$ and Taylor expansion of $pit+12δt$ are different from each other. Because there are two momenta in this integration, we can define kinetic energy at t, K(t), in two ways. One is the kinetic energy using p_i(t),

We name it the full-time step kinetic energy because momentum at the full time step is used for the evaluation of kinetic energy. The other is obtained from $p̃it+12δt$⁠,

We name it the half-time step kinetic energy. At any time t, K_half is always greater than K_full with the following relationship:

To derive the accurate temperature estimation based on Eq. (2), the first step is to understand the relationship between the time derivatives of momentum and force. The accurate time derivative of momentum can be obtained from the Taylor expansion of momentum, and force is directly involved in momentum updates in MD integration. To connect these, let us see the difference between $p̃it+12δt$ and $p̃it−12δt$ into two ways, i.e., Taylor expansion of $p̃it+12δt$ and $p̃it−12δt$ in terms of $p̃it$⁠, and MD integration leads the following relationship between $p̃̇it$ and F_i(t):

Similarly, the Taylor expansion and MD integration expression for the sum of $p̃it+12δt$ and $p̃it−12δt$ relates p_i(t) and $p̃i(t)$ as follows:

By applying the time derivatives in Eq. (9) and canceling out $p̃̇it$⁠, one can obtain the relationship between $ṗi(t)$ and F_i(t) as follows:

In other words, a time derivative of momentum can be expressed in terms of force as the accuracy of the first order of δt.

Optimal temperature evaluation is obtained from the relationship between time derivatives of position vector and momentum, and that between time derivatives of momentum and force. For the understanding of the relationship between time derivatives of position vector and momentum, we express the difference between r_i(t + δt) and r_i(t − δt) using Taylor expansion and MD integration expressions,

By canceling out r_i(t + δt) − r_i(t − δt) and applying the time average after multiplying p_i(t), temperature is expressed as

Therefore, temperature evaluated from K_full is accurate up to the δt term and does not guarantee accurate temperature for a large time step. $ri(3)t$ is equivalent to $pi(2)(t)m−16ri(5)tδt2+Oδt4$ from Eq. (11), so the δt² term in Eq. (12) is

By making use of $pi(t)⋅pi(2)t=−ṗi(t)2$ as well as the relationship between $ṗi(t)$ and F_i(t) in Eq. (10), N_fk_BT is expressed in terms of momentum and force,

From the relationship shown in Eq. (7), one can finally derive the optimal temperature using both ⟨K_full⟩ and ⟨K_half⟩,

The relationship between kinetic energy and virial is obtained from the temperature evaluation in Eq. (3) and the relationship between time derivative of momentum and force shown in Eq. (10),

The δt² term in Eq. (16) is again expressed as

It can be replaced by K_half − K_full, and finally, one can obtain

without external forces. This means that the virial theorem is satisfied up to the third order of δt only when evaluating the kinetic energy as K_half. Pressure evaluation should be based on the relationship between virial and kinetic energy in Eq. (18) and can be expressed as

for the isotropic case using half-time step kinetic energy.

Under the NPT (isobaric-isothermal) condition, the equation of motion for isotropic volume fluctuations is written by introducing additional dynamical variables related to the barostat and thermostat, p_ε and p_η. These variables scale position vectors and momenta,¹⁹ 

where W, V, P_ext, K, p_ε, and p_η are barostat mass, volume, target pressure, instantaneous kinetic energy of the system, barostat momentum, and thermostat momentum for particles/barostat, respectively. Temperature is usually obtained from the sum of kinetic energies of particles and a barostat. There are also cases where thermostat momenta for particles and barostat are separated¹⁹ or thermostat momentum for the barostat is neglected.^8,p_η is mainly decided from the difference between instantaneous and target temperatures, and the detailed ways of update are varied according to thermostat schemes. For example, in the case of Nose–Hoover chain (NHC) thermostat,¹⁶ there are additional equations for chains. In the Berendsen thermostat, the rescaling factor is obtained just by the difference between the target and instantaneous temperature divided by the damping constant.¹⁴ Bussi’s stochastic velocity rescaling (SVR) thermostat is similar to the Berendsen thermostat, but it includes the stochastic term to generate statistically meaningful temperature distributions.²⁰ From the equation of motions, the phase space is propagated by

where

where iL_R and iL_P can further be decomposed into

and

respectively. In the case of anisotropic or semi-isotropic case, p_ε is no longer scalar, but a matrix p_ε. With the rectangular system size, all off-diagonal terms of p_ε become zero. In addition, further restrictions based on the rectangular system size can also be assigned. For example, p_ε,xx = p_ε,yy and p_ε,xx = p_ε,yy = 0 conditions are assigned for semi-isotropic and NPAT conditions, respectively.

Considering the accurate estimation of pressure and temperature described in the previous sections, instantaneous temperature and pressure at time t should be evaluated as

Let us name the instantaneous temperature and pressure defined in Eq. (25) the optimal temperature and half-step pressure. For both temperature and pressure evaluations at time t, $p̃it+12δt$ is required because of K_half(t) evaluation. If there is no holonomic bond constraint, K_half(t) can be evaluated without propagating momentum from p_i(t) to $p̃it+12δt$ but just using the following symmetry:

K_half(t) is then evaluated as

If there are constraints, however, propagation from p_i(t) to $p̃it+12δt$ is necessary to evaluate K_half(t) because $pi(t)−p̃it−12δt=p̃it+12δt−pi(t)$ is not generally true. In this case, the temperature and pressure evaluations, and followed propagations of positions and momenta should be done in an iterative way. First, we make an initial guess of K_half(t) and estimate instantaneous temperature T_opt(t) and pressure P_half(t). Based on the estimated T_opt(t) and P_half(t), thermostat and barostat related variables, p_ε and p_η, are updated and followed by position vector and momenta updates. Updated position vectors and momenta are again used for K_half(t). This is done until the temperature and pressure values are converged.

The overall integration scheme can be simplified by applying group temperature and pressure (group T/P). The group pressure approach or even molecular pressure approaches are used in several MD programs,^8,9 while the same approach was not applied to the evaluation of temperature. Let us assume that the particles are partitioned into N_G disjoint groups, each of which consists of one heavy atom and bonded hydrogen atoms. If a heavy atom is not bonded to any hydrogen atom, the heavy atom itself consists of one group. Let us define mass, position, momentum, and force of each group as follows:

Let us also define the relative positions and momenta as follows:

In terms of group coordinate R_g and momentum P_g, temperature evaluated as $NfkBT=∑ipit⋅ṙit$ is expressed as

Let us assume that N₁ and N₂ are the numbers of degrees of freedoms related to the overall translation of groups and the relative motion in each group, respectively. The sum of N₁ and N₂ is equal to N_f, and temperatures evaluated by $∑gPg⋅Ṙg$ and $∑g∑i∈gpgi⋅ṙgi$ are equal to each other from the generalized equipartition theorem,

Therefore, we can evaluate temperature just considering the group approach,

In velocity Verlet integration, $∑gPg⋅Ṙg$ is equal to

where

Accordingly, we update the thermostat by evaluating instantaneous temperature as

In the same way, the instantaneous group pressure is evaluated as

In each group, constraint does not change the center of mass motion because the sum of constraint forces is zero, so constraint does not change the group temperature and pressure values. With these, half-time step kinetic energy is obtained without propagating momentum from p_i(t) to $p̃it+12Δt$ but just using

and the iteration procedure is not needed for the evaluation of K_half,g. It should be noted that the group approach temperature evaluation can be applied only when the equipartition theorem is satisfied, i.e., temperature is evaluated in an accurate way.

In a thermostat procedure, we update p_η from the difference between instantaneous and target temperatures, which is followed by the momentum scaling of particles. With group temperature, we scale only the center of mass momenta of each group,

where g is the group index that includes the ith particle and $vη=pηQ$⁠. Without group pressure, for given values of $vε=pεW$ and $α=1+3Nf$⁠, the update of momenta and position vectors in the first half-time step integration, corresponding to $eiLPΔt/2$ and $eiLRΔt$⁠, is done by using hyperbolic sine functions of αv_εδt/4 and v_εδt/2.¹ With the group pressure approach, we only scale the center of mass position vectors and momenta of each group as suggested by Lippert et al.⁸ The update of momenta and position vectors is then expressed as

With the multiple time step (MTS) integration, some propagators are applied less frequently with a larger time step. One of the most famous ways in MTS integration is to split forces into two: those related to the slow and fast motions: f_i = f_i,fast + f_i,slow.

If fast- and slow-motion forces are evaluated with the time steps of δt and Δt, slow-motion forces are evaluated at every nth step, Δt = nδt, and a time propagator at the increment of Δt under the NVE condition is¹³ 

where

Evaluations of T/P values required for thermostat and barostat updates are not straightforward because of K_half(t) evaluation at time t, which requires propagation of momenta and position vectors. Evaluation of K_half(t) requires propagation of momentum with a half time step. With MTS integration, there are two time steps, δt and Δt, and momentum at the half time step cannot be obtained. We suggested one solution in our previous work, but it requires self-consistent iterations in thermostat and barostat updates, making the integration process inefficient. With the group approach, T/P are evaluated without any propagation of position vectors and momenta because a bond constraint does not affect instantaneous T/P values. As a result, we can avoid such inefficient self-consistent iteration processes. Without MTS integration, $Pg(t)−P̃gt−12δt$ is replaced by

and F_g can be the sum of slow- and fast-motion forces: F_g = F_long,g + F_short,g. Therefore, the half-time step kinetic energy expressed in the velocity Verlet equation can be reformulated as

and this is used for the evaluations of temperature and pressure in MTS integration. In MTS integration, we can also consider the calculations of thermostat or barostat less frequently with a much larger time step. At every time, we can evaluate temperature and pressure as follows:

where Δt_baro, Δt_therm, and Δt are related by Δt_baro = bΔt_therm and Δt_therm = aδt with integers a and b greater than or equal to 1. When splitting forces with slow and fast motions, the one cycle propagator is again represented as

A pure water box consists of 9250 water molecules in a cubic box. The interactions are described with the CHARMM TIP3P model.²¹ Two DPPC lipid bilayer systems are also prepared, consisting of 80 and 160 DPPC (1,2-dipalmitoyl-sn-phosphatidylcholine) molecules. They are solvated with 2701 and 5521 CHARMM TIP3P water molecules, respectively. The initial models of the DPPC lipid bilayers were constructed using the CHARMM-GUI membrane builder.²² CHARMM36 force fields are used for lipid molecules.^23,24 In all tests, water molecules were treated as rigid using SETTLE constraints.

There are mainly three test conditions in this study: δt = 2 fs with velocity Verlet integration (VV_2 fs), δt = 5 fs with velocity Verlet integration (VV_5 fs), and MTS integration with δt = 3.5 fs and Δt = 7.0 fs (MTS_3.5 fs). For VV_2 fs and VV_5 fs, we examined both atomic and group temperature/pressure (T/P) evaluations. In MTS integration, thermostat/barostat is applied at every 21 fs, i.e., Δt_tb = 21 fs, and only group T/P is tested. Slow motions are originated from the electrostatic reciprocal-space interactions in the particle mesh Ewald approximation.²⁵ To test temperature and pressure evaluations, we define instantaneous temperature and pressure with the full time step,

as well as our suggestions, T_opt(t) and P_half(t).

In this test, we observed density of water molecules by changing the time step and T/P evaluations under isothermal–isobaric conditions. To understand the dependency of density on T/P evaluations, we adopted four T/P evaluations: T_opt/P_half, T_opt/P_full, T_full/P_half, and T_full/P_full. In all cases, stochastic velocity rescaling (SVR) thermostats¹⁵ with the MTK barostat in isotropic pressure coupling¹⁹ were assigned. MD simulation time in each case was 20 ns, and the first 2 ns simulations were not included in the density evaluation. According to Table I, T_full/P_half and T_opt/P_full have non-negligible errors with underestimation and overestimation of densities, respectively, even with δt = 2 fs when atomic T/P is used. By increasing the time step up to δt = 5 fs, the amount of underestimation and overestimation increases significantly, and all T/P evaluations except T_full/P_half have non-negligible errors. Group T/P reduces the errors for all time steps, and all T/P evaluations are acceptable when δt = 2 fs is used. T_opt/P_half has consistent results irrespective of group T/P or time steps. T_full/P_full also generates consistent density results even with δt = 5 fs when group T/P is applied, but this is a compensation of temperature underestimation from T_full and pressure overestimation from P_full. It can be understood from the MD simulations of heterogeneous system with anisotropic pressure evaluation, as shown in Sec. III. D 2.

Here, we prepared two working conditions. The first is 20 ns NVT simulations to compare temperatures of the solvent and solute and those of atomic and group evaluations. The second is 300 ns NPT simulations to examine numerical stability by comparing various properties. The target temperature and pressure are 323.15 K (for both NVT and NPT) and 1 bar (for NPT). In the NVT simulations, the first 2 ns simulation data were considered as equilibrium processes and the last 18 ns simulation data were analyzed. In the case of NPT, we analyzed the data from the last 250 ns trajectories by assuming that the first 50 ns simulations are involved in equilibrium. The initial box sizes of 80 DPPC and 160 DPPC are (52.141 Å, 52.141 Å, and 66.444 Å) and (68.789 Å, 68.789 Å, and 76.441 Å). In these simulations, we assigned the semi-isotropic condition to observe the effect of pressure evaluation in each dimension. In NPT MD simulations of DPPC systems, we observed the statistics of area per lipid (APL), volume, order parameter of the first chain, and lipid thickness. The APL was calculated by dividing the system area in the xy dimension by the number of lipids in each leaflet. Lipid thickness was calculated from the difference of z coordinates of phosphorous atoms in the upper and lower leaflets. Average and standard errors are from the block averages of 25 ns simulation length. The order parameter for each chain was calculated by

where $θ$ is the angle between the CH-bond vector and the bilayer normal (z-axis). APL, volume, lipid thickness, and order parameters are the structural properties that can be compared to the experimental values and are sensitive to the lipid state under the working conditions. If NPT MD simulations were not accurate, these quantities would be affected significantly, suggesting the deformation of the lipid bilayer. In all cases, we assigned SHAKE/RATTLE/SETTLE constraints.^11,12,26

From the NVT simulations, we first observed temperatures of the solvent and solute in the 80 DPPC system. The results are written in Table II. In all thermostats, irrespective of temperature style, temperature evaluation using T_full results in the hot-solvent/cold-solute phenomenon, about 6 K and 3.5 K differences using atomic and group temperatures, respectively. While group temperature reduces the temperature difference, the difference is still not negligible.

To understand the effect of group T/P in the accuracy point of view, we first investigate physical and structural properties of 80 DPPC with δt = 2 fs with three combinations of temperature and pressure. From the comparison between T_full/P_full and T_full/P_half in Fig. 1 and Table III, we observed that the group pressure approach reduces the errors from inaccurate pressure evaluations, and there is no significant difference in structural properties from different pressure evaluations if δt = 2 fs. By comparing T_full/P_half and T_opt/P_half, one can understand that the group temperature approach reduces the errors of temperature evaluations, too. When evaluating temperature and pressure in conventional ways, δt² term errors are mainly due to the vibrational motions of the chemical groups involving hydrogen atoms. In the group T/P approach, we can neglect the vibrational motions in each hydrogen group and reduce errors. From the observations, we suggest that any choice of T/P evaluations is okay for δt = 2 fs. Group T/P also has the effect to increase the integration speed by skipping iterations in thermostat and barostat updates. With atomic T/P, the MD integration time including coordinate, position, barostat, and thermostat updates for 300 ns was 27 040 s on a single node of Intel central processing unit (CPU) (E5-2690), and it was reduced to 24 240 s with group T/P when we use T_opt/P_half evaluations.

To understand the effect of T/P evaluations in MD integration with a large time step, we performed MD simulations of 80 DPPC with δt = 2 fs and δt = 5 fs and compared the structural and physical properties of a DPPC lipid bilayer, which are described in Fig. 2 and Table III. If T_opt/P_half evaluations are used, there are no significant differences of structural properties, namely, the order parameters of carbon atoms in the first chain of a DPPC molecule, the membrane volume, the area per lipid (APL), and the membrane thickness, between δt = 2 fs and δt = 5 fs regardless of the atomic and group approaches. In contrast, if T_full/P_half or T_full/P_full evaluations are used, non-negligible errors in the structural properties are observed, in particular, from MD simulations with δt = 5 fs. The errors in MD simulations using atomic T/P are much greater than those using group T/P evaluations, suggesting the advantages of the group approach. One point is that different errors can be compensated to each other by chance, and we should avoid deciding only from one structural parameter. For example, with T_full/P_full evaluations, the atomic approach seems to be better than the group approach only considering the distribution of volume [Fig. 2(b)]. In fact, the atomic approach with T_full/P_full evaluations overestimates the APL and underestimates thickness, and the accurate volume estimation is obtained just by compensating the over- and under-estimated values.³ 

To understand the effect of different thermostat algorithms using the same T/P evaluations, we compared three thermostats: stochastic velocity rescaling (SVR),¹⁵ Berendsen,¹⁴ and Nose–Hoover chain¹⁶ thermostats. We found that there are no significant changes in physical and structural properties (Fig. 3 and Table IV) when we used SVR and Nose–Hoover thermostats. In contrast, the Berendsen thermostat has a different temperature distribution from those using SVR and Nose–Hoover ones. Only considering the average structural properties, accuracy of T/P evaluations could give a more high impact than thermostat styles.