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.

## I. INTRODUCTION

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 *δ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^{11} or RATTLE^{12} 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.

## II. METHODS

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.

### A. Derivation of optimal temperature and pressure evaluations in MD

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

where $\beta =1kBT$, *k*_{B} is the Boltzmann constant, and *H* is the Hamiltonian. By replacing $\u2202H\u2202pie\u2212\beta H$ with $\u22121\beta \u2202e\u2212\beta H\u2202pi=\u2212kBT\u2202e\u2212\beta H\u2202pi$ and applying the integration by parts, $pi\u22c5r\u0307i$ and $ri\u22c5p\u0307i$ 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\u0303it+12\delta t$ is the intermediate momentum of the *i*th particle propagated from **p**_{i}(*t*) by $12\delta t$. We do not write it as $pit+12\delta t$ because $p\u0303it+12\delta t$ and Taylor expansion of $pit+12\delta 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\u0303it+12\delta 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\u0303it+12\delta t$ and $p\u0303it\u221212\delta t$ into two ways, i.e., Taylor expansion of $p\u0303it+12\delta t$ and $p\u0303it\u221212\delta t$ in terms of $p\u0303it$, and MD integration leads the following relationship between $p\u0303\u0307it$ and **F**_{i}(*t*):

Similarly, the Taylor expansion and MD integration expression for the sum of $p\u0303it+12\delta t$ and $p\u0303it\u221212\delta t$ relates **p**_{i}(*t*) and $p\u0303i(t)$ as follows:

By applying the time derivatives in Eq. (9) and canceling out $p\u0303\u0307it$, one can obtain the relationship between $p\u0307i(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\u221216ri(5)t\delta t2+O\delta t4$ from Eq. (11), so the *δt*^{2} term in Eq. (12) is

By making use of $pi(t)\u22c5pi(2)t=\u2212p\u0307i(t)2$ as well as the relationship between $p\u0307i(t)$ and **F**_{i}(*t*) in Eq. (10), *N*_{f}*k*_{B}*T* 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*^{2} 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.

### B. MD integration under the NVT/NPT conditions using optimal temperature and pressure evaluations

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}

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^{19} 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

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\u0303it+12\delta 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\u0303it+12\delta 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\u0303it+12\delta t$ is necessary to evaluate *K*_{half}(*t*) because $pi(t)\u2212p\u0303it\u221212\delta t=p\u0303it+12\delta t\u2212pi(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.

### C. MD integration with group pressure and temperature

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=\u2211ipit\u22c5r\u0307it$ is expressed as

Let us assume that *N*_{1} and *N*_{2} 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*_{1} and *N*_{2} is equal to *N*_{f}, and temperatures evaluated by $\u2211gPg\u22c5R\u0307g$ and $\u2211g\u2211i\u2208gpgi\u22c5r\u0307gi$ are equal to each other from the generalized equipartition theorem,

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

In velocity Verlet integration, $\u2211gPg\u22c5R\u0307g$ 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\u0303it+12\Delta 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.

### D. MD integration with group pressure and temperature

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 *i*th particle and $v\eta =p\eta Q$. Without group pressure, for given values of $v\epsilon =p\epsilon W$ and $\alpha =1+3Nf$, the update of momenta and position vectors in the first half-time step integration, corresponding to $eiLP\Delta t/2$ and $eiLR\Delta 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

### E. Multiple time step integration with group pressure and temperature

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 *n*th step, Δ*t* = *nδt*, and a time propagator at the increment of Δ*t* under the NVE condition is^{13}

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)\u2212P\u0303gt\u221212\delta 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

## III. RESULTS

### A. Test systems: A pure water box and DPPC lipid bilayer systems (80 DPPC and 160 DPPC)

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.

### B. Test conditions in MD simulations

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.^{25} 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*).

### C. Test of a pure water box

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^{15} with the MTK barostat in isotropic pressure coupling^{19} 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.

### D. Test of DPPC

#### 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

where $\theta $ 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 *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.

#### 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 *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*^{2} 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.

##### 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 *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.^{3}

##### 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 chain^{16} 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.

##### 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 *T*_{opt}/*P*_{half}. 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.

## IV. DISCUSSION AND CONCLUSION

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: *T*_{opt}/*P*_{half} 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.

## ACKNOWLEDGMENTS

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

## DATA AVAILABILITY

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