Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

Dominik Sidler; Marc Lehner; Simon Frasch; Michael Cristófol-Clough; Sereina Riniker

doi:10.12688/f1000research.16715.1

Home Browse Density artefacts at interfaces caused by multiple time-step effects...

ALL Metrics

Views

Downloads

Get PDF

Get XML

Export

▬

✚

Research Article

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

[version 1; peer review: 1 approved, 1 approved with reservations]

Dominik Sidler¹, Marc Lehner¹, Simon Frasch¹, Michael Cristófol-Clough¹, Sereina Riniker ¹

Dominik Sidler¹, Marc Lehner¹, [...] Simon Frasch¹, Michael Cristófol-Clough¹, Sereina Riniker ¹

PUBLISHED 05 Nov 2018

Author details Author details

¹ Laboratory of Physical Chemistry, ETH Zürich, Zurich, 8093, Switzerland

Dominik Sidler
Roles: Conceptualization, Data Curation, Investigation, Methodology, Validation, Visualization, Writing – Original Draft Preparation, Writing – Review & Editing

Marc Lehner
Roles: Investigation, Methodology, Validation, Visualization

Simon Frasch
Roles: Investigation, Methodology, Validation

Michael Cristófol-Clough
Roles: Investigation, Methodology, Validation

Sereina Riniker
Roles: Conceptualization, Funding Acquisition, Methodology, Project Administration, Resources, Supervision, Writing – Review & Editing

OPEN PEER REVIEW

REVIEWER STATUS

This article is included in the Cheminformatics gateway.

This article is included in the Mathematical, Physical, and Computational Sciences collection.

Abstract

Background: Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts.
Methods: The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution.
Results: In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces. The presence of the observed artefacts was found to be independent of the interface shape and the thermostatting method used. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs.
Conclusion: A SR scheme avoids numerical artefacts and offers an interesting alternative to TR RESPA with respect to performance optimization.

Keywords

Multiple Time-Step Integration, Molecular Dynamics, Resonance Artefacts, RESPA, Twin-range, Liquid Phase Interfaces

Corresponding author: Sereina Riniker

Competing interests: No competing interests were disclosed.

Grant information: The authors gratefully acknowledge financial support by
the Swiss National Science Foundation [200021-178762]
and by Eidgenössische Technische Hochschule (ETH) Zürich [ETH-34 17-2].

Copyright: © 2018 Sidler D et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Data associated with the article are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).

