ALL Metrics
-
Views
-
Downloads
Get PDF
Get XML
Cite
Export
Track
Research Article
Revised

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

[version 3; peer review: 2 approved]
PUBLISHED 08 Mar 2019
Author details Author details
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 for weak coupling and Nosé-Hoover (chain) thermostats. 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 or a stochastic dynamics thermostat, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs.
Conclusion: TR schemes are not advisable to use in combination with weak coupling or Nosé-Hoover (chain) thermostats due to the occurrence of significant numerical artifacts at interfaces.

Keywords

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

Revised Amendments from Version 2

We have added short discussions of the issue of NVT versus NPT simulations (including Table S1 in Supplementary File 1), and the limit of lowering the pairlist-update frequency in the SR scheme, as well as a more detailed explanation of Figure S6 in Supplementary File 1.

See the authors' detailed response to the review by Bojan Zagrovic and Daniel Braun
See the authors' detailed response to the review by Nico. F. A. van der Vegt

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 subproblems1,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 SHAKE3, SETTLE4, M-SHAKE5, LINCS6, or FLEXSHAKE7. Alternatively, multiple time-step (MTS) integration is an efficient approach to accommodate motions on different time scales. In the twin-range (TR) scheme810, the pairwise nonbonded forces finb acting on particle i are split into a short-range contribution from particles j within the cutoff Rs , and a mid-range contribution from particles j between Rs and the long-range cutoff Rl (see Figure 1),

00e60184-e594-4462-8011-c25cc7328b0c_figure1.gif

Figure 1. Simulated Systems.

Left: Schematic illustration of the twin-range (TR) setup with the inner cutoff Rs and the outer cutoff Rl. Interactions beyond Rl 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.

finb(t)=finb,s(t)+finb,m(t).(1)

with

finb,s(t)=jN(rijRs)fijele+fijvdw,(2)

finb,m(t)=jN(Rs<rijRl)fijele+fijvdw(3)

The short-range forces finb,s are updated every time step ∆t (inner time step), while the mid-range forces finb,m are updated every nm-th time step, nmt with nm ≥ 1 (outer time step). Long-range electrostatic interactions beyond the cutoff radius Rl can be included using the (an)isotropic reaction-field methods (RF or ARF)11,12, or lattice-sum methods13,14. Although theoretically possible, the TR scheme is generally not used together with lattice-sum methods. Simulations with (A)RF typically use larger cutoffs Rl than simulations with PME (e.g. 1.4 nm for the GROMOS force field15 compared to 0.9 nm for the AMBER force field16), 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 np∆t should be equal to the period of updating the mid-range forces nmt. In the following, the TR scheme is defined by Rs < Rl and nm = np > 1, and the single-range (SR) scheme by Rs = Rl , nm = 1 and np ≥ 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) finb,m(CFA)(t) at every inner time step ∆t, or (ii) using a reversible reference system propagator algorithm (RESPA)810, which is based on the Trotter factorisation of the Liouville propagator17. In the case of RESPA, scaled mid-range forces finb,m(RESPA)(t) are typically applied impulse-wise, i.e. only at every outer time step nmt. Note that in this work the leap-frog scheme18 is used for the integration of the equations of motion, which is similar19 to velocity-Verlet20 as typically employed with RESPA. The resulting discrete integration scheme can be written as follows,

