Recently, we proposed novel temperature and pressure evaluations in molecular dynamics (MD) simulations to preserve the accuracy up to the third order of a time step, δt [J. Jung, C. Kobayashi, and Y. Sugita, J. Chem. Theory Comput. 15, 84–94 (2019); J. Jung, C. Kobayashi, and Y. Sugita, J. Chem. Phys. 148, 164109 (2018)]. These approaches allow us to extend δt of MD simulations under an isothermal–isobaric condition up to 5 fs with a velocity Verlet integrator. Here, we further improve the isothermal–isobaric MD integration by introducing the group-based evaluations of system temperature and pressure to our previous approach. The group-based scheme increases the accuracy even using inaccurate temperature and pressure evaluations by neglecting the high-frequency vibrational motions of hydrogen atoms. It also improves the overall performance by avoiding iterations in thermostat and barostat updates and by allowing a multiple time step integration such as r-RESPA (reversible reference system propagation algorithm) with our proposed high-precision evaluations of temperature and pressure. Now, the improved integration scheme conserves physical properties of lipid bilayer systems up to δt = 5 fs with velocity Verlet as well as δt = 3.5 fs for fast motions in r-RESPA, respectively.

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.1 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.2 This unexpected result from the inaccurate T/P evaluations may reduce the validation of MD simulations.3 

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.7 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 δt5,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 SHAKE11 or RATTLE12 constraints are applied to the bond lengths. Second, a multiple time step (MTS) integration, such as reversible reference system propagation algorithm (r-RESPA),13 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,14 Bussi,15 and Nose–Hoover chain (NHC)16], 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.6 For an N particle system with momentum p ≡ {p1, p2, …, pN} and position vector r ≡ {r1, r2, …, rN}, the canonical ensemble partition function Z(N,V,T) at temperature T is given by

Z(N,V,T)=1Cdp1dpNdr1drNeβH,
(1)

where β=1kBT, kB is the Boltzmann constant, and H is the Hamiltonian. By replacing HpieβH with 1βeβHpi=kBTeβHpi and applying the integration by parts, piṙi and riṗi result in

piṙi=riṗi=dkBT,
(2)

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:

K=i=1Npi22mi=12NfkBT,Vir=i=1NriFi=NfkBT,
(3)

where Nf 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

p̃it+12δt=pit+12Fitδt,rit+δt=rit+p̃it+12δtδtmi,pi(t+δt)=p̃it+12δt+12Fi(t+δt)δt,
(4)

where p̃it+12δt is the intermediate momentum of the ith particle propagated from pi(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 pi(t),

Kfull(t)=i=1Npi(t)22mi.
(5)

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,

Khalf(t)=12i=1Np̃it12δt22mi+p̃it+12δt22mi.
(6)

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

Khalf(t)=Kfull(t)+δt28i=1NFi(t)2mi.
(7)

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̃it12δt into two ways, i.e., Taylor expansion of p̃it+12δt and p̃it12δt in terms of p̃it, and MD integration leads the following relationship between p̃̇it and Fi(t):

Fi(t)=p̃̇i(t)+124p̃i(3)tδt2+Oδt4.
(8)

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

pi(t)=p̃i(t)+18p̃i(2)tδt2+Oδt4.
(9)

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

ṗi(t)=Fi(t)+112p̃i(3)tδt2+Oδt4=Fi(t)+Oδt2.
(10)

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 ri(t + δt) and ri(tδt) using Taylor expansion and MD integration expressions,

ri(t+δt)ri(tδt)=2ṙi(t)δt+13ri(3)tδt3+Oδt5,ri(t+δt)ri(tδt)=2pi(t)miδt.
(11)

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

NfkBT=i=1Npi(t)2mδt26i=1Npi(t)ri(3)t+Oδt4.
(12)

Therefore, temperature evaluated from Kfull 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)m16ri(5)tδt2+Oδt4 from Eq. (11), so the δt2 term in Eq. (12) is

δt26i=1Npi(t)pi(2)tm+Oδt4.
(13)

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

NfkBT=i=1Npi(t)2m+δt26i=1NFi(t)2m+Oδt4.
(14)

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

NfkBT=23Kfull+43Khalf+Oδt4.
(15)

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

NfkBT=i=1Nri(t)Fi(t)δt212i=1Nri(t)pi(3)t+Oδt4.
(16)

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

δt212pi(t)pi(2)tmiOδt2=δt212ṗi(t)2miOδt2=δt212Fi(t)2mi+Oδt4.
(17)

It can be replaced by KhalfKfull, and finally, one can obtain

i=1Nri(t)Fi(t)+2Khalf+Oδt4=0,
(18)

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 Khalf. Pressure evaluation should be based on the relationship between virial and kinetic energy in Eq. (18) and can be expressed as

P=13Viri(t)Fi(t)+2Khalf(t)
(19)

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

ṙi=pimi+pεWri,ṗi=FipεW1+3Nfpi+pηQpi,V̇=3VpεW,ṗε=VPtPext+2KNf+pηQpε=Fε+vηpε,
(20)