How to cite: Sidler D, Lehner M, Frasch S et al. Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations [version 1; peer review: 1 approved, 1 approved with reservations]. F1000Research 2018, 7:1745 (https://doi.org/10.12688/f1000research.16715.1) First published: 05 Nov 2018, 7:1745 (https://doi.org/10.12688/f1000research.16715.1) Latest published: 08 Mar 2019, 7:1745 (https://doi.org/10.12688/f1000research.16715.3)

Introduction

To describe the dynamic processes of biomolecules, such as proteins, DNA, carbohydrates and membranes, molecular dynamics (MD) simulations have been proven to be a powerful technique to provide insights at the atomic level. In MD simulations, physical processes that involve a wide range of time scales have to be considered appropriately. However, the accurate and efficient treatment of various time scales presents a notoriously difficult problem due to resonance artefacts arising from separation into distinct numerical subproblems^1,2. Over the past decades, various methods have been proposed to mediate the issue and to speed up the simulation substantially. For example, fast vibrational modes (e.g. covalent bonds) can be removed completely by constraining the motion, e.g. using SHAKE³, SETTLE⁴, M-SHAKE⁵, LINCS⁶, or FLEXSHAKE⁷. Alternatively, multiple time-step (MTS) integration is an efficient approach to accommodate motions on different time scales. In the twin-range (TR) scheme^8–10, the pairwise nonbonded forces $f_{i}^{nb}$ acting on particle i are split into a short-range contribution from particles j within the cutoff R_s, and a mid-range contribution from particles j between R_s and the long-range cutoff R_l (see Figure 1),

Figure 1. Simulated Systems.

Left: Schematic illustration of the twin-range (TR) setup with the inner cutoff R_s and the outer cutoff R_l. Interactions beyond R_l are not calculated explicitly. Middle: Schematic illustration of the planar layer systems. Right: Schematic illustration of the systems with a droplet in a water box.

f_{i}^{nb} (t) = f_{i}^{nb, s} (t) + f_{i}^{nb, m} (t) . (1)

with

f_{i}^{nb, s} (t) = \sum_{j \in N (r_{i j} \leq R_{s})} f_{i j}^{ele} + f_{i j}^{vdw}, (2)

f_{i}^{nb, m} (t) = \sum_{j \in N (R_{s} < r_{i j} \leq R_{l})} f_{i j}^{ele} + f_{i j}^{vdw} \cdot (3)

The short-range forces $f_{i}^{nb, s}$ are updated every time step ∆t (inner time step), while the mid-range forces $f_{i}^{nb, m}$ are updated every n_m-th time step, n_m∆t with n_m ≥ 1 (outer time step). Long-range electrostatic interactions beyond the cutoff radius R_l can be included using the (an)isotropic reaction-field methods (RF or ARF)^11,12, or lattice-sum methods^13,14. Although theoretically possible, the TR scheme is generally not used together with lattice-sum methods. Simulations with (A)RF typically use larger cutoffs R_l than simulations with PME (e.g. 1.4 nm for the GROMOS force field¹⁵ compared to 0.9 nm for the AMBER force field¹⁶), which in the past encouraged the use of a TR scheme to reduce the computational cost.

In practice, for each particle i two pairlists are generated, which track the neighbouring short-range and mid-range particles j. To be consistent, the period of updating both pairlists n_p ∆t should be equal to the period of updating the mid-range forces n_m∆t. In the following, the TR scheme is defined by R_s < R_l and n_m = n_p > 1, and the single-range (SR) scheme by R_s = R_l, n_m = 1 and n_p ≥ 1, i.e. there is only one pairlist for SR.

Splitting the forces into different contributions leaves some freedom of how to consider them in the integration. There are two main strategies for the inclusion of the outer contributions: (i) constant mid-range forces application (CFA) $f_{i}^{m}$ at every inner time step ∆t, or (ii) using a reversible reference system propagator algorithm (RESPA)^8–10, which is based on the Trotter factorisation of the Liouville propagator¹⁷. In the case of RESPA, scaled mid-range forces n_m $f_{i}^{m}$ are typically applied impulse-wise, i.e. only at every outer time step n_m∆t. Note that in this work the leap-frog scheme¹⁸ is used for the integration of the equations of motion, which is similar¹⁹ to velocity-Verlet²⁰ as typically employed with RESPA.

Although computationally efficient, MTS algorithms have the disadvantage of introducing numerical artefacts. Due to nonlinearity, it is notoriously difficult to formulate general physical or chemical conditions, which favour or suppress the build-up of these undesired resonances. Instead, specialised thermostats have been proposed as a remedy, which suppress resonance artefacts from MTS integration, for example by imposing isokinetic constraints on every degree of freedom^21–23, or by selectively targeting the high-frequency modes using a colored noise thermostat²⁴. In practice, however, standard thermostats like weak-coupling²⁵, Nosé-Hoover (chain)^26–28, or stochastic thermostats^29–31 remain widely used in combination with a moderate integration time-step difference between fast and slow motions, as this should limit resonance artefacts to an acceptable value.

In this work, artefacts due to MTS integration are investigated with the main focus on a water-chloroform layer system using different thermostats (Figure 1). This system represents a simplified version of a biological membrane, for which an accurate simulation is of high interest in biology. The impact of splitting the nonbonded interactions into a short-range and mid-range contribution on the density profile normal to the interface and the kinetic energy distribution is investigated. In addition, a planar system consisting of cyclohexane and water as well as corresponding droplet systems are investigated. These additional setups allow to assess the impact of different solvents as well as different interface geometries on the formation of MTS artefacts. Overall, the artefacts are larger than those observed for more well-behaved, homogeneous systems^32,33.

Methods

For all simulations, the GROMOS software package³⁴ version 1.4.0 was used with the GROMOS 53A6 force field¹⁵. Periodic boundary conditions were applied. Newton’s equations of motion were integrated using the leap-frog scheme¹⁸ with an inner time step of 2 fs. Different outer time steps for the pairlist update or mid-range force evaluation were used (i.e. n_m, n_p ∈ [1, 20]). If not indicated otherwise, the SR scheme corresponds to a midrange force evaluation every 2 fs (n_m = 1) and a pairlist update every 10 fs (n_p = 5). The TR scheme typically corresponds to a mid-range force evaluation in combination with a pairlist update every 10 fs (n_m = n_p = 5). If not stated otherwise, the CFA method was used for the mid-range force contributions. The inner cutoff radius R_s was set to 0.8 nm and the outer cutoff radius R_l to 1.4 nm. For all simulations with RF, a charge-group based pairlist algorithm and cutoff were used (each solvent molecule is a single charge group), whereas an atomic cutoff was used for the PME simulations. The center of mass motion was stopped every 1 ps. The coordinates were saved every 0.2 ps and the energies every 0.1 ps. Bond lengths were constrained using the SHAKE algorithm with a geometric tolerance of 10⁻⁴. Three different thermostats were tested, i.e. weak coupling²⁵ (WC), Nosé-Hoover^26,27 (NH,) and Nosé-Hoover chain²⁸ (NHchain). Different temperature baths were used for different solvents. The coupling time τ was set to 0.1 ps for all three thermostatting methods. A chain length of three was employed for the NHchain thermostat. A standard homogeneous RF approach¹¹ was used to mimic long-range electrostatic interactions beyond a charge-group based cutoff R_l . The dielectric permittivity was set to 78.5, which corresponds to the experimentally measured value of water³⁵.

Systems

The MTS investigation was performed mainly using a water-chloroform layer, in combination with three auxiliary systems consisting of a water-cyclohexane layer, a chloroform sphere in a water box, and a cyclohexane sphere in a water box. An illustration of the setups is shown in Figure 1. With these additional systems, different geometries of the interface as well as different chemical compositions were tested. An overview of the simulated systems can be found in Table 1. Corresponding snapshots are provided in Figure S1 (Supplementary File 1).

Table 1. Overview of the system dimensions and simulation times.

Setup	Phase 1 [#molecules]	Phase 2 [#molecules]	Box (x, y, z) [nm]	Time [ns]
Layer	2002 CHCl₃	8129 H₂O	(4.1, 4.1, 30.2)	4
Layer	6750 C₆H₁₂	39’366 H₂O	(8.5, 8.5, 33.9)	1
Droplet	847 CHCl₃	59’484 H₂O	(12.2, 12.2, 12.2)	1
Droplet	619 C₆H₁₂	59’018 H₂O	(12.2, 12.2, 12.2)	1

The main system consisted of 2002 chloroform (CHCl₃)³⁶ molecules and 8129 SPC³⁷ water molecules within an elongated rectangular box of 30.2 nm length in z-direction. Identical initial coordinates were used for all simulations, which were equilibrated for 1 ns using the SR scheme under NVT conditions. The simulation length of the production simulations was 1–4 ns, performed under NVT conditions. For the water-cyclohexane layer system, 6750 cyclohexane (C₆H₁₂) and 39’366 water molecules were used in an elongated rectangular box of 33.85 nm length in z-direction. This system was simulated under NVT as well as NPT conditions (using a Berendsen barostat²⁵ with a coupling time τ_b = 0.5 ps). The droplet systems consisted of a droplet of about 3 nm radius, solvated in a box of 12.2 nm length. This setup led to a droplet with 847 CHCl₃ molecules in 59’484 SPC water molecules, or a droplet with 619 C₆H₁₂ molecules solvated in 59’018 SPC water molecules, respectively. As the convergence towards the numerical artefacts was found to be very fast, all auxiliary systems were simulated without initial equilibration for 1 ns.

Analysis

For the spatially resolved density analysis, a program was developed using the data structures provided by the GROMOS++ collection of programs³⁸. The energy distributions were calculated by post-processing the output of the GROMOS++ program ene_ana accordingly. For the density calculations of the spherical interfaces, a slightly more sophisticated analysis approach was necessary. In order to obtain a radial density profile, the molecules were assigned to equidistant concentric spherical shells with respect to the center of geometry of the droplet. However, for the droplet it was difficult to define consistently a distance to the surface, because the interface may deviate substantially from its ideal spherical shape during simulation. In our analysis, the center of mass was chosen as an approximation of the center of geometry, which was then shifted (if necessary) in order to account for periodic boundary conditions, i.e. to ensure that the droplet remained at the center of the periodic box. Note that the volume of the spherical shells scales cubically with the radius, which increases the uncertainty of the estimator substantially close to the center of the droplet.

Results & discussion

The density profile of the water-chloroform system was determined along the z-axis normal to the layer. As can be seen in Figure 2, the splitting of the forces into a short-range and a less frequently updated mid-range contribution (given in Equation (1)) using the TR CFA scheme with n_m = n_p = 5 led to a significant density increase for chloroform close to the interface. The density increase was accompanied by a density decrease further away from the interface due to the constant volume conditions. The presence of an increase was found to be independent of the chosen thermostatting algorithm. An analysis of the density resonance pattern with respect to the outer time step n_m∆t is shown in Figure 3. Note that the density artefacts are not perfectly symmetric with respect to the layer norm due to finite simulation time. It can be seen that the TR CFA scheme leads to large resonance artefacts beyond n_m∆t ≈ 20 fs independent of the thermostatting method. On the other hand, evaluating the mid-range forces every step while keeping the pairlist update every fifth step (i.e. SR scheme with n_m = 1 and n_p = 5) resolved the issue almost entirely for all three investigated thermostats (Figure 2 and Figure 3, data not shown for NH and NHchain). The minor differences between SR schemes with different update frequencies can only be seen in density difference plots (Figure 4) with respect to the SR reference run (n_m = n_p = 1).

Figure 2. Profile of the density ρ in the water-chloroform layer system with respect to the layer normal z.

For the four setups with the twin-range (TR) scheme, the mid-range forces were updated with the same frequency as the pairlist (i.e. n_m = n_p = 5). For the single-range (SR) scheme, the mid-range forces were evaluated every time step ∆t (i.e. n_m = 1), whereas the pairlist was updated every 5th step (i.e. n_p = 5). The weak-coupling (WC), Nosé-Hoover (NH) or Nosé-Hoover chain (NHchain) thermostat was used.

Figure 3. Density profile with respect to the layer normal z and the pairlist update period n_p∆t.

For the single-range (SR) scheme, the mid-range forces were evaluated every time stepc ∆t (i.e. n_m = 1). For the three setups with the twin-range (TR) scheme, the mid-range forces were updated with the same period as the pairlist (i.e. n_p = n). The TR scheme was applied in a constant (CFA) or impulse-wise (RESPA) manner. The weak-coupling (WC) or Nosé-Hoover (NH) thermostat was used. Note that the density range is chosen to visualize the major artefacts occurring for n_m∆t > 15 fs. A visualisation of the minor artefacts around n_m∆t ≈ 10 fs can be seen in Figure 4.

Figure 4. Spacial density difference profiles with respect to a SR reference simulation (n_m = n_p = 1) under NVT conditions for layer systems (left) and droplet systems (right) with water-chloroform composition (top) or water-cyclohexane composition (bottom).

The Nosé-Hoover (NH) thermostat was used in all setups. The single-range (SR) scheme (red and blue lines) and the twin-range scheme with constant-force application (TR CFA, orange and green lines) and two different update frequencies (n_p = 5, 10) are compared. Running averages with a window of seven data points are shown.

These findings indicate that the observed density increase is mainly caused by resonance artefacts due to the CFA MTS integration. This argumentation is supported by the pairlist error analysis shown in Figure S2, Supplementary File 1. The change in the pairlist after 10 fs is only on the order of 1%, and, even more importantly, it scales sublinearly with n_m. Thus, the error from the pairlist cannot explain the substantial density increases seen in Figure 3. Instead, resonance effects lead to the build-up of the observed density artefacts over time, until a local equilibrium configuration is reached at the interface. This was found to be accompanied by a local change in the diffusion of the molecules at the interface. The mean residence time with respect to the layer normal z shows that regions of high density imply a locally reduced diffusivity of the molecules and vice versa (Figure S3, Supplementary File 1).

The numerical artefacts stem likely from an increasing mismatch between the forces and the position of the atoms the forces are applied to. In line with this interpretation, the density artefacts were suppressed for a typical value of the outer time step n_m∆t = 10 fs when the mid-range forces were applied impulse-wise (TR RESPA scheme, blue line in Figure 2), i.e. the direction of the forces and the atom positions matched at the time point when the forces were applied. However, the TR RESPA scheme does not eliminate the resonance artefacts but only shifts them to longer outer time steps in this context (Figure 3). Thus, the partitioning of the system into different time scales together with the choice of the time steps in the MTS integration remain delicate aspects of current MD implementations³⁹, also for impulse-wise MTS algorithms^1,2,40.

On the contrary, only small effects on the density profile could be observed for the SR scheme with pairlist update periods up to n_p∆t = 40 fs (Figure 3 and Figure 4). Due to the large cutoff R_l used here (i.e. 1.4 nm), the changes in the pairlist are small between updates and scale sub-linearly with n_m (Figure S2, Supplementary File 1). This observation is in line with the study of Krieger et al.⁴¹, where they found a negligible energy drift when using a similarly low pair-list update frequency in a homogeneous system. A rough performance analysis of the run time showed that for a realistic setup approximately half of the computational performance loss due to evaluating the mid-range forces at every time step could be compensated by a less frequent update of the pairlist (Table 2).

Table 2. Comparison of run times for three realistic simulation setups with different nonbonded force treatment: single-range scheme (SR) and twin-range scheme with constant-force application (TR CFA).

All systems were simulated for 0.4 ns on 1 CPU with 4 cores. A time step of 2 fs was used. The pairlist was generated using a charge-group based algorithm.

Force splitting	n_m	n_p	Run time [h]
TR CFA	5	5	6.1
SR	1	5	11.4
SR	1	20	8.8

Chemical composition and constant pressure

In order to exclude that the observed artefacts are specific to the water-chloroform system, the densities at a water-cyclohexane interface in a layer system were investigated. The corresponding density difference plots with respect to a SR reference are given in Figure 4. They show a similar behaviour of the TR scheme compared to the one observed for the water-chloroform systems, with slightly less pronounced artefacts. In addition, imposing constant pressure (NPT) conditions in the simulation setup does not reduce the observed density artefacts (Figure 5).

Figure 5. Density profile of the water-cyclohexane layer system under NVT conditions (left) or NPT conditions (right).

A semi-aniostropic pressure coupling was applied in x/y-direction using a Berendsen barostat, while the box length in z-direction remained fixed. A Nosé-Hoover (NH) thermostat was used in all setups. The single-range (SR) scheme (red, blue and orange lines) and the twin-range scheme with constant-force application (TR CFA, green and purple lines) and three different update frequencies (n_p = 1, 5, 10) are compared. Running averages with a window of seven data points are shown.

Geometry of the interface

In addition to planar interfaces, systems with a droplet of chloroform or cyclohexane in a water box were investigated. As can be seen in the right panels of Figure 4, density artefacts occur also at curved interfaces. Note that the corresponding plots contain more noise at small radii compared to the planar setups due to surface fluctuations, which influence the position of the centre of mass used as a reference. In addition, the surface fluctuations also give rise to a broadening of interface in terms of densities compared to the planar setup as can be seen in Figure S4 (Supplementary File 1).

Thermostatting method

The impact of the chosen thermostatting algorithm as well as the MTS integrator, and the energy distributions was investigated for the water-chloroform layer system. The WC thermostat does not give a canonical distribution of the kinetic energy by construction⁴², whereas the NH and NHchain thermostats do. In Figure 6, the kinetic energy distributions within water and chloroform are shown for the different setups. For chloroform, the difference in the width of the distribution between WC and NH(chain) thermostats is clearly visible, whereas the average kinetic energy remains the same, despite the significant density gradient in this solvent. In water, on the other hand, a clear shift in the kinetic energy distribution towards a higher average value was observed for the TR CFA scheme compared to the SR scheme with the WC thermostat. This shift was less pronounced with the NHchain thermostat, and absent with the NH thermostat. These observations indicate that a temperature shift may occur due to MTS resonances. However, as the density increase occurred with all three thermostatting methods, but the temperature shift only with two of them, the latter cannot be the cause for the former. Therefore, one can conclude that the resonance artefacts caused by the TR scheme cannot be removed with thermostats that control the velocities globally. For this, one would have to control the velocities selectively, e.g. by imposing isokinetic constraints or restraints on each degree of freedom^21–23.

Figure 6. Histogram of the kinetic energy for chloroform (left) and water (right) using a weak coupling (WC, black), Nosé-Hoover (NH, red) and Nosé-Hoover chain (NHchain, green) thermostat (bin size 1 K).

Note that a different number of molecules is simulated for chloroform (left) and water (right), which leads to a different width of the canonical distribution functions. For the twin-range (TR) scheme (solid lines), the mid-range forces were updated with the same period as the pairlist (i.e. n_m = n_p = 5). For the single-range (SR) scheme (dashed lines), the mid-range forces were evaluated every time step Δt (i.e. n_m = 1), whereas the pairlist was updated every 5th step (i.e. n_p = 5).

Dataset 1.Zip containing simulation input files.

Conclusions

The effect of using the TR scheme with a constant force application (CFA) at solvent-solvent interfaces was investigated with layer and droplet systems of different composition (water-chloroform and water-cyclohexane). These systems represent simplified versions of a biological membrane and a protein in solution. The force splitting was found to lead to substantial MTS resonance artefacts at the interfaces. It could be shown that the presence of the artificial density increase at the interfaces is independent of the thermostatting method, the interface geometry as well as the chemical composition of the system although the magnitude could vary. In addition, the diffusivity of the solvent molecules at the interface changed due to the density artefacts.

The resonance effects likely stem from a mismatch between the forces and the position of the atoms the forces are applied to. This hypothesis is supported by the fact that the pairlist changes only in the order of 1 % between updates and that the TR scheme with RESPA (where the forces are applied impulse-wise) does not eliminate the numerical artefacts but shifts them to longer time steps. Based on these findings, the use of TR CFA is not recommended for heterogeneous systems. Instead, the SR scheme should be used, which did not show any significant density errors in the entire range of pairlist updates (analysed up to n_p∆t = 40 fs). Thus, the increase in computational cost for SR can be partially compensated by updating the pairlist less frequently (i.e. n_p∆t > 10 fs). In addition, using only a single pairlist reduces the complexity with respect to implementation, as the additional difficulties arising from efficient memory management of two (or more) pairlists can be avoided.

Data availability

We are unable to provide the output files from the simulations directly due to file size. We instead provide the input files that when used as input with the GROMOS software package will generate the same results -

Dataset 1: Zip containing simulation input files 10.5256/f1000research.16715.d222949⁴³

Grant information

The authors gratefully acknowledge financial support by the Swiss National Science Foundation [200021-178762] and by Eidgenössische Technische Hochschule (ETH) Zürich [ETH-34 17-2].

Acknowledgements

The authors thank Philippe Hünenberger and Wilfred van Gunsteren for helpful discussions.

Supplementary material

Supplementary File 1: File containing the following supplementary figures:

Click here to access data

Figure S1: Snapshots of different simulation setups.
Figure S2: Change in mid-range pairlist with respect to update frequency.
Figure S3: Spatially resolved mean residence time.
Figure S4: Droplet density profiles.

Faculty Opinions recommended

References

1. Bishop TC, Skeel RD, Schulten K: Difficulties with multiple time stepping and fast multipole algorithm in molecular dynamics. J Comput Chem. 1997; 18(14): 1785–1791. Publisher Full Text
2. Ma Q, Izaguirre JA, Skeel RD: Verlet-I/r-RESPA/Impulse is limited by nonlinear instabilities. SIAM J Sci Comput. 2003; 24(6): 1951–1973. Publisher Full Text
3. Ryckaert JP, Ciccotti G, Berendsen HJ: Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J Comput Phys. 1977; 23(3): 327–341. Publisher Full Text
4. Miyamoto S, Kollman PA: SETTLE: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. J Comput Chem. 1992; 13(8): 952–962. Publisher Full Text
5. Kräutler V, van Gunsteren WF, Hünenberger PH: A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J Comput Chem. 2001; 22(5): 501–508. Publisher Full Text
6. Hess B, Bekker H, Berendsen HJC, et al.: LINCS: A linear constraint solver for molecular simulations. J Comput Chem. 1997; 18(12): 1463–1472. Publisher Full Text
7. Christen M, van Gunsteren WF: An approximate but fast method to impose flexible distance constraints in molecular dynamics simulations. J Chem Phys. 2005; 122(14): 144106. PubMed Abstract | Publisher Full Text
8. Berendsen HJ, Van Gunsteren WF, Zwinderman HR, et al.: Simulations of proteins in water. Ann N Y Acad Sci. 1986; 482(1): 269–286. PubMed Abstract | Publisher Full Text
9. Van Gunsteren WF, Berendsen HJ, Geurtsen RG, et al.: A molecular dynamics computer simulation of an eight-base-pair DNA fragment in aqueous solution: comparison with experimental two-dimensional NMR data. Ann N Y Acad Sci. 1986; 482(1): 287–303. PubMed Abstract | Publisher Full Text
10. van Gunsteren WF, Berendsen HJ: Computer simulation of molecular dynamics: Methodology, applications, and perspectives in chemistry. Angew Chem Int Ed. 1990; 29(9): 992–1023. Publisher Full Text
11. Tironi IG, Sperb R, Smith PE, et al.: A generalized reaction field method for molecular dynamics simulations. J Chem Phys. 1995; 102: 5451–5459. Publisher Full Text
12. Sidler D, Frasch S, Cristòfol-Clough M, et al.: Anisotropic reaction field correction for long-range electrostatic interactions in molecular dynamics simulations. J Chem Phys. 2018; 148(23): 234105. PubMed Abstract | Publisher Full Text
13. Hockney RW, Eastwood JW: Computer Simulation Using Particles. McGraw-Hill, New York, NY, 1981. Reference Source
14. Darden T, York D, Pedersen L: Particle mesh Ewald: An N log (N) method for Ewald sums in large systems. J Chem Phys. 1993; 98(12): 10089–10092. Publisher Full Text
15. Oostenbrink C, Villa A, Mark AE, et al.: A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. J Comput Chem. 2004; 25(13): 1656–1676. PubMed Abstract | Publisher Full Text
16. Wang J, Cieplak P, Kollman PA: How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J Comput Chem. 2000; 21(12): 1049–1074. Publisher Full Text
17. Tuckerman M, Berne BJ, Martyna GJ: Reversible multiple time scale molecular dynamics. J Chem Phys. 1992; 97: 1990–2001. Publisher Full Text
18. Hockney RW: Potential calculation and some applications. Technical report, Langley Research Center, Hampton, Va., 1970. Reference Source
19. Tuckerman M, Berne BJ, Martyna GJ: Reply to comment on: Reversible multiple time scale molecular dynamics. J Chem Phys. 1993; 99: 2278–2279. Publisher Full Text
20. Swope WC, Andersen HC, Berens PH, et al.: A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. J Chem Phys. 1982; 76: 637–649. Publisher Full Text
21. Minary P, Tuckerman M, Martyna GJ: Long time molecular dynamics for enhanced conformational sampling in biomolecular systems. Phys Rev Lett. 2004; 93(15): 150201. PubMed Abstract | Publisher Full Text
22. Leimkuhler B, Margul DT, Tuckerman ME: Stochastic, resonance-free multiple time-step algorithm for molecular dynamics with very large time steps. Mol Phys. 2013; 111(22–23): 3579–3594. Publisher Full Text
23. Margul DT, Tuckerman ME: A Stochastic, Resonance-Free Multiple Time-Step Algorithm for Polarizable Models That Permits Very Large Time Steps. J Chem Theory Comput. 2016; 12(5): 2170–2180. PubMed Abstract | Publisher Full Text
24. Morrone JA, Markland TE, Ceriotti M, et al.: Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances. J Chem Phys. 2011; 134(1): 014103. PubMed Abstract | Publisher Full Text
25. Berendsen HJC, Postma JPM, van Gunsteren WF, et al.: Molecular dynamics with coupling to an external bath. J Chem Phys. 1984; 81(8): 3684–3690. Publisher Full Text
26. Nosé S: A unified formulation of the constant temperature molecular dynamics methods. J Chem Phys. 1984; 81(1): 511–519. Publisher Full Text
27. Hoover WG: Canonical dynamics: Equilibrium phase-space distributions. Phys Rev A Gen Phys. 1985; 31(3): 1695–1697. PubMed Abstract | Publisher Full Text
28. Martyna GJ, Klein ML, Tuckerman M: Nosé–hoover chains: the canonical ensemble via continuous dynamics. J Chem Phys. 1992; 97(4): 2635–2643. Publisher Full Text
29. Schneider T, Stoll E: Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions. Phys Rev B. 1978; 17(3): 1302. Publisher Full Text
30. Andersen HC: Molecular dynamics simulations at constant pressure and/or temperature. J Chem Phys. 1980; 72(4): 2384–2393. Publisher Full Text
31. Bussi G, Donadio D, Parrinello M: Canonical sampling through velocity rescaling. J Chem Phys. 2007; 126(1): 014101. PubMed Abstract | Publisher Full Text
32. Schlick T, Mandziuk M, Skeel RD, et al.: Nonlinear resonance artifacts in molecular dynamics simulations. J Chem Phys. 1998; 140(1): 1–29. Publisher Full Text
33. Sandu A, Schlick T: Masking resonance artifacts in force-splitting methods for biomolecular simulations by extrapolative Langevin dynamics. J Chem Phys. 1999; 151(1): 74–113. Publisher Full Text
34. Schmid N, Christ CD, Christen M, et al.: Architecture, implementation and parallelisation of the GROMOS software for biomolecular simulation. Comput Phys Commun. 2012; 183(4): 890–903. Publisher Full Text
35. Haynes WM: CRC Handbook of Chemistry and Physics. CRC Press, Boca Raton, FL, 2014. Reference Source
36. Tironi IG, van Gunsteren WF: A molecular dynamics simulation study of chloroform. Mol Phys. 1994; 83(2): 381–403. Publisher Full Text
37. Berendsen HJC, Postma JPM, van Gunsteren WF, et al.: Interaction models for water in relation to protein hydration. In Intermolecular forces. Springer, 1981; 331–342. Publisher Full Text
38. Eichenberger AP, Allison JR, Dolenc J, et al.: GROMOS++ software for the analysis of biomolecular simulation trajectories. J Chem Theory Comput. 2011; 7(10): 3379–3390. Publisher Full Text
39. Reißer S, Poger D, Stroet M, et al.: Real Cost of Speed: The Effect of a Time-Saving Multiple-Time-Stepping Algorithm on the Accuracy of Molecular Dynamics Simulations. J Chem Theory Comput. 2017; 13(6): 2367–2372. PubMed Abstract | Publisher Full Text
40. Morrone JA, Zhou R, Berne BJ: Molecular Dynamics with Multiple Time Scales: How to Avoid Pitfalls. J Chem Theory Comput. 2010; 6(6): 1798–1804. PubMed Abstract | Publisher Full Text
41. Krieger E, Vriend G: New ways to boost molecular dynamics simulations. J Comput Chem. 2015; 36(13): 996–1007. PubMed Abstract | Publisher Full Text
42. Morishita T: Fluctuation formulas in molecular-dynamics simulations with the weak coupling heat bath. J Chem Phys. 2000; 113(8): 2976. Publisher Full Text
43. Sidler D, Lehner M, Frasch S, et al.: Dataset 1 in: Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations. F1000Research. 2018. http://www.doi.org/10.5256/f1000research.16715.d222949

Comments on this article Comments (0)

Version 3

VERSION 3 PUBLISHED 05 Nov 2018

Author details Author details

¹ Laboratory of Physical Chemistry, ETH Zürich, Zurich, 8093, Switzerland

Dominik Sidler
Roles: Conceptualization, Data Curation, Investigation, Methodology, Validation, Visualization, Writing – Original Draft Preparation, Writing – Review & Editing

Marc Lehner
Roles: Investigation, Methodology, Validation, Visualization

Simon Frasch
Roles: Investigation, Methodology, Validation

Michael Cristófol-Clough
Roles: Investigation, Methodology, Validation

Sereina Riniker
Roles: Conceptualization, Funding Acquisition, Methodology, Project Administration, Resources, Supervision, Writing – Review & Editing

Competing interests

No competing interests were disclosed.

Grant information

The authors gratefully acknowledge financial support by
the Swiss National Science Foundation [200021-178762]
and by Eidgenössische Technische Hochschule (ETH) Zürich [ETH-34 17-2].

Article Versions (3)

version 3

Revised

Published: 08 Mar 2019, 7:1745

https://doi.org/10.12688/f1000research.16715.3

version 2

Revised

Published: 18 Jan 2019, 7:1745

https://doi.org/10.12688/f1000research.16715.2

version 1

Published: 05 Nov 2018, 7:1745

https://doi.org/10.12688/f1000research.16715.1

© 2018 Sidler D et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Data associated with the article are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).

Download

Export To

metrics

	Views	Downloads
F1000Research	-	-
PubMed Central Data from PMC are received and updated monthly.	-	-

Citations

SEE MORE DETAILS

CITE

how to cite this article

Sidler D, Lehner M, Frasch S et al. Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations [version 1; peer review: 1 approved, 1 approved with reservations]. F1000Research 2018, 7:1745 (https://doi.org/10.12688/f1000research.16715.1)

NOTE: If applicable, it is important to ensure the information in square brackets after the title is included in all citations of this article.

track

receive updates on this article

Track an article to receive email alerts on any updates to this article.

Open Peer Review

Version 1

VERSION 1

PUBLISHED 05 Nov 2018

Views

Reviewer Report 03 Dec 2018

Nico. F. A. van der Vegt, Eduard-Zintl-Institute of Inorganic and Physical Chemistry, TU Darmstadt, Darmstadt, Germany

Approved

https://doi.org/10.5256/f1000research.18272.r40235

This is an interesting and well-conducted study on density artefacts observed in molecular dynamics simulations that use multiple time-step integration algorithms to study liquid-liquid interfacial systems. The authors nicely show how depending on choices made for pairlist updates and thermostatting algorithms twin-range schemes can lead to severe density artefacts at (water-chloroform and water-cyclohexane) liquid-liquid interfaces. It is shown that these density artefacts cannot be remedied when widely used thermostats (WC, NH, NH-chain) are used that control the velocities globally. Do the density artefacts observed in these systems disappear when a standard, stochastic thermostat is used instead?

Is the work clearly and accurately presented and does it cite the current literature?

Yes
Is the study design appropriate and is the work technically sound?

Yes
Are sufficient details of methods and analysis provided to allow replication by others?

Yes
If applicable, is the statistical analysis and its interpretation appropriate?

Yes
Are all the source data underlying the results available to ensure full reproducibility?

Yes
Are the conclusions drawn adequately supported by the results?

Yes

Competing Interests: No competing interests were disclosed.

Reviewer Expertise: molecular simulation, physical chemistry

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

CITE

Report a concern

Author Response 10 Jan 2019

Sereina Riniker, Laboratory of Physical Chemistry, ETH Zürich, Zurich, 8093, Switzerland

10 Jan 2019

Author Response

We thank the reviewer for this suggestion. We have performed additional simulations with a stochastic thermostat, which show that the artifacts disappear with this thermostat (see updated Figures 2 and ... Continue reading We thank the reviewer for this suggestion. We have performed additional simulations with a stochastic thermostat, which show that the artifacts disappear with this thermostat (see updated Figures 2 and 3). We have adapted the figures, text and conclusions accordingly.
We thank the reviewer for this suggestion. We have performed additional simulations with a stochastic thermostat, which show that the artifacts disappear with this thermostat (see updated Figures 2 and 3). We have adapted the figures, text and conclusions accordingly.
Competing Interests: No competing interests were disclosed. Close
Report a concern
Respond or Comment

COMMENTS ON THIS REPORT

Author Response 10 Jan 2019

Sereina Riniker, Laboratory of Physical Chemistry, ETH Zürich, Zurich, 8093, Switzerland

10 Jan 2019

Author Response

We thank the reviewer for this suggestion. We have performed additional simulations with a stochastic thermostat, which show that the artifacts disappear with this thermostat (see updated Figures 2 and ... Continue reading We thank the reviewer for this suggestion. We have performed additional simulations with a stochastic thermostat, which show that the artifacts disappear with this thermostat (see updated Figures 2 and 3). We have adapted the figures, text and conclusions accordingly.
We thank the reviewer for this suggestion. We have performed additional simulations with a stochastic thermostat, which show that the artifacts disappear with this thermostat (see updated Figures 2 and 3). We have adapted the figures, text and conclusions accordingly.
Competing Interests: No competing interests were disclosed. Close
Report a concern

Views

Reviewer Report 23 Nov 2018

Bojan Zagrovic, Department of Structural and Computational Biology, Max F. Perutz Laboratories, University of Vienna, Vienna, Austria

Daniel Braun, Department of Structural and Computational Biology, Max F. Perutz Laboratories, University of Vienna, Vienna, Austria

Approved with Reservations

https://doi.org/10.5256/f1000research.18272.r40231

In the article "Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations" by Sidler et al., the authors investigate the possible artefacts introduced by multiple time-step (MTS) integration methods for molecular dynamics (MD) simulations with respect to the local density and diffusivity at simulated interfaces between polar and apolar media. The simulated systems include chloroform-water and cyclohexane-water mixtures, both arranged in either planar-layer or droplet conformations. For all simulated systems, the authors report artificial density increases at the interface of the two liquids, whereby the degree of the observed effect depends on the MTS algorithm used. The density artefacts are studied mainly in the NVT ensemble, but are shown to also appear in the NPT ensemble. Moreover, some of the investigated thermostats introduce a shift in the distribution of kinetic energies in the water component. Based on these findings, the authors recommend the usage of the single-range treatment of the non-bonded interactions with lower pairlist update frequencies, which are less sensitive with respect to the discussed artefacts as compared to the widely used twin-range schemes. The study makes a valuable contribution towards the improvement of the existing MD methods when it comes to both accuracy and computational performance. The observed artefacts are certainly pronounced and cannot be ignored. On the other hand, certain aspects of the study lack the necessary explanations or depth of exploration and thus have to be improved accordingly.

Major comments

The length of individual simulations is relatively short by state-of-the-art standards. This, of course, is not a major issue if the phenomena at hand converge comparatively quickly, as indeed implied by the authors. However, for completeness, the convergence should be quantitatively demonstrated.
For all NVT simulations, the authors should perform preliminary runs of several 100 ps in the NPT ensemble for density equilibration, as well as several 100 ps of equilibration in the studied ensemble before the production runs. The overall density in the present NVT simulations is determined by the authors in the course of the setup of the simulation boxes. It is conceivable that the ideal simulated density would vary with respect to the integration method, time-step or thermostating method. As a consequence, some simulated NVT systems could exhibit too-low or too-high overall density, which on its own could introduce or at least influence some of the observed artefacts.
The microscopic origin(s) of the density artefacts should be investigated in greater detail. In particular, an in-depth investigation of the molecular structure at the interfaces would be instructive. For example, are the relative orientations of the chloroform/cyclohexane molecules in the areas of higher density different from those in the bulk? How about water? If so, how can this be related to the algorithmic differences? This analysis should, inter alia, involve a decomposition of density profiles into chloroform (or cyclohexane) and water components, but also more elaborate structural analyses (for an example of such analysis, please refer to¹).
The density profiles in Figure 2 exhibit curious short-range fluctuation patterns, which even appear somewhat periodic. The authors should analyze the origin of these patterns and potentially link them with the main analysis at hand.
The authors should analyze and comment more extensively on the asymmetric appearance of the TR-CFA-WC and TR-CFA-NH density profiles in Figure 3 with respect to the planar layers at high values of the pairlist update period. They comment that the asymmetry may be due to insufficient sampling, but given how extensive and relevant for the central argument of the paper it may be, it should be analyzed in more detail. In an ideal case (i.e. if largely symmetric) the density profiles could be averaged over the two sides for greater statistical power.

Minor comments

The CFA and RESPA approaches need to be explained in greater detail in the introduction. This seems especially important since RESPA seems to be much less prone to density artefacts than CFA.
The colors are referenced wrong in the caption of Figure 4.
The caption of Figure 3 contains an "n" without an index.
The authors should explicitly describe how the density was calculated even for planar systems (bin size, averaging volume etc.).
While the advantages of updating the Pairlist less frequently are clearly illustrated, in the conclusions section the authors should also discuss in more detail the potential downsides of such a procedure.

Is the work clearly and accurately presented and does it cite the current literature?

Yes
Is the study design appropriate and is the work technically sound?

Yes
Are sufficient details of methods and analysis provided to allow replication by others?

Partly
If applicable, is the statistical analysis and its interpretation appropriate?

Yes
Are all the source data underlying the results available to ensure full reproducibility?

Yes
Are the conclusions drawn adequately supported by the results?

Yes

References

1. Hantal G, Sega M, Kantorovich S, Schröder C, et al.: Intrinsic Structure of the Interface of Partially Miscible Fluids: An Application to Ionic Liquids. The Journal of Physical Chemistry C. 2015; 119 (51): 28448-28461 Publisher Full Text

Competing Interests: No competing interests were disclosed.

We confirm that we have read this submission and believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however we have significant reservations, as outlined above.

CITE

Report a concern

Author Response 10 Jan 2019

Sereina Riniker, Laboratory of Physical Chemistry, ETH Zürich, Zurich, 8093, Switzerland

10 Jan 2019

Author Response
Major comments:
1. We have checked the convergence, which is reached within 100 ps (new Figure S6 in the Supporting Information), well below the simulation time.
2. The individual
... Continue reading
Major comments:

We have checked the convergence, which is reached within 100 ps (new Figure S6 in the Supporting Information), well below the simulation time.

The individual solvent boxes were equilibrated in the NPT ensemble, before combining them in the layer system. In the absence of artefacts (e.g. with a single-range scheme), the density profiles remain perfectly constant in the NVT ensemble. We consider this to be a more suited starting configuration than the layer system equilibrated additionally in the NPT ensemble with the twin-range scheme, because it shows the full extent of the artefacts introduced by the latter.

We thank the reviewers for this suggestion. We have investigated the orientation of the solvent molecules at the interface, which is affected by the artefacts from the MTS algorithm. The results are presented in the new Figure S5 in the Supporting Information and discussed in the text. A decomposition of the density profiles is shown in the new Figure S2.

These "fluctuations" stem from the binning of the density profiles. We have changed Figure 2 to have a consistent bin size for all curves.

These fluctuations arise mainly due to the finite box size. We have adapted the text accordingly.

Minor comments:

We have added additional explanations of the MTS algorithms with equations.

We thank the reviewer for spotting this. We have corrected the caption of Figure 4 accordingly.

We thank the reviewer for spotting this. We have corrected the caption of Figure 3 accordingly.

We have added these details in the Methods section.

Due to the new results with the stochastic dynamics thermostat, the conclusions have been adapted.
Major comments:

We have checked the convergence, which is reached within 100 ps (new Figure S6 in the Supporting Information), well below the simulation time.

The individual solvent boxes were equilibrated in the NPT ensemble, before combining them in the layer system. In the absence of artefacts (e.g. with a single-range scheme), the density profiles remain perfectly constant in the NVT ensemble. We consider this to be a more suited starting configuration than the layer system equilibrated additionally in the NPT ensemble with the twin-range scheme, because it shows the full extent of the artefacts introduced by the latter.

We thank the reviewers for this suggestion. We have investigated the orientation of the solvent molecules at the interface, which is affected by the artefacts from the MTS algorithm. The results are presented in the new Figure S5 in the Supporting Information and discussed in the text. A decomposition of the density profiles is shown in the new Figure S2.

These "fluctuations" stem from the binning of the density profiles. We have changed Figure 2 to have a consistent bin size for all curves.

These fluctuations arise mainly due to the finite box size. We have adapted the text accordingly.

Minor comments:

We have added additional explanations of the MTS algorithms with equations.

We thank the reviewer for spotting this. We have corrected the caption of Figure 4 accordingly.

We thank the reviewer for spotting this. We have corrected the caption of Figure 3 accordingly.

We have added these details in the Methods section.

Due to the new results with the stochastic dynamics thermostat, the conclusions have been adapted.
Competing Interests: No competing interests were disclosed. Close
Report a concern
Respond or Comment

COMMENTS ON THIS REPORT

Author Response 10 Jan 2019

Sereina Riniker, Laboratory of Physical Chemistry, ETH Zürich, Zurich, 8093, Switzerland

10 Jan 2019

Author Response
Major comments:
1. We have checked the convergence, which is reached within 100 ps (new Figure S6 in the Supporting Information), well below the simulation time.
2. The individual
... Continue reading
Major comments:

We have checked the convergence, which is reached within 100 ps (new Figure S6 in the Supporting Information), well below the simulation time.

The individual solvent boxes were equilibrated in the NPT ensemble, before combining them in the layer system. In the absence of artefacts (e.g. with a single-range scheme), the density profiles remain perfectly constant in the NVT ensemble. We consider this to be a more suited starting configuration than the layer system equilibrated additionally in the NPT ensemble with the twin-range scheme, because it shows the full extent of the artefacts introduced by the latter.

We thank the reviewers for this suggestion. We have investigated the orientation of the solvent molecules at the interface, which is affected by the artefacts from the MTS algorithm. The results are presented in the new Figure S5 in the Supporting Information and discussed in the text. A decomposition of the density profiles is shown in the new Figure S2.

These "fluctuations" stem from the binning of the density profiles. We have changed Figure 2 to have a consistent bin size for all curves.

These fluctuations arise mainly due to the finite box size. We have adapted the text accordingly.

Minor comments:

We have added additional explanations of the MTS algorithms with equations.

We thank the reviewer for spotting this. We have corrected the caption of Figure 4 accordingly.

We thank the reviewer for spotting this. We have corrected the caption of Figure 3 accordingly.

We have added these details in the Methods section.

Due to the new results with the stochastic dynamics thermostat, the conclusions have been adapted.
Major comments:

We have checked the convergence, which is reached within 100 ps (new Figure S6 in the Supporting Information), well below the simulation time.

The individual solvent boxes were equilibrated in the NPT ensemble, before combining them in the layer system. In the absence of artefacts (e.g. with a single-range scheme), the density profiles remain perfectly constant in the NVT ensemble. We consider this to be a more suited starting configuration than the layer system equilibrated additionally in the NPT ensemble with the twin-range scheme, because it shows the full extent of the artefacts introduced by the latter.

We thank the reviewers for this suggestion. We have investigated the orientation of the solvent molecules at the interface, which is affected by the artefacts from the MTS algorithm. The results are presented in the new Figure S5 in the Supporting Information and discussed in the text. A decomposition of the density profiles is shown in the new Figure S2.

These "fluctuations" stem from the binning of the density profiles. We have changed Figure 2 to have a consistent bin size for all curves.

These fluctuations arise mainly due to the finite box size. We have adapted the text accordingly.

Minor comments:

We have added additional explanations of the MTS algorithms with equations.

We thank the reviewer for spotting this. We have corrected the caption of Figure 4 accordingly.

We thank the reviewer for spotting this. We have corrected the caption of Figure 3 accordingly.

We have added these details in the Methods section.

Due to the new results with the stochastic dynamics thermostat, the conclusions have been adapted.
Competing Interests: No competing interests were disclosed. Close
Report a concern

Comments on this article Comments (0)

Version 3

VERSION 3 PUBLISHED 05 Nov 2018

Open Peer Review

Reviewer Status

Reviewer Reports

	Invited Reviewers
	1	2
Version 3 (revision) 08 Mar 19	read
Version 2 (revision) 18 Jan 19	read
Version 1 05 Nov 18	read	read

Bojan Zagrovic, University of Vienna, Vienna, Austria

Daniel Braun, University of Vienna, Vienna, Austria
Nico. F. A. van der Vegt, TU Darmstadt, Darmstadt, Germany

Comments on this article

All Comments(0)

Add a comment

Browse by related subjects

Back to all reports

Reviewer Report

13 Views

22 Mar 2019 | for Version 3

Bojan Zagrovic, Department of Structural and Computational Biology, Max F. Perutz Laboratories, University of Vienna, Vienna, Austria

Daniel Braun, Department of Structural and Computational Biology, Max F. Perutz Laboratories, University of Vienna, Vienna, Austria

13 Views Cite this report Responses(0)

Approved

The authors have adequately addressed our comments.

Competing Interests

No competing interests were disclosed.

We confirm that we have read this submission and believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

Respond to this report

Responses (0)

Back to all reports

Reviewer Report

28 Views

29 Jan 2019 | for Version 2

Bojan Zagrovic, Department of Structural and Computational Biology, Max F. Perutz Laboratories, University of Vienna, Vienna, Austria

Daniel Braun, Department of Structural and Computational Biology, Max F. Perutz Laboratories, University of Vienna, Vienna, Austria

28 Views Cite this report Responses(1)

Approved With Reservations

The authors have adequately addressed most of the concerns raised in the original review. However, there are still some open questions that need to be clarified before the manuscript is indexed. The numbers below refer to the numbering in the original review.

Major comments:

Comment 1. It is not fully clear what the chloroform density in Figure S6 is and how it is calculated. Assuming chloroform density describes the density of the chloroform component, how was the respective volume calculated? Does the plot really show data for each frame, as stated, or is the data averaged in some way? Secondly, our original concern related to the convergence of the phenomena at hand, i.e. the buildup of the density artefacts. In the main text, it is stated that a density increase at the interface "was accompanied by a density decrease further away from the interface due to the constant volume conditions" and not by a change in density of the whole component (as also shown in Figure S2). Therefore, it is not clear how the quantity given in Figure S6 relates to density artefacts.

Comment 2. Regarding the box size in constant volume simulations, we believe that there is an issue that needs to be analyzed in more detail. In general, it is common practice to use constant pressure simulations to determine the appropriate box size for constant volume simulations after finalizing the MD setup, including the definition of integrator, thermostat, pairlist specifications, etc. We acknowledge that in a study which aims at comparing different MD setups, the researcher faces a dilemma. On the one hand, equilibrating the box size for each MD setup entails a variation in the overall density. On the other hand, comparing all MD setups at the same box size entails a variation in pressure, with a potentially even more severe influence on a number of properties, e.g. system dynamics and, as relevant here, also on the extent of the observed artefacts in the present study. The authors should address this topic in more detail. In the Methods section, this should include a detailed explanation of the problem at hand and the choices that the authors have made (variation in pressure vs. variation in density). In particular, a list of the overall values of density/pressure for each system should be included. Furthermore, the authors should provide an estimate of the influence of their choice of how to set up the NVT simulations on system properties and foremost on the density artefacts in their paper.

Comment 4. Figure 2 still shows somewhat systematically conserved fluctuations (i.e. present across different simulation systems). This was even more obvious in the original figure. Independent of bin size, where do these (apparently) conserved fluctuations come from?

Minor comments:

Comment 5. The downsides of SR increased pairlist update intervals should be discussed.

Competing Interests

No competing interests were disclosed.

Respond to this report

Responses (1)

Author Response

04 Mar 2019

Sereina Riniker, Laboratory of Physical Chemistry, ETH Zürich, Zurich, 8093, Switzerland

Major comments:
(1) We have adapted the caption of Figure S6 and added Equation (1), which was used in the calculation, in the Supplementary Material.
(2) As the density artefacts were observed in both NVT and NPT simulations, we can rule out this factor as a cause for the occurrence of the artefacts. Nevertheless, we have added a sentence in the Methods section describing the issue and our rationale for choosing the equilibration procedure. In addition, we have added Table S1 in the Supplementary Material with the pressure and box volume in the NVT and NPT simulations of the water-cyclohexane system. The data is discussed in the Section "Chemical composition and constant pressure".
(4) The apparent fluctuations are due to a combination of the relative orientation of the molecules at the interface (see Figure S5), the granularity of the atoms in a classical force field, and the binning.

Minor comment:
(5) We have added a short discussion of this issue in the Conclusions.

View more View less

Competing Interests

No competing interests were disclosed.

Back to all reports

Reviewer Report

30 Views

03 Dec 2018 | for Version 1

Nico. F. A. van der Vegt, Eduard-Zintl-Institute of Inorganic and Physical Chemistry, TU Darmstadt, Darmstadt, Germany

30 Views Cite this report Responses(1)

Approved

Is the work clearly and accurately presented and does it cite the current literature?

Yes
Is the study design appropriate and is the work technically sound?

Yes
Are sufficient details of methods and analysis provided to allow replication by others?

Yes
If applicable, is the statistical analysis and its interpretation appropriate?

Yes
Are all the source data underlying the results available to ensure full reproducibility?

Yes
Are the conclusions drawn adequately supported by the results?

Yes

Competing Interests

No competing interests were disclosed.

Reviewer Expertise

molecular simulation, physical chemistry

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

Respond to this report

Responses (1)

Back to all reports

Reviewer Report

43 Views

23 Nov 2018 | for Version 1

Bojan Zagrovic, Department of Structural and Computational Biology, Max F. Perutz Laboratories, University of Vienna, Vienna, Austria

Daniel Braun, Department of Structural and Computational Biology, Max F. Perutz Laboratories, University of Vienna, Vienna, Austria

43 Views Cite this report Responses(1)

Approved With Reservations

The length of individual simulations is relatively short by state-of-the-art standards. This, of course, is not a major issue if the phenomena at hand converge comparatively quickly, as indeed implied by the authors. However, for completeness, the convergence should be quantitatively demonstrated.
For all NVT simulations, the authors should perform preliminary runs of several 100 ps in the NPT ensemble for density equilibration, as well as several 100 ps of equilibration in the studied ensemble before the production runs. The overall density in the present NVT simulations is determined by the authors in the course of the setup of the simulation boxes. It is conceivable that the ideal simulated density would vary with respect to the integration method, time-step or thermostating method. As a consequence, some simulated NVT systems could exhibit too-low or too-high overall density, which on its own could introduce or at least influence some of the observed artefacts.
The microscopic origin(s) of the density artefacts should be investigated in greater detail. In particular, an in-depth investigation of the molecular structure at the interfaces would be instructive. For example, are the relative orientations of the chloroform/cyclohexane molecules in the areas of higher density different from those in the bulk? How about water? If so, how can this be related to the algorithmic differences? This analysis should, inter alia, involve a decomposition of density profiles into chloroform (or cyclohexane) and water components, but also more elaborate structural analyses (for an example of such analysis, please refer to¹).
The density profiles in Figure 2 exhibit curious short-range fluctuation patterns, which even appear somewhat periodic. The authors should analyze the origin of these patterns and potentially link them with the main analysis at hand.
The authors should analyze and comment more extensively on the asymmetric appearance of the TR-CFA-WC and TR-CFA-NH density profiles in Figure 3 with respect to the planar layers at high values of the pairlist update period. They comment that the asymmetry may be due to insufficient sampling, but given how extensive and relevant for the central argument of the paper it may be, it should be analyzed in more detail. In an ideal case (i.e. if largely symmetric) the density profiles could be averaged over the two sides for greater statistical power.

Minor comments

The CFA and RESPA approaches need to be explained in greater detail in the introduction. This seems especially important since RESPA seems to be much less prone to density artefacts than CFA.
The colors are referenced wrong in the caption of Figure 4.
The caption of Figure 3 contains an "n" without an index.
The authors should explicitly describe how the density was calculated even for planar systems (bin size, averaging volume etc.).
While the advantages of updating the Pairlist less frequently are clearly illustrated, in the conclusions section the authors should also discuss in more detail the potential downsides of such a procedure.

Is the work clearly and accurately presented and does it cite the current literature?

Yes
Is the study design appropriate and is the work technically sound?

Yes
Are sufficient details of methods and analysis provided to allow replication by others?

Partly
If applicable, is the statistical analysis and its interpretation appropriate?

Yes
Are all the source data underlying the results available to ensure full reproducibility?

Yes
Are the conclusions drawn adequately supported by the results?

Yes

References

Competing Interests

No competing interests were disclosed.

Respond to this report

Responses (1)

Author Response

10 Jan 2019

Sereina Riniker, Laboratory of Physical Chemistry, ETH Zürich, Zurich, 8093, Switzerland

Major comments:

We have checked the convergence, which is reached within 100 ps (new Figure S6 in the Supporting Information), well below the simulation time.
The individual solvent boxes were equilibrated in the NPT ensemble, before combining them in the layer system. In the absence of artefacts (e.g. with a single-range scheme), the density profiles remain perfectly constant in the NVT ensemble. We consider this to be a more suited starting configuration than the layer system equilibrated additionally in the NPT ensemble with the twin-range scheme, because it shows the full extent of the artefacts introduced by the latter.
We thank the reviewers for this suggestion. We have investigated the orientation of the solvent molecules at the interface, which is affected by the artefacts from the MTS algorithm. The results are presented in the new Figure S5 in the Supporting Information and discussed in the text. A decomposition of the density profiles is shown in the new Figure S2.
These "fluctuations" stem from the binning of the density profiles. We have changed Figure 2 to have a consistent bin size for all curves.
These fluctuations arise mainly due to the finite box size. We have adapted the text accordingly.

Minor comments:

We have added additional explanations of the MTS algorithms with equations.
We thank the reviewer for spotting this. We have corrected the caption of Figure 4 accordingly.
We thank the reviewer for spotting this. We have corrected the caption of Figure 3 accordingly.
We have added these details in the Methods section.
Due to the new results with the stochastic dynamics thermostat, the conclusions have been adapted.

View more View less

Competing Interests

No competing interests were disclosed.

Alongside their report, reviewers assign a status to the article:

Approved - the paper is scientifically sound in its current form and only minor, if any, improvements are suggested

Approved with reservations - A number of small changes, sometimes more significant revisions are required to address specific details and improve the papers academic merit.

Not approved - fundamental flaws in the paper seriously undermine the findings and conclusions

Click here to access the data.

Downloaded data do not display as expected? Download the data

[1] 1. Bishop TC, Skeel RD, Schulten K: Difficulties with multiple time stepping and fast multipole algorithm in molecular dynamics. J Comput Chem. 1997; 18(14): 1785–1791. Publisher Full Text

[2] 2. Ma Q, Izaguirre JA, Skeel RD: Verlet-I/r-RESPA/Impulse is limited by nonlinear instabilities. SIAM J Sci Comput. 2003; 24(6): 1951–1973. Publisher Full Text

[3] 3. Ryckaert JP, Ciccotti G, Berendsen HJ: Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J Comput Phys. 1977; 23(3): 327–341. Publisher Full Text

[4] 4. Miyamoto S, Kollman PA: SETTLE: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. J Comput Chem. 1992; 13(8): 952–962. Publisher Full Text

[5] 5. Kräutler V, van Gunsteren WF, Hünenberger PH: A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J Comput Chem. 2001; 22(5): 501–508. Publisher Full Text

[6] 6. Hess B, Bekker H, Berendsen HJC, et al.: LINCS: A linear constraint solver for molecular simulations. J Comput Chem. 1997; 18(12): 1463–1472. Publisher Full Text

[7] 7. Christen M, van Gunsteren WF: An approximate but fast method to impose flexible distance constraints in molecular dynamics simulations. J Chem Phys. 2005; 122(14): 144106. PubMed Abstract | Publisher Full Text

[8] 8. Berendsen HJ, Van Gunsteren WF, Zwinderman HR, et al.: Simulations of proteins in water. Ann N Y Acad Sci. 1986; 482(1): 269–286. PubMed Abstract | Publisher Full Text

[9] 9. Van Gunsteren WF, Berendsen HJ, Geurtsen RG, et al.: A molecular dynamics computer simulation of an eight-base-pair DNA fragment in aqueous solution: comparison with experimental two-dimensional NMR data. Ann N Y Acad Sci. 1986; 482(1): 287–303. PubMed Abstract | Publisher Full Text

[10] 10. van Gunsteren WF, Berendsen HJ: Computer simulation of molecular dynamics: Methodology, applications, and perspectives in chemistry. Angew Chem Int Ed. 1990; 29(9): 992–1023. Publisher Full Text

[11] 11. Tironi IG, Sperb R, Smith PE, et al.: A generalized reaction field method for molecular dynamics simulations. J Chem Phys. 1995; 102: 5451–5459. Publisher Full Text

[12] 12. Sidler D, Frasch S, Cristòfol-Clough M, et al.: Anisotropic reaction field correction for long-range electrostatic interactions in molecular dynamics simulations. J Chem Phys. 2018; 148(23): 234105. PubMed Abstract | Publisher Full Text

[13] 13. Hockney RW, Eastwood JW: Computer Simulation Using Particles. McGraw-Hill, New York, NY, 1981. Reference Source

[14] 14. Darden T, York D, Pedersen L: Particle mesh Ewald: An N log (N) method for Ewald sums in large systems. J Chem Phys. 1993; 98(12): 10089–10092. Publisher Full Text

[15] 15. Oostenbrink C, Villa A, Mark AE, et al.: A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. J Comput Chem. 2004; 25(13): 1656–1676. PubMed Abstract | Publisher Full Text

[16] 16. Wang J, Cieplak P, Kollman PA: How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J Comput Chem. 2000; 21(12): 1049–1074. Publisher Full Text

[17] 17. Tuckerman M, Berne BJ, Martyna GJ: Reversible multiple time scale molecular dynamics. J Chem Phys. 1992; 97: 1990–2001. Publisher Full Text

[18] 18. Hockney RW: Potential calculation and some applications. Technical report, Langley Research Center, Hampton, Va., 1970. Reference Source

[19] 19. Tuckerman M, Berne BJ, Martyna GJ: Reply to comment on: Reversible multiple time scale molecular dynamics. J Chem Phys. 1993; 99: 2278–2279. Publisher Full Text

[20] 20. Swope WC, Andersen HC, Berens PH, et al.: A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. J Chem Phys. 1982; 76: 637–649. Publisher Full Text

[21] 21. Minary P, Tuckerman M, Martyna GJ: Long time molecular dynamics for enhanced conformational sampling in biomolecular systems. Phys Rev Lett. 2004; 93(15): 150201. PubMed Abstract | Publisher Full Text

[22] 22. Leimkuhler B, Margul DT, Tuckerman ME: Stochastic, resonance-free multiple time-step algorithm for molecular dynamics with very large time steps. Mol Phys. 2013; 111(22–23): 3579–3594. Publisher Full Text

[23] 23. Margul DT, Tuckerman ME: A Stochastic, Resonance-Free Multiple Time-Step Algorithm for Polarizable Models That Permits Very Large Time Steps. J Chem Theory Comput. 2016; 12(5): 2170–2180. PubMed Abstract | Publisher Full Text

[24] 24. Morrone JA, Markland TE, Ceriotti M, et al.: Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances. J Chem Phys. 2011; 134(1): 014103. PubMed Abstract | Publisher Full Text

[25] 25. Berendsen HJC, Postma JPM, van Gunsteren WF, et al.: Molecular dynamics with coupling to an external bath. J Chem Phys. 1984; 81(8): 3684–3690. Publisher Full Text

[26] 26. Nosé S: A unified formulation of the constant temperature molecular dynamics methods. J Chem Phys. 1984; 81(1): 511–519. Publisher Full Text

[27] 27. Hoover WG: Canonical dynamics: Equilibrium phase-space distributions. Phys Rev A Gen Phys. 1985; 31(3): 1695–1697. PubMed Abstract | Publisher Full Text

[28] 28. Martyna GJ, Klein ML, Tuckerman M: Nosé–hoover chains: the canonical ensemble via continuous dynamics. J Chem Phys. 1992; 97(4): 2635–2643. Publisher Full Text

[29] 29. Schneider T, Stoll E: Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions. Phys Rev B. 1978; 17(3): 1302. Publisher Full Text

[30] 30. Andersen HC: Molecular dynamics simulations at constant pressure and/or temperature. J Chem Phys. 1980; 72(4): 2384–2393. Publisher Full Text

[31] 31. Bussi G, Donadio D, Parrinello M: Canonical sampling through velocity rescaling. J Chem Phys. 2007; 126(1): 014101. PubMed Abstract | Publisher Full Text

[32] 32. Schlick T, Mandziuk M, Skeel RD, et al.: Nonlinear resonance artifacts in molecular dynamics simulations. J Chem Phys. 1998; 140(1): 1–29. Publisher Full Text

[33] 33. Sandu A, Schlick T: Masking resonance artifacts in force-splitting methods for biomolecular simulations by extrapolative Langevin dynamics. J Chem Phys. 1999; 151(1): 74–113. Publisher Full Text

[34] 34. Schmid N, Christ CD, Christen M, et al.: Architecture, implementation and parallelisation of the GROMOS software for biomolecular simulation. Comput Phys Commun. 2012; 183(4): 890–903. Publisher Full Text

[35] 35. Haynes WM: CRC Handbook of Chemistry and Physics. CRC Press, Boca Raton, FL, 2014. Reference Source

[36] 36. Tironi IG, van Gunsteren WF: A molecular dynamics simulation study of chloroform. Mol Phys. 1994; 83(2): 381–403. Publisher Full Text

[37] 37. Berendsen HJC, Postma JPM, van Gunsteren WF, et al.: Interaction models for water in relation to protein hydration. In Intermolecular forces. Springer, 1981; 331–342. Publisher Full Text

[38] 38. Eichenberger AP, Allison JR, Dolenc J, et al.: GROMOS++ software for the analysis of biomolecular simulation trajectories. J Chem Theory Comput. 2011; 7(10): 3379–3390. Publisher Full Text

[39] 39. Reißer S, Poger D, Stroet M, et al.: Real Cost of Speed: The Effect of a Time-Saving Multiple-Time-Stepping Algorithm on the Accuracy of Molecular Dynamics Simulations. J Chem Theory Comput. 2017; 13(6): 2367–2372. PubMed Abstract | Publisher Full Text

[40] 40. Morrone JA, Zhou R, Berne BJ: Molecular Dynamics with Multiple Time Scales: How to Avoid Pitfalls. J Chem Theory Comput. 2010; 6(6): 1798–1804. PubMed Abstract | Publisher Full Text

[41] 41. Krieger E, Vriend G: New ways to boost molecular dynamics simulations. J Comput Chem. 2015; 36(13): 996–1007. PubMed Abstract | Publisher Full Text

[42] 42. Morishita T: Fluctuation formulas in molecular-dynamics simulations with the weak coupling heat bath. J Chem Phys. 2000; 113(8): 2976. Publisher Full Text

[43] 43. Sidler D, Lehner M, Frasch S, et al.: Dataset 1 in: Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations. F1000Research. 2018. http://www.doi.org/10.5256/f1000research.16715.d222949

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

Abstract

Keywords

Introduction

Figure 1. Simulated Systems.

Methods

Systems

Table 1. Overview of the system dimensions and simulation times.

Analysis

Results & discussion

Figure 2. Profile of the density ρ in the water-chloroform layer system with respect to the layer normal z.

Figure 3. Density profile with respect to the layer normal z and the pairlist update period np∆t.

Figure 4. Spacial density difference profiles with respect to a SR reference simulation (nm = np = 1) under NVT conditions for layer systems (left) and droplet systems (right) with water-chloroform composition (top) or water-cyclohexane composition (bottom).

Table 2. Comparison of run times for three realistic simulation setups with different nonbonded force treatment: single-range scheme (SR) and twin-range scheme with constant-force application (TR CFA).

Chemical composition and constant pressure

Figure 5. Density profile of the water-cyclohexane layer system under NVT conditions (left) or NPT conditions (right).

Geometry of the interface

Thermostatting method

Figure 6. Histogram of the kinetic energy for chloroform (left) and water (right) using a weak coupling (WC, black), Nosé-Hoover (NH, red) and Nosé-Hoover chain (NHchain, green) thermostat (bin size 1 K).

Conclusions

Data availability

Grant information

Acknowledgements

Supplementary material

References

Comments on this article Comments (0)

Open Peer Review

Comments on this article Comments (0)

Open Peer Review

Reviewer Status

Reviewer Reports

Comments on this article

Browse by related subjects

The problem

How to fix it

Competing Interests Policy

Stay Updated

Figure 3. Density profile with respect to the layer normal z and the pairlist update period n_p∆t.

Figure 4. Spacial density difference profiles with respect to a SR reference simulation (n_m = n_p = 1) under NVT conditions for layer systems (left) and droplet systems (right) with water-chloroform composition (top) or water-cyclohexane composition (bottom).