xi(t+Δt)=xi(t)+vi(t+12Δt)Δt,vi(t+12Δt)=vi(t12Δt)+1mi(finb,s(t)+finb,m()(t))Δt,finb,m(CFA)(t):=finb,m(ttmodnmΔt),finb,m(RESPA)(t):={nmfinb,m(t),tmodnmΔt=00,else.

Note that for simplicity, the notation is restricted to the microcanonical case and it is implicitly assumed that t mod ∆t: = ∆t(n mod 1) = 0, i.e. the time t = nt is discretized with respect to ∆t, where mod denotes the modulo operation.

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 freedom2123, or by selectively targeting the high-frequency modes using a colored noise thermostat24. In practice, however, standard thermostats like weak-coupling25, Nosé-Hoover (chain)2628, or stochastic thermostats2931 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 systems32,33.

Methods

For all simulations, the GROMOS software package34 version 1.4.0 was used with the GROMOS 53A6 force field15. Periodic boundary conditions were applied. Newton’s equations of motion were integrated using the leap-frog scheme18 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. nm , np ∈ [1, 20]). If not indicated otherwise, the SR scheme corresponds to a midrange force evaluation every 2 fs (nm = 1) and a pairlist update every 10 fs (np = 5). The TR scheme typically corresponds to a mid-range force evaluation in combination with a pairlist update every 10 fs (nm = np = 5). If not stated otherwise, the CFA method was used for the mid-range force contributions. The inner cutoff radius Rs was set to 0.8 nm and the outer cutoff radius Rl to 1.4 nm. 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−4. Four different thermostats were tested, i.e. weak coupling25 (WC), Nosé-Hoover26,27 (NH), Nosé-Hoover chain28 (NHchain), and stochastic dynamics (SD) i.e. a Langevin thermostat35. Different temperature baths were used for different solvents. The coupling time τ was set to 0.1 ps for all non-stochastic thermostatting methods. A chain length of three was employed for the NHchain thermostat. For SD, a friction coefficient γ = 10 ps-1 was used. A standard homogeneous RF approach11 with a charge-group based pairlist algorithm was used (each solvent molecule is a single charge group) to mimic long-range electrostatic interactions beyond the cutoff Rl. The dielectric permittivity was set to 78.5, which corresponds to the experimentally measured value of water36.

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.