where W, V, Pext, 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 separated19 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,16 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.14 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.20 From the equation of motions, the phase space is propagated by

eiLΔt=eiLB2δt/2eiLTδt/2eiLPδt/2eiLB1δteiLRδteiLPδt/2eiLTδt/2eiLB2δt/2,
(21)

where

iLR=ipimi+pεWriri,iLP=ifipεW1+3Nfpipi,iLB1=3VpεWV,iLB2=Fεpε,iLT=pηQipipi+pηQpεpε,
(22)

where iLR and iLP can further be decomposed into

iLR1=ipimiri,iLR2=ipεWriri
(23)

and

iLP1=ifipi,iLP2=pεW1+3Nfipipi,
(24)

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

Topt(t)=1NfkB23Kfull(t)+43Khalf(t),Phalf(t)=13Viri(t)Fi(t)+2Khalf(t).
(25)

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 Khalf(t) evaluation. If there is no holonomic bond constraint, Khalf(t) can be evaluated without propagating momentum from pi(t) to p̃it+12δt but just using the following symmetry:

pi(t)p̃it12δt=p̃it+12δtpi(t).
(26)

Khalf(t) is then evaluated as

Khalf(t)=i=1N12mipi(t)2+pi(t)p̃it12δt2.
(27)

If there are constraints, however, propagation from pi(t) to p̃it+12δt is necessary to evaluate Khalf(t) because pi(t)p̃it12δt=p̃it+12δtpi(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 Khalf(t) and estimate instantaneous temperature Topt(t) and pressure Phalf(t). Based on the estimated Topt(t) and Phalf(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 Khalf(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 NG 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:

Mg=igmi,Rg=igmiri/Mg,Pg=igpi,Fg=igfi.
(28)

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

rgi=riRg,pgi=pimiMgPg.
(29)

In terms of group coordinate Rg and momentum Pg, temperature evaluated as NfkBT=ipitṙit is expressed as

NfkBT=gPgṘg+gigpgiṙgi.
(30)

Let us assume that N1 and N2 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 N1 and N2 is equal to Nf, and temperatures evaluated by gPgṘg and gigpgiṙgi are equal to each other from the generalized equipartition theorem,

T=1N1kBgPgṘg=1N2kBgigpgiṙgi.
(31)

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

T=1N1kBgPgṘg.
(32)

In velocity Verlet integration, gPgṘg is equal to

gPg(t)Ṙg(t)=23Kfull,g(t)+43Khalf,g(t),
(33)

where

Kfull,g(t)=gPg(t)22Mg,Khalf,g(t)=12gP̃gt12δt22Mg+P̃gt+12δt22Mg.
(34)

Accordingly, we update the thermostat by evaluating instantaneous temperature as

T(t)=1N1kB23Kfull,g(t)+43Khalf,g(t).
(35)

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

P(t)=13VgRgFg+2Khalf,g(t).
(36)

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 pi(t) to p̃it+12Δt but just using

Khalf,g(t)=g12MgPg(t)2+Pg(t)P̃gt12δt2,
(37)

and the iteration procedure is not needed for the evaluation of Khalf,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,

pipi+evηδt1miMgpg,
(38)

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.1 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.8 The update of momenta and position vectors is then expressed as

pit+12δt=pit+eαvεδt/21miMgPg(t)+δt2fit,rit+δt=rit+evεδt1Rg(t)+δt2mipit+12δt.
(39)

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: fi = fi,fast + fi,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 is13 

eiLslowΔt/2eiLfastδtneiLslowΔt/2,
(40)

where

iLslow=ifi,slowpi,iLfast=ifi,fastpi+pimiri.
(41)

Evaluations of T/P values required for thermostat and barostat updates are not straightforward because of Khalf(t) evaluation at time t, which requires propagation of momenta and position vectors. Evaluation of Khalf(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̃gt12δt is replaced by

Pg(t)P̃gt12δt=δt2Fg(t)
(42)

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

Khalf,g(t)=g12MgPg(t)2+δt24Flong,g(t)+Fshort,g(t)2,
(43)

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:

eiLΔt=eiLB2Δtbaro/2eiLTΔttherm/2eiLPδt/2eiLB1δt×eiLRδteiLPδt/2aeiLTΔttherm/2beiLB2Δtbaro/2,
(44)

where Δtbaro, Δttherm, and Δt are related by Δtbaro = bΔttherm and Δttherm = 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

eiLΔt=eiLB2Δtbaro/2eiLTΔttherm/2eiLslowΔt/2eiLB1δteiLfastδtn×eiLslowΔt/2aeiLTΔttherm/2beiLB2Δtbaro/2.
(45)

A pure water box consists of 9250 water molecules in a cubic box. The interactions are described with the CHARMM TIP3P model.21 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.22 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., Δttb = 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.25 To test temperature and pressure evaluations, we define instantaneous temperature and pressure with the full time step,

Tfull(t)=Kfull(t)NfkB,Pfull(t)=13Viri(t)Fi(t)+2Kfull(t),
(46)

as well as our suggestions, Topt(t) and Phalf(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: Topt/Phalf, Topt/Pfull, Tfull/Phalf, and Tfull/Pfull. In all cases, stochastic velocity rescaling (SVR) thermostats15 with the MTK barostat in isotropic pressure coupling19 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, Tfull/Phalf and Topt/Pfull 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 Tfull/Phalf 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. Topt/Phalf has consistent results irrespective of group T/P or time steps. Tfull/Pfull 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 Tfull and pressure overestimation from Pfull. It can be understood from the MD simulations of heterogeneous system with anisotropic pressure evaluation, as shown in Sec. III. D 2.

TABLE I.

Density of a pure water box in isothermal–isobaric MD simulations.

T/P evaluationT/P styleδt (fs)Density (kg/m3)
Tfull/Pfull Atomic T/P 1009.73 ± 0.15 
Tfull/Pfull Atomic T/P 1008.37 ± 0.20 
Tfull/Phalf Atomic T/P 1008.84 ± 0.17 
Tfull/Phalf Atomic T/P 1002.45 ± 0.16 
Topt/Pfull Atomic T/P 1010.81 ± 0.14 
Topt/Pfull Atomic T/P 1015.02 ± 0.14 
Topt/Phalf Atomic T/P 1009.84 ± 0.17 
Topt/Phalf Atomic T/P 1009.30 ± 0.17 
Tfull/Pfull Group T/P 1009.81 ± 0.21 
Tfull/Pfull Group T/P 1009.27 ± 0.13 
Tfull/Phalf Group T/P 1009.65 ± 0.23 
Tfull/Phalf Group T/P 1007.95 ± 0.24 
Topt/Pfull Group T/P 1009.90 ± 0.24 
Topt/Pfull Group T/P 1009.65 ± 0.19 
Topt/Phalf Group T/P 1009.94 ± 0.19 
Topt/Phalf Group T/P 1009.10 ± 0.22 
T/P evaluationT/P styleδt (fs)Density (kg/m3)
Tfull/Pfull Atomic T/P 1009.73 ± 0.15 
Tfull/Pfull Atomic T/P 1008.37 ± 0.20 
Tfull/Phalf Atomic T/P 1008.84 ± 0.17 
Tfull/Phalf Atomic T/P 1002.45 ± 0.16 
Topt/Pfull Atomic T/P 1010.81 ± 0.14 
Topt/Pfull Atomic T/P 1015.02 ± 0.14 
Topt/Phalf Atomic T/P 1009.84 ± 0.17 
Topt/Phalf Atomic T/P 1009.30 ± 0.17 
Tfull/Pfull Group T/P 1009.81 ± 0.21 
Tfull/Pfull Group T/P 1009.27 ± 0.13 
Tfull/Phalf Group T/P 1009.65 ± 0.23 
Tfull/Phalf Group T/P 1007.95 ± 0.24 
Topt/Pfull Group T/P 1009.90 ± 0.24 
Topt/Pfull Group T/P 1009.65 ± 0.19 
Topt/Phalf Group T/P 1009.94 ± 0.19 
Topt/Phalf Group T/P 1009.10 ± 0.22 

1. Test conditions

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

Scd=3cos2θ12,
(47)

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

2. Temperature estimation from NVT simulations

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 Tfull 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.

TABLE II.

Solute and solvent temperatures of 80 DPPC in NVT MD simulations (unit is K).

TemperatureTemperature
EvaluationStyleThermostat(solvent)(solute)
Tfull Atomic Berendsen 326.99 ± 0.18 320.54 ± 0.10 
Tfull Atomic NHC 327.03 ± 0.32 320.49 ± 0.22 
Tfull Atomic SVR 326.98 ± 0.34 320.59 ± 0.26 
Tfull Group Berendsen 325.32 ± 0.19 321.55 ± 0.14 
Tfull Group NHC 325.24 ± 0.22 321.52 ± 0.27 
Tfull Group SVR 325.29 ± 0.40 321.67 ± 0.38 
Topt Atomic Berendsen 323.05 ± 0.16 323.13 ± 0.10 
Topt Atomic NHC 322.96 ± 0.28 323.03 ± 0.28 
Topt Atomic SVR 323.03 ± 0.36 323.16 ± 0.25 
Topt Group Berendsen 322.97 ± 0.23 323.12 ± 0.15 
Topt Group NHC 323.16 ± 0.21 323.09 ± 0.24 
Topt Group SVR 323.12 ± 0.40 323.13 ± 0.33 
TemperatureTemperature
EvaluationStyleThermostat(solvent)(solute)
Tfull Atomic Berendsen 326.99 ± 0.18 320.54 ± 0.10 
Tfull Atomic NHC 327.03 ± 0.32 320.49 ± 0.22 
Tfull Atomic SVR 326.98 ± 0.34 320.59 ± 0.26 
Tfull Group Berendsen 325.32 ± 0.19 321.55 ± 0.14 
Tfull Group NHC 325.24 ± 0.22 321.52 ± 0.27 
Tfull Group SVR 325.29 ± 0.40 321.67 ± 0.38 
Topt Atomic Berendsen 323.05 ± 0.16 323.13 ± 0.10 
Topt Atomic NHC 322.96 ± 0.28 323.03 ± 0.28 
Topt Atomic SVR 323.03 ± 0.36 323.16 ± 0.25 
Topt Group Berendsen 322.97 ± 0.23 323.12 ± 0.15 
Topt Group NHC 323.16 ± 0.21 323.09 ± 0.24 
Topt Group SVR 323.12 ± 0.40 323.13 ± 0.33 

3. Test results of NPT simulations of DPPC systems

a. Effect of the group T/P approach.

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 Tfull/Pfull and Tfull/Phalf 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 Tfull/Phalf and Topt/Phalf, one can understand that the group temperature approach reduces the errors of temperature evaluations, too. When evaluating temperature and pressure in conventional ways, δt2 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 Topt/Phalf evaluations.

FIG. 1.

(a) Order parameter of the first chain and (b) area per lipid (APL) with δt = 2 fs of the velocity Verlet Integrator.

FIG. 1.

(a) Order parameter of the first chain and (b) area per lipid (APL) with δt = 2 fs of the velocity Verlet Integrator.

Close modal
TABLE III.

Area per lipid (APL, unit is Å2) and membrane thickness (unit is Å) of the 80 DPPC system in isothermal–isobaric MD simulations with the SVR thermostat and MTK barostat.

T/P evaluationT/P styleδt (fs)APLThickness
Tfull/Pfull Atomic T/P 64.57 ± 0.39 37.93 ± 0.17 
Tfull/Pfull Atomic T/P 76.52 ± 0.50 33.00 ± 0.17 
Tfull/Phalf Atomic T/P 62.43 ± 0.47 39.17 ± 0.22 
Tfull/Phalf Atomic T/P 65.16 ± 0.25 38.35 ± 0.13 
Topt/Phalf Atomic T/P 61.44 ± 0.57 39.53 ± 0.26 
Topt/Phalf Atomic T/P 61.21 ± 0.83 39.57 ± 0.39 
Tfull/Pfull Group T/P 61.96 ± 0.71 39.27 ± 0.35 
Tfull/Pfull Group T/P 66.67 ± 0.40 37.17 ± 0.16 
Tfull/Phalf Group T/P 61.88 ± 0.68 39.36 ± 0.33 
Tfull/Phalf Group T/P 64.09 ± 0.39 38.71 ± 0.20 
Topt/Phalf Group T/P 61.27 ± 0.26 39.59 ± 0.12 
Topt/Phalf Group T/P 61.11 ± 0.57 39.63 ± 0.27 
T/P evaluationT/P styleδt (fs)APLThickness
Tfull/Pfull Atomic T/P 64.57 ± 0.39 37.93 ± 0.17 
Tfull/Pfull Atomic T/P 76.52 ± 0.50 33.00 ± 0.17 
Tfull/Phalf Atomic T/P 62.43 ± 0.47 39.17 ± 0.22 
Tfull/Phalf Atomic T/P 65.16 ± 0.25 38.35 ± 0.13 
Topt/Phalf Atomic T/P 61.44 ± 0.57 39.53 ± 0.26 
Topt/Phalf Atomic T/P 61.21 ± 0.83 39.57 ± 0.39 
Tfull/Pfull Group T/P 61.96 ± 0.71 39.27 ± 0.35 
Tfull/Pfull Group T/P 66.67 ± 0.40 37.17 ± 0.16 
Tfull/Phalf Group T/P 61.88 ± 0.68 39.36 ± 0.33 
Tfull/Phalf Group T/P 64.09 ± 0.39 38.71 ± 0.20 
Topt/Phalf Group T/P 61.27 ± 0.26 39.59 ± 0.12 
Topt/Phalf Group T/P 61.11 ± 0.57 39.63 ± 0.27 
b. Effect of T/P evaluations in MD integration with a large time step.

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 Topt/Phalf 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 Tfull/Phalf or Tfull/Pfull 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 Tfull/Pfull 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 Tfull/Pfull evaluations overestimates the APL and underestimates thickness, and the accurate volume estimation is obtained just by compensating the over- and under-estimated values.3 

FIG. 2.

Distributions of various physical properties of the 80 DPPC system on the basis of three different T/P evaluations (Tfull/Pfull, Tfull/Phalf, and Topt/Phalf). In each T/P evaluations, four simulation conditions (VV_2 fs, VV_2 fs (+group), VV_5 fs_HMR (+atomic), and VV_5 fs_HMR) are tested. (a) Order parameter of carbon atoms in the first chain of DPPC, (b) volume, (c) APL, and (d) membrane thickness.

FIG. 2.

Distributions of various physical properties of the 80 DPPC system on the basis of three different T/P evaluations (Tfull/Pfull, Tfull/Phalf, and Topt/Phalf). In each T/P evaluations, four simulation conditions (VV_2 fs, VV_2 fs (+group), VV_5 fs_HMR (+atomic), and VV_5 fs_HMR) are tested. (a) Order parameter of carbon atoms in the first chain of DPPC, (b) volume, (c) APL, and (d) membrane thickness.

Close modal
c. Effect of thermostat algorithms.

To understand the effect of different thermostat algorithms using the same T/P evaluations, we compared three thermostats: stochastic velocity rescaling (SVR),15 Berendsen,14 and Nose–Hoover chain16 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.

FIG. 3.

Temperature distributions of the 80 DPPC system in NPT MD simulations using atomic (left) and group (right) T/P. Three thermostats, namely, stochastic velocity rescaling (SVR), Nose–Hoover chain (NHC), and Berendsen thermostats are combined with the MTK barostat in isothermal–isobaric MD simulations. For NHC and Berendsen, thermostats are marked as “(+NHC)” and “(+Berendsen),” respectively.

FIG. 3.

Temperature distributions of the 80 DPPC system in NPT MD simulations using atomic (left) and group (right) T/P. Three thermostats, namely, stochastic velocity rescaling (SVR), Nose–Hoover chain (NHC), and Berendsen thermostats are combined with the MTK barostat in isothermal–isobaric MD simulations. For NHC and Berendsen, thermostats are marked as “(+NHC)” and “(+Berendsen),” respectively.

Close modal
TABLE IV.

Area per lipid (APL, unit: Å2) and membrane thickness (unit is Å) of the 80 DPPC system with three different thermostats.

δtΔttb
ThermostatT/P style(fs)(fs)APLThickness
SVR Atomic T/P 61.44 ± 0.57 39.53 ± 0.26 
SVR Group T/P 61.27 ± 0.26 39.59 ± 0.12 
SVR Atomic T/P 61.21 ± 0.83 39.57 ± 0.39 
SVR Group T/P 61.11 ± 0.57 39.63 ± 0.27 
SVR Atomic T/P 20 61.00 ± 0.48 39.68 ± 0.21 
SVR Group T/P 20 61.11 ± 0.32 39.64 ± 0.16 
NHC Atomic T/P 61.18 ± 0.62 39.60 ± 0.29 
NHC Group T/P 61.43 ± 0.30 39.50 ± 0.14 
NHC Atomic T/P 61.30 ± 0.49 39.54 ± 0.24 
NHC Group T/P 60.93 ± 0.48 39.71 ± 0.25 
NHC Atomic T/P 20 61.36 ± 0.36 39.49 ± 0.19 
NHC Group T/P 20 61.23 ± 0.57 39.59 ± 0.29 
Berendsen Atomic T/P 61.77 ± 0.53 39.36 ± 0.26 
Berendsen Group T/P 61.71 ± 0.37 39.36 ± 0.16 
Berendsen Atomic T/P 61.15 ± 0.62 39.60 ± 0.31 
Berendsen Group T/P 60.89 ± 0.77 39.76 ± 0.38 
Berendsen Atomic T/P 20 60.86 ± 0.64 39.76 ± 0.30 
Berendsen Group T/P 20 61.04 ± 0.53 39.66 ± 0.25 
δtΔttb
ThermostatT/P style(fs)(fs)APLThickness
SVR Atomic T/P 61.44 ± 0.57 39.53 ± 0.26 
SVR Group T/P 61.27 ± 0.26 39.59 ± 0.12 
SVR Atomic T/P 61.21 ± 0.83 39.57 ± 0.39 
SVR Group T/P 61.11 ± 0.57 39.63 ± 0.27 
SVR Atomic T/P 20 61.00 ± 0.48 39.68 ± 0.21 
SVR Group T/P 20 61.11 ± 0.32 39.64 ± 0.16 
NHC Atomic T/P 61.18 ± 0.62 39.60 ± 0.29 
NHC Group T/P 61.43 ± 0.30 39.50 ± 0.14 
NHC Atomic T/P 61.30 ± 0.49 39.54 ± 0.24 
NHC Group T/P 60.93 ± 0.48 39.71 ± 0.25 
NHC Atomic T/P 20 61.36 ± 0.36 39.49 ± 0.19 
NHC Group T/P 20 61.23 ± 0.57 39.59 ± 0.29 
Berendsen Atomic T/P 61.77 ± 0.53 39.36 ± 0.26 
Berendsen Group T/P 61.71 ± 0.37 39.36 ± 0.16 
Berendsen Atomic T/P 61.15 ± 0.62 39.60 ± 0.31 
Berendsen Group T/P 60.89 ± 0.77 39.76 ± 0.38 
Berendsen Atomic T/P 20 60.86 ± 0.64 39.76 ± 0.30 
Berendsen Group T/P 20 61.04 ± 0.53 39.66 ± 0.25 
d. Effect of integrators (velocity Verlet vs MTS integrators).

For more extensive tests, we compared the properties of a larger DPPC system: 160 DPPC based on Topt/Phalf. Because of possibilities of changed average results, we performed several MD simulations under the same conditions. The structural properties are written in Table V. We found that δt = 5 fs with the conventional integrator or δt = 3.5 fs with the MTS integrator gives almost similar results to δt = 2 fs. Based on accurate T/P evaluations, group T/P or HMR does not give any effect in average properties. δt = 4 fs with the MTS integrator was also tested, and we observed that it slightly underestimates the APL and overestimates lipid thickness. The amount of error is not significant, but we suggest using δt = 3.5 fs for better reliability.

TABLE V.

Area per lipid (APL, unit: Å2) and lipid thickness (unit is Å) of the 160 DPPC system with various simulation conditions.

Test schemeRun nos.T/P styleδt (fs)APLThickness
VV_2 fs Atomic 61.44 ± 0.30 39.50 ± 0.16 
 Atomic 61.46 ± 0.21 39.48 ± 0.09 
 Atomic 61.65 ± 0.26 39.38 ± 0.18 
 Atomic 61.72 ± 0.27 39.35 ± 0.13 
 Atomic 61.57 ± 0.16 39.44 ± 0.08 
VV_2 fs (+group) Group 61.54 ± 0.29 39.45 ± 0.15 
 Group 61.36 ± 0.33 39.51 ± 0.17 
VV_2 fs_HMR Group 61.67 ± 0.29 39.37 ± 0.15 
 Group 61.46 ± 0.42 39.47 ± 0.22 
 Group 61.43 ± 0.38 39.47 ± 0.18 
 Group 61.32 ± 0.41 39.53 ± 0.20 
 Group 61.52 ± 0.23 39.45 ± 0.11 
VV_5 fs_HMR (+atomic) Atomic 61.28 ± 0.42 39.53 ± 0.21 
 Atomic 61.43 ± 0.37 39.49 ± 0.17 
VV_5 fs_HMR Group 61.13 ± 0.24 39.60 ± 0.12 
 Group 61.08 ± 0.42 39.61 ± 0.22 
 Group 61.15 ± 0.29 39.57 ± 0.14 
 Group 61.58 ± 0.45 39.36 ± 0.21 
 Group 61.27 ± 0.43 39.52 ± 0.21 
 Group 60.92 ± 0.40 39.69 ± 0.20 
MTS_3.5 fs_HMR Group 3.5 60.92 ± 0.42 39.73 ± 0.20 
 Group 3.5 60.78 ± 0.39 39.79 ± 0.20 
 Group 3.5 60.74 ± 0.40 39.80 ± 0.21 
 Group 3.5 60.82 ± 0.51 39.78 ± 0.26 
 Group 3.5 60.87 ± 0.42 39.74 ± 0.21 
MTS_4 fs_HMR Group 60.39 ± 0.48 39.97 ± 0.25 
Test schemeRun nos.T/P styleδt (fs)APLThickness
VV_2 fs Atomic 61.44 ± 0.30 39.50 ± 0.16 
 Atomic 61.46 ± 0.21 39.48 ± 0.09 
 Atomic 61.65 ± 0.26 39.38 ± 0.18 
 Atomic 61.72 ± 0.27 39.35 ± 0.13 
 Atomic 61.57 ± 0.16 39.44 ± 0.08 
VV_2 fs (+group) Group 61.54 ± 0.29 39.45 ± 0.15 
 Group 61.36 ± 0.33 39.51 ± 0.17 
VV_2 fs_HMR Group 61.67 ± 0.29 39.37 ± 0.15 
 Group 61.46 ± 0.42 39.47 ± 0.22 
 Group 61.43 ± 0.38 39.47 ± 0.18 
 Group 61.32 ± 0.41 39.53 ± 0.20 
 Group 61.52 ± 0.23 39.45 ± 0.11 
VV_5 fs_HMR (+atomic) Atomic 61.28 ± 0.42 39.53 ± 0.21 
 Atomic 61.43 ± 0.37 39.49 ± 0.17 
VV_5 fs_HMR Group 61.13 ± 0.24 39.60 ± 0.12 
 Group 61.08 ± 0.42 39.61 ± 0.22 
 Group 61.15 ± 0.29 39.57 ± 0.14 
 Group 61.58 ± 0.45 39.36 ± 0.21 
 Group 61.27 ± 0.43 39.52 ± 0.21 
 Group 60.92 ± 0.40 39.69 ± 0.20 
MTS_3.5 fs_HMR Group 3.5 60.92 ± 0.42 39.73 ± 0.20 
 Group 3.5 60.78 ± 0.39 39.79 ± 0.20 
 Group 3.5 60.74 ± 0.40 39.80 ± 0.21 
 Group 3.5 60.82 ± 0.51 39.78 ± 0.26 
 Group 3.5 60.87 ± 0.42 39.74 ± 0.21 
MTS_4 fs_HMR Group 60.39 ± 0.48 39.97 ± 0.25 

MD simulations of biological molecules have become more popular tools than before, for investigating the structure–dynamics–function relationship. In recent applications, longer time scales of motions for larger biomolecules, such as membrane proteins in an explicit lipid bilayer, large protein/nucleic acid complexes such as a ribosome, RNA polymerase, nucleosomes, and even a part of chromatins, have become the target of atomistic MD simulations. To pursue the realism of cellular environments has also been a new trend in MD simulations using the latest supercomputers. However, in all cases, insufficient conformational sampling of biomolecules in these simulations have caused serious problems if one aims at comparing the simulation results with experimental data. As we discussed in this work, MD simulation lengths depend on a time step, δt and the number of iterations in the time propagator. So, stable and reliable MD simulations under realistic conditions, namely, isothermal and isobaric conditions, are essential for all the biological applications. Because of this consideration, we have focused on the development of novel MD integrators, which can be used with a larger time step while keeping numerical stability and reproducing physical properties of target molecular systems.

In this work, we introduced the group-based approach for evaluating the instantaneous temperature and pressure of a simulation system and combined this with our previous MD integrator where the accuracy is kept up to the third order of a time step, δt. As we discussed in the Introduction, system temperature and pressure affect structural and dynamical properties of target simulation systems as well as the numerical stability of MD integration with a larger time step, δt. Since we derive the mathematical formulation starting from the Boltzmann distribution, the MD integrator proposed in the work has a solid theoretical base compared to the conventional ones. In addition, the extensive MD simulations of 80 DPPC and 160 DPPC systems in various simulation conditions suggest how robust the proposed MD integrator is and how sensitive the conventional ones are.

The tested conditions in the works are different T/P definitions, atomic vs group T/P, different time steps (δt = 2.0 fs, or 5.0 fs), different integrators (velocity Verlet or r-RESPA), different thermostat algorithms (Berendsen,14 SVR,15 and NHC,16), and with/without HMR.17 As shown in Tables III–V, a large number of combinations of different simulation conditions were tested for 80 DPPC and 160 DPPC systems. We observed that, in particular, the evaluation forms of temperature and pressure, as we derived previously,5,7 affect APL and membrane thickness significantly: Topt/Phalf produces the accurate APL and thickness, which are robust for the other conditions, such as atomic vs group T/P and a short vs longer time step. The other choices are more sensitive to these conditions.

By introducing the group-based T/P evaluations, an MTS integrator is now available in the proposed MD integrator for isothermal and isobaric conditions. The usage of an MTS integrator implies several practical advantages in MD simulations. For instance, if one uses hybrid CPU + GPU computers using GENESIS MD software, the performance of MD simulations with the MTS integrator is much better.27,28 In GENESIS, slow forces related to the reciprocal-space non-bonded calculations are assigned to the CPU, while the real-space non-bonded interactions are computed using graphics processing unit (GPU). In this implementation, the CPU calculations could be a bottleneck in the total cost, while MD simulation based on MTS is able to avoid the CPU computation every other step. So, the current development allows us not only to do the accurate MD simulations but also to perform them more efficiently. Such computational bases are important for longer time scale dynamics of large biomolecular systems.

We will further test the MD integrator for other biological systems, such as membrane proteins with an explicit lipid bilayer, nucleic acids, and so on. In addition, for free-energy calculations based on MD simulations, we need to apply enhanced conformational sampling algorithms, such as replica-exchange MD (T-REMD),29 replica-exchange umbrella sampling (REUS),30 Gaussian accelerated replica-exchange umbrella sampling (GaREUS),31 replica exchange with solute tempering (REST/REST2/gREST),32–34 and so on. Such comprehensive tests are carried out as a separated work and will be published elsewhere.

This research used in-house computer resources and the RIKEN HOKUSAI supercomputer. The computer resources of Oakforest-PACS were also provided through the HPCI System Research project (Project IDs: hp190097, hp200129, and hp200135). The research was supported in part by the MEXT as “FLAGSHIP 2020 project,” “Priority Issue on Post-K computer” (Building Innovative Drug Discovery Infrastructure Through Functional Control of Biomolecular Systems), “Program for Promoting Researches on the Supercomputer Fugaku” (Biomolecular dynamics in a living cell/MD-driven Precision Medicine), and MEXT/KAKENHI (Grant No. 19H05645) (to Y.S.).

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

1.
M. E.
Tuckerman
,
Statistical Mechanics: Theory and Molecular Simulation
(
Oxford University Press
,
New York
,
2010
).
2.
A. L.
Cheng
and
K. M.
Merz
,
J. Phys. Chem.
100
,
1927
(
1996
).
3.
W. F.
van Gunsteren
,
X.
Daura
,
N.
Hansen
,
A. E.
Mark
,
C.
Oostenbrink
,
S.
Riniker
, and
L. J.
Smith
,
Angew. Chem., Int. Ed.
57
,
884
(
2018
).
4.
M. P.
Eastwood
,
K. A.
Stafford
,
R. A.
Lippert
,
M. Ø.
Jensen
,
P.
Maragakis
,
C.
Predescu
,
R. O.
Dror
, and
D. E.
Shaw
,
J. Chem. Theory Comput.
6
,
2045
(
2010
).
5.
J.
Jung
,
C.
Kobayashi
, and
Y.
Sugita
,
J. Chem. Theory Comput.
15
,
84
94
(
2019
).
6.
7.
J.
Jung
,
C.
Kobayashi
, and
Y.
Sugita
,
J. Chem. Phys.
148
,
164109
(
2018
).
8.
R. A.
Lippert
,
C.
Predescu
,
D. J.
Ierardi
,
K. M.
Mackenzie
,
M. P.
Eastwood
,
R. O.
Dror
, and
D. E.
Shaw
,
J. Chem. Phys.
139
,
164106
(
2013
).
9.
M.
Di Pierro
,
R.
Elber
, and
B.
Leimkuhler
,
J. Chem. Theory Comput.
11
,
5624
(
2015
).
10.
M.
Ikeguchi
,
J. Comput. Chem.
25
,
529
(
2004
).
11.
J.-P.
Ryckaert
,
G.
Ciccotti
, and
H. J. C.
Berendsen
,
J. Comput. Phys.
23
,
327
(
1977
).
12.
H. C.
Andersen
,
J. Comput. Phys.
52
,
24
(
1983
).
13.
M.
Tuckerman
,
B. J.
Berne
, and
G. J.
Martyna
,
J. Chem. Phys.
97
,
1990
(
1992
).
14.
H. J. C.
Berendsen
,
J. P. M.
Postma
,
W. F.
van Gunsteren
,
A.
Dinola
, and
J. R.
Haak
,
J. Chem. Phys.
81
,
3684
(
1984
).
15.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
,
J. Chem. Phys.
126
,
014101
(
2007
).
16.
G. J.
Martyna
,
M. L.
Klein
, and
M.
Tuckerman
,
J. Chem. Phys.
97
,
2635
(
1992
).
17.
K. A.
Feenstra
,
B.
Hess
, and
H. J. C.
Berendsen
,
J. Comput. Chem.
20
,
786
(
1999
).
18.
C. W.
Hopkins
,
S.
Le Grand
,
R. C.
Walker
, and
A. E.
Roitberg
,
J. Chem. Theory Comput.
11
,
1864
(
2015
).
19.
G. J.
Martyna
,
M. E.
Tuckerman
,
D. J.
Tobias
, and
M. L.
Klein
,
Mol. Phys.
87
,
1117
(
1996
).
20.
G.
Bussi
,
T.
Zykova-Timan
, and
M.
Parrinello
,
J. Chem. Phys.
130
,
074101
(
2009
).
21.
D. J.
Price
and
C. L.
Brooks
,
J. Chem. Phys.
121
,
10096
(
2004
).
22.
E. L.
Wu
,
X.
Cheng
,
S.
Jo
,
H.
Rui
,
K. C.
Song
,
E. M.
Dávila-Contreras
,
Y. F.
Qi
,
J. M.
Lee
,
V.
Monje-Galvan
,
R. M.
Venable
,
J. B.
Klauda
, and
W.
Im
,
J. Comput. Chem.
35
,
1997
(
2014
).
23.
J. B.
Klauda
,
R. M.
Venable
,
J. A.
Freites
,
J. W.
O’Connor
,
D. J.
Tobias
,
C.
Mondragon-Ramirez
,
I.
Vorobyov
,
A. D.
MacKerell
, and
R. W.
Pastor
,
J. Phys. Chem. B
114
,
7830
(
2010
).
24.
J.
Huang
and
A. D.
MacKerell
,
J. Comput. Chem.
34
,
2135
(
2013
).
25.
U.
Essmann
,
L.
Perera
,
M. L.
Berkowitz
,
T.
Darden
,
H.
Lee
, and
L. G.
Pedersen
,
J. Chem. Phys.
103
,
8577
(
1995
).
26.
S.
Miyamoto
and
P. A.
Kollman
,
J. Comput. Chem.
13
,
952
(
1992
).
27.
J.
Jung
,
A.
Naurse
,
C.
Kobayashi
, and
Y.
Sugita
,
J. Chem. Theory Comput.
12
,
4947
(
2016
).
28.
C.
Kobayashi
,
J.
Jung
,
Y.
Matsunaga
,
T.
Mori
,
T.
Ando
,
K.
Tamura
,
M.
Kamiya
, and
Y.
Sugita
,
J. Comput. Chem.
38
,
2193
(
2017
).
29.
Y.
Sugita
and
Y.
Okamoto
,
Chem. Phys. Lett.
314
,
141
(
1999
).
30.
Y.
Sugita
,
A.
Kitao
, and
Y.
Okamoto
,
J. Chem. Phys.
113
,
6042
(
2000
).
31.
H.
Oshima
,
S. Y.
Re
, and
Y.
Sugita
,
J. Chem. Theory Comput.
15
,
5199
(
2019
).
32.
P.
Liu
,
B.
Kim
,
R. A.
Friesner
, and
B. J.
Berne
,
Proc. Natl. Acad. Sci. U. S. A.
102
,
13749
(
2005
).
33.
M.
Kamiya
and
Y.
Sugita
,
J. Chem. Phys.
149
,
072304
(
2018
).
34.
L.
Wang
,
R. A.
Friesner
, and
B. J.
Berne
,
J. Phys. Chem. B
115
,
11305
(
2011
).