SetupPhase 1 [#molecules]Phase 2 [#molecules]Box (x, y, z) [nm]Time [ns]
Layer2002 CHCl38129 H2O(4.1, 4.1, 30.2)4
Layer6750 C6H1239’366 H2O(8.5, 8.5, 33.9)1
Droplet847 CHCl359’484 H2O(12.2, 12.2, 12.2)1
Droplet619 C6H1259’018 H2O(12.2, 12.2, 12.2)1

The main system consisted of 2002 chloroform (CHCl3)37 molecules and 8129 SPC38 water molecules within an elongated rectangular box of 30.2 nm length in z-direction. The initial coordinates were generated by combining pre-equilibrated (under NPT conditions) boxes of chloroform and water. The combined system was equilibrated for 1 ns using the SR scheme under NVT conditions. Although this meant a slightly too small box size for the NVT conditions, we chose this procedure to avoid effects from the MTS scheme and/or the barostatting method. The simulation length of the production simulations was 1–4 ns, performed under NVT conditions. For the water-cyclohexane layer system, 6750 cyclohexane (C6H12) and 39’366 water molecules were used in an elongated rectangular box of 33.85 nm length in z-direction. Coordinates were generated by combining pre-equilibrated boxes of cyclohexane and water. This system was simulated under NVT as well as NPT conditions (using a Berendsen barostat25 with a coupling time τb = 0.5 ps). The same starting coordinates were used for all production runs to allow direct comparison, followed by 200 ps of equilibration under the given conditions.

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 CHCl3 molecules in 59’484 SPC water molecules, or a droplet with 619 C6H12 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 programs39. For the planar systems, the densities were obtained from 100 equidistant bins along the z-axis, e.g. for the main chloroform-water system the averaging bin volume was 5.0 nm3. 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 quadratically 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 nm = np = 5 led to a significant density increase in chloroform close to the interface for WC, NH and NHchain thermostats. A corresponding graph with the separate densities of the two phases is given in Figure S2, Supplementary File 1. The small density fluctuations in Figure 2 are due to the relative orientation of the molecules at the interface in combination with the granularity of the molecules and the binning. On the contrary, no significant density increase was found for SD for pairlist-update periods up to 40 fs as shown in Figure 3. However, SD appears to cause slightly more pronounced local density fluctuations in the chloroform phase compared to the non-stochastic thermostats for the chosen simulation time of 1 ns, independent of the pairlist-update period.

00e60184-e594-4462-8011-c25cc7328b0c_figure2.gif

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

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

00e60184-e594-4462-8011-c25cc7328b0c_figure3.gif

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

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

For the non-stochastic thermostats, the density increase of the TR setup was accompanied by a density decrease further away from the interface due to the constant volume conditions. An analysis of the density resonance pattern with respect to the outer time step nmt is shown in Figure 3. Note that the density artefacts are not perfectly symmetric with respect to the layer normal due to NVT conditions and finite layer thickness. The asymmetry for large outer time steps is due to the small thickness of the layers, which enables aggregation at one interface and depletion at the other interface. It can be seen that the TR CFA scheme leads to large resonance artefacts beyond nmt ≈ 20 fs. On the other hand, evaluating the mid-range forces every step while keeping the pairlist update every fifth step (i.e. SR scheme with nm = 1 and np = 5) resolved the issue almost entirely for all investigated thermostats (Figure 2 and Figure 3, data not shown for 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 (nm = np = 1).

00e60184-e594-4462-8011-c25cc7328b0c_figure4.gif

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

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

These findings indicate that the observed density increase among non-stochastic thermostats is mainly caused by resonance artefacts due to the CFA MTS integration. This argumentation is supported by the pairlist error analysis shown in Figure S3, 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 nm. 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 S4, Supplementary File 1). In addition, it could be shown that the erroneous density at the interface leads to a slight deviation of the molecular orientations (Figure S5, Supplementary File 1). Note that the build-up of the density artefacts for realistic MTS setups evolves on a time scale much shorter than the duration of performed simulations (Figure S6, 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 nmt = 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 of non-stochastic thermostats, 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 implementations40, also for impulse-wise MTS algorithms1,2,41.

On the contrary, only small effects on the density profile could be observed for the SR scheme with pairlist update periods up to npt = 40 fs (Figure 3 and Figure 4). Due to the large cutoff Rl used here (i.e. 1.4 nm), the changes in the pairlist are small between updates and scale sub-linearly with nm (Figure S3, Supplementary File 1). This observation is in line with the study of Krieger et al.42, where they found a negligible energy drift when using a similarly low pair-list update frequency in a homogeneous system. When insisting on a non-stochastic thermostat, 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 partially compensated by a less frequent update of the pairlist (Table 2).

Table 2. Comparison of run times for three simulation setups with different nonbonded force treatment: single-range scheme (SR) and twin-range scheme with constant-force application (TR CFA) using a Nosé-Hoover (NH) thermostat.

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 the standard charge-group based algorithm of the GROMOS program. The simulated system consists of 2002 chloroform and 8129 water molecules.

Force splittingnmnpRun time [h]
TR CFA556.1
SR1511.4
SR1208.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). The average pressure in the NVT simulations and the average volume in the NPT simulations are listed in Table S1 (Supplementary File 1). Note that under NVT conditions, the observed density artefacts for the TR scheme lead to a reduction of the pressure as a function of the pairlist update frequency. As mentioned the box volume is slightly too small due to the chosen setup. Under NPT conditions and using the SR scheme, the box volume remains relatively constant. However, the semi-anisotropic pressure coupling in the x,y-direction with fixed z-dimension did not allow to preserve the system pressure at 1 atm for the elongated bilayer system. Separate scaling of the z-dimension did not improve the pressure control of the given system (data not shown). Using the TR scheme, the corresponding box volumes decrease as a function of the pairlist update frequency, in line with the NVT results.

00e60184-e594-4462-8011-c25cc7328b0c_figure5.gif

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 (purple, red, and orange lines) and the twin-range scheme with constant-force application (TR CFA, light and dark blue lines) and three different update frequencies (np = 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 S7 (Supplementary File 1)

Thermostatting method

The impact of the chosen non-stochastic 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 construction43, 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 or perturb the velocities selectively as done in SD or as suggested in the literature by imposing isokinetic constraints or restraints on each degree of freedom2123.

00e60184-e594-4462-8011-c25cc7328b0c_figure6.gif

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. nm = np = 5). For the single-range (SR) scheme (dashed lines), the mid-range forces were evaluated every time step Δt (i.e. nm = 1), whereas the pairlist was updated every 5th step (i.e. np = 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 for WC, NH and NHchain thermostats. It could be shown that the presence of the artificial density increase at the interfaces is independent of 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. The randomly applied stochastic forces of a SD thermostat seem to disrupt the MTS resonance build-up and thus successfully prevent the occurrence of density artefacts at the interfaces. Based on these findings, the use of TR CFA is not recommended for heterogeneous systems with WC, NH or NHchain thermostats. 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 npt = 40 fs). Thus, the increase in computational cost for SR can be partially compensated by updating the pairlist less frequently (i.e. npt > 10 fs). Of course, there is a limit to the latter, as other errors are introduced by less frequent pairlist updates. By using only a single pairlist, the complexity of the implementation can be reduced, 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.d22294944

Comments on this article Comments (0)

Version 3
VERSION 3 PUBLISHED 05 Nov 2018
Comment
Author details Author details
Competing interests
Grant information
Copyright
Download
 
Export To
metrics
Views Downloads
F1000Research - -
PubMed Central
Data from PMC are received and updated monthly.
- -
Citations
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 3; peer review: 2 approved]. F1000Research 2019, 7:1745 (https://doi.org/10.12688/f1000research.16715.3)
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

Current Reviewer Status: ?
Key to Reviewer Statuses VIEW
ApprovedThe 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 approvedFundamental flaws in the paper seriously undermine the findings and conclusions
Version 3
VERSION 3
PUBLISHED 08 Mar 2019
Revised
Views
13
Cite
Reviewer Report 22 Mar 2019
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
VIEWS 13
The authors have ... Continue reading
CITE
CITE
HOW TO CITE THIS REPORT
Zagrovic B and Braun D. Reviewer Report For: Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations [version 3; peer review: 2 approved]. F1000Research 2019, 7:1745 (https://doi.org/10.5256/f1000research.20224.r45443)
NOTE: it is important to ensure the information in square brackets after the title is included in all citations of this article.
Version 2
VERSION 2
PUBLISHED 18 Jan 2019
Revised
Views
28
Cite
Reviewer Report 29 Jan 2019
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
VIEWS 28
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 ... Continue reading
CITE
CITE
HOW TO CITE THIS REPORT
Zagrovic B and Braun D. Reviewer Report For: Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations [version 3; peer review: 2 approved]. F1000Research 2019, 7:1745 (https://doi.org/10.5256/f1000research.19532.r43226)
NOTE: it is important to ensure the information in square brackets after the title is included in all citations of this article.
  • Author Response 04 Mar 2019
    Sereina Riniker, Laboratory of Physical Chemistry, ETH Zürich, Zurich, 8093, Switzerland
    04 Mar 2019
    Author Response
    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 ... Continue reading
COMMENTS ON THIS REPORT
  • Author Response 04 Mar 2019
    Sereina Riniker, Laboratory of Physical Chemistry, ETH Zürich, Zurich, 8093, Switzerland
    04 Mar 2019
    Author Response
    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 ... Continue reading
Version 1
VERSION 1
PUBLISHED 05 Nov 2018
Views
30
Cite
Reviewer Report 03 Dec 2018
Nico. F. A. van der Vegt, Eduard-Zintl-Institute of Inorganic and Physical Chemistry, TU Darmstadt, Darmstadt, Germany 
Approved
VIEWS 30
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 ... Continue reading
CITE
CITE
HOW TO CITE THIS REPORT
van der Vegt NFA. Reviewer Report For: Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations [version 3; peer review: 2 approved]. F1000Research 2019, 7:1745 (https://doi.org/10.5256/f1000research.18272.r40235)
NOTE: it is important to ensure the information in square brackets after the title is included in all citations of this article.
  • 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
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
Views
43
Cite
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
VIEWS 43
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 ... Continue reading
CITE
CITE
HOW TO CITE THIS REPORT
Zagrovic B and Braun D. Reviewer Report For: Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations [version 3; peer review: 2 approved]. F1000Research 2019, 7:1745 (https://doi.org/10.5256/f1000research.18272.r40231)
NOTE: it is important to ensure the information in square brackets after the title is included in all citations of this article.
  • 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
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

Comments on this article Comments (0)

Version 3
VERSION 3 PUBLISHED 05 Nov 2018
Comment
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
Sign In
If you've forgotten your password, please enter your email address below and we'll send you instructions on how to reset your password.

The email address should be the one you originally registered with F1000.

Email address not valid, please try again

You registered with F1000 via Google, so we cannot reset your password.

To sign in, please click here.

If you still need help with your Google account password, please click here.

You registered with F1000 via Facebook, so we cannot reset your password.

To sign in, please click here.

If you still need help with your Facebook account password, please click here.

Code not correct, please try again
Email us for further assistance.
Server error, please try again.