ACS Publications. Most Trusted. Most Cited. Most Read
Quick View

OR SEARCH CITATIONS

My Activity
Publications
CONTENT TYPES

Figure 1Loading Img
Download Hi-Res ImageDownload to MS-PowerPointCite This:J. Chem. Theory Comput. 2023, 19, 23, 8919-8929
RETURN TO ISSUEPREVBiomolecular SystemsNEXT

Neighbor List Artifacts in Molecular Dynamics Simulations

  • Hyuntae Kim
    Hyuntae Kim
    Department of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany
    International Max Planck Research School on Cellular Biophysics, Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany
    More by Hyuntae Kim
    Orcidhttps://orcid.org/0009-0009-1446-1354
  • , 
  • Balázs Fábián
    Balázs Fábián
    Department of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany
    More by Balázs Fábián
    Orcidhttps://orcid.org/0000-0002-6881-716X
  • , and 
  • Gerhard Hummer*
    Gerhard Hummer
    Department of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany
    Institute of Biophysics, Goethe University Frankfurt, Max-von-Laue-Straße 1, 60438 Frankfurt am Main, Germany
    *Email: [email protected]
    More by Gerhard Hummer
    Orcidhttps://orcid.org/0000-0001-7768-746X
Cite this: J. Chem. Theory Comput. 2023, 19, 23, 8919–8929
Publication Date (Web):November 30, 2023
https://doi.org/10.1021/acs.jctc.3c00777

Copyright © 2023 The Authors. Published by American Chemical Society. This publication is licensed under

CC-BY 4.0.
  • Open Access

Article Views

4879

Altmetric

-

Citations

1
LEARN ABOUT THESE METRICS
PDF (10 MB)
Get e-Alerts
Supporting Info (1)»
SUBJECTS:

Abstract

High Resolution Image
Download MS PowerPoint Slide

Molecular dynamics (MD) simulations are widely used in biophysical research. To aid nonexpert users, most simulation packages provide default values for key input parameters. In MD simulations using the GROMACS package with default parameters, we found large membranes to deform under the action of a semi-isotropically coupled barostat. As the primary cause, we identified overly short outer cutoffs and infrequent neighbor list updates that resulted in missed nonbonded interactions. Small but systematic imbalances in the apparent pressure tensor then induce unphysical asymmetric box deformations that crumple the membrane. We also observed rapid oscillations in averages of the instantaneous pressure tensor components and traced these to the use of a dual pair list with dynamic pruning. We confirmed that similar effects are present in MD simulations of neat water in atomistic and coarse-grained representations. Whereas the slight pressure imbalances likely have minimal impact in most current atomistic MD simulations, we expect their impact to grow in studies of ever-larger systems with coarse-grained representation, in particular, in combination with anisotropic pressure coupling. We present measures to diagnose problems with missed interactions and guidelines for practitioners to avoid them, including estimates for appropriate values for the outer cutoff rl and the number of time steps nstlist between neighbor list updates.

This publication is licensed under

CC-BY 4.0.
  • cc licence
  • by licence

1. Introduction

ARTICLE SECTIONS
Jump To
Molecular dynamics (MD) simulations are a powerful tool to probe molecular processes at a level of detail not currently accessible to experiments. (1) The GROMACS molecular dynamics simulations package (2) is widely used, in particular for applications in biophysics, chemistry, and soft-matter science. It is computationally efficient (2) and easy to use with a wide range of atomistic and coarse-grained force fields quantifying the energetics of molecular interactions. (3,4) Central to its high performance are the nearly linear scaling of the computational cost with system size and its efficient parallelization over multiple computational nodes. (5,6) A key factor for the computational efficiency is the use of neighbor lists containing the pairs of interacting particles. To avoid costly neighbor list updates at every time step, the Verlet scheme includes a buffer of particles between the actual cutoff distance for pair interactions, rc, and an outer cutoff rl > rc. The neighbor list is updated at time intervals chosen so that crossing from distances r > rl to r < rc by ballistic motion is highly unlikely. For the construction of neighbor lists on single-instruction multiple-data (SIMD) hardware architectures, GROMACS implements the MxN algorithm, (7,8) which minimizes internode communication and memory footprint. (2,9) The grouping of particles into spatial clusters by the MxN algorithm enables an efficient evaluation of the real-space pair interactions.
Here, we show that the use of default simulation parameters (10−12) can cause artificial pressure oscillations and broken spatial isotropy. As a consequence, large membrane systems can undergo drastic deformations in the form of unrealistic buckling (Figure 1). We analyze the temporal evolution of lipid bilayers in the NPT ensemble with constant particle number N, pressure P, and temperature T; and of neat solvents in both NPT and NVT ensembles, the latter fixing the volume V instead of the pressure P. As the primary cause, we identify the infrequent construction of the neighbor list as a result of a somewhat too large update interval of nstlist time steps and a somewhat too short outer cutoff distance rl. Consequently, nonbonded interactions are occasionally missed in the force evaluation. The missed interactions cause errors in the elements of the instantaneous pressure tensor. In the NPT ensemble, these errors in the pressure lead to incorrect box rescaling by the barostat, both with the weak-coupling (Berendsen) barostat (13) and the Parrinello–Rahman (PR) barostat. (14) We conclude by providing tools that practitioners can use to detect such problems and guidance for minimizing their impact or avoiding them altogether.

Figure 1

Figure 1. Large Martini POPC bilayer crumples in MD simulations with default simulation parameters, yet stays flat with frequent neighbor list updates. (A) Snapshots of the membrane (phosphate groups in gold) in MD simulation with default parameters. The Verlet-buffer-tolerance, which denotes the maximally allowed energy drift per particle between neighbor list updates due to missed nonbonded interactions, was set to VBT = 0.005 kJ·mol–1·ps–1; the outer cutoff was rl = 1.269 nm; the number of time steps between neighbor list updates was nstlist = 25; and dual pair list was enabled. (B) Snapshots in an MD simulation with neighbor list updates enforced at every time step (nstlist = 1). Snapshots are at time points 50, 150, 250, and 350 ns (left to right). Simulation boxes are indicated as blue lines, and box heights Lz are listed.

High Resolution Image
Download MS PowerPoint Slide
The problems identified here may have afflicted earlier simulation studies. Pointedly, several studies of large membrane systems prevented excessive membrane undulations by restraining the vertical movement of certain lipid head groups with harmonic (15−17) or flat-bottom (18−20) potentials. One can also restrain the box with a weak harmonic potential, for example, by using the plumed software package (21) as a GROMACS plug-in. However, the introduction of such external potentials is unsatisfactory, motivating our efforts to identify and correct the underlying issues.

2. Methods

ARTICLE SECTIONS
Jump To

2.1. Neighbor List and Missed Interactions

In MD simulations, nonelectrostatic nonbonded pair interactions are usually truncated beyond a given distance cutoff rc. (22) Interactions beyond this cutoff are usually estimated analytically (23) assuming a uniform density of particles outside the cutoff sphere, but can be evaluated in Fourier space for power–law potentials in a periodic system using the Ewald method as implemented, e.g., in the particle mesh Ewald (PME) algorithm. (24,25) Without truncation of the real-space interactions, the computational cost of evaluating pairwise forces would scale with the square of the particle number N in the system. With a fixed cutoff rc, the cost scales roughly linearly with N. For the PME algorithm, evaluating the remaining long-range contributions results in N log N scaling.
The neighbor list of a particle contains the indices of its neighboring particles for which the pairwise interactions are explicitly evaluated in real space. In the Verlet scheme, the neighbor list is constructed by searching for neighbors within the cutoff radius rl, with rl ≥ rc. The spherical shell between rl and rc provides a buffer so that neighbor list updates are not required at every time step. Neighbor search requires an evaluation of the pairwise distances, and hence, its computational cost scales at least linearly with the system size. (26) Furthermore, the neighbor search requires internode communications, which can be a major bottleneck for modern hardware architectures. (10,27) If neighbor list updates are performed only every nstlist time steps of length Δt = dt, we expect that some pair interactions are missed because particle pairs move from distances r > rl to r < rc within the time interval nstlist × Δt. For point particles of mass m uniformly distributed in space with number density ρ and moving with velocities following a Maxwell–Boltzmann distribution, we estimate (see Supporting Information (SI) text) the probability that a particular particle misses an interaction as
p missed ≈ 2 π ρ σ 3 ( r c 2 + r l r c + r l 2 ) e − ( r l − r c ) 2 / 2 σ 2 3 ( r l − r c ) 2
(1)
where we assumed that r l − r c ≫ σ = nstlist × Δ t 2 k B T / m , with kB Boltzmann’s constant. For the analyses conducted in this study and in the figures, we instead used eq S8. In the SI Text, we extend the model from point particles to rigid and near-rigid molecules, such as TIP3P water, with an approximate treatment of rigid-body rotations.

2.2. GROMACS Input Parameters

In this study, we critically examine the following GROMACS input parameters: nstlist, nstenergy, nstcalcenergy, nstpcouple, nsttcouple, verlet-buffer-tolerance, and rlist, also denoted as rl. The parameters nstlist, nstenergy, nstcalcenergy, nstpcouple, and nsttcouple denote the number of time steps between neighbor list updates, energy sampling, energy evaluation, barostatting, and thermostatting, respectively. The default values recommended by the developers can be found in the manual: (12) nstlist = 10, nstenergy = 1000 and nstcalcenergy = 100.
The Verlet-buffer-tolerance (VBT) denotes the maximally allowed energy drift per particle between neighbor list updates due to missed nonbonded interactions. Its default value is 0.005 kJ·mol–1·ps–1. (2,12) Changes in VBT result in adjustments of rl and nstlist. In standard GROMACS runs, the values of rl and nstlist are therefore not only system-dependent, but there is also no guarantee that the values are constant throughout a trajectory. In particular, the possible values of nstlist are 20, 25, 40, 50, 60, 80, and 100. The values of the adjusted rl and nstlist can be found in the output log file. To ensure a constant value of nstlist, whether it is user-defined or the default value of 10, VBT must be disabled by setting VBT = −1.
In GROMACS, the maximally allowed energy drift, VBT, determines the frequency of neighbor list updates. By contrast, in LAMMPS (28) the neighbor list is updated when any particle travels more than half the buffer thickness. As the system size increases, the time interval between updates shrinks to the point of forcing an update at every time step. (27)

2.3. MxN Algorithm and Dual Pair List

In SIMD hardware architectures, GROMACS employs the MxN algorithm for a grid-based neighbor search. (7,8) The algorithm clusters a fixed number of particles by gridding the xy plane and binning along the z axis. The clusters with insufficient numbers of particles are filled with dummy particles. The implementation of the algorithm promises a high computational performance. Moreover, the clusters act as another layer of buffer on top of the predefined rl – rc shell, enabling a further increase in nstlist. (8)
The performance can be further improved by implementing a dual pair-list algorithm, (27) using a long outer and a short inner list cutoff. The inner neighbor lists are generated from a pool of particles within the outer list and, hence, updated more frequently. The implementation of the dual pair-list algorithm reduces the overall computational cost of the neighbor search. The update frequencies and the cutoff radii for both the outer and inner list are by default controlled by VBT, and their exact values can be found in the log file. The dual pair-list algorithm can be disabled by setting VBT to −1. When the dual pair-list algorithm is enabled, rl becomes the cutoff radius for the outer neighbor list. For GPUs, dynamic pruning is used to take advantage of their typically large execution width during the neighbor search. (27)

2.4. Membrane Bending Free Energy

The bending energy E associated with elastic deformations of a fluid and incompressible membrane can be estimated as an integral of the squared local mean curvature H over the membrane surface A (29)
E bend = 2 κ ∫ d A H 2
(2)
where κ is the bending rigidity of the membrane. Here, we ignored the contribution of the Gaussian curvature, which is invariant for a given topology. We evaluated the local mean curvature H of the membrane systems using the MemCurv program (30) and then integrated it numerically over the xy plane of the box, thus ignoring curvature corrections to the area element dA.

2.5. Simulation Code

Two versions of GROMACS were examined, namely, 2020.3 and 2023. All numerical analyses were performed using GROMACS 2020.3, the version for which the artifacts were initially observed. However, all of the system types described in this section were also simulated using GROMACS 2023, the latest version available. All the major artifacts caused by the use of inadequate combinations of rl and nstlist were also observed in GROMACS 2023 runs with default parameters. These include the unphysical distortion of the large membrane systems, oscillations in the average instantaneous pressure, and broken spatial isotropy.

2.6. Simulation of Large and Small Martini Membranes

A large membrane system, consisting of 33,282 1-palmitoyl-2-oleoyl-glycero-3-phosphocholine (POPC) lipids and 6,721,594 water particles, was built using the insane.py (31) script and the Martini force field (version 2.2). (4) NaCl salt was added at a concentration of 0.15 M, and 10% of the water particles were replaced with the antifreeze beads (WF). (4) The initial dimensions of the system were 100 nm × 100 nm × 80 nm. The system was equilibrated first in the NVT ensemble for 150 ns and then in the NPT ensemble with semi-isotropic pressure coupling for another 150 ns, using rl = 1.422 nm and nstlist = 20. Production runs of 1 μs length were then performed using the new-rf (11) simulation parameters with rc = 1.1 nm and a 20 fs time step. The system was coupled to a v-rescale thermostat (32) at 310 K and semi-isotropically coupled to a PR barostat with a target pressure of 1 bar (τP = 12 ps). Also, note that nstcalcenergy = 1 was used for all the simulations here unless specified otherwise.
Similarly, a smaller Martini membrane system was built, consisting of 722 POPC lipids and 10,732 water particles. The corresponding initial box dimensions were 15 nm × 15 nm × 10 nm. All preparation procedures and input parameters were identical to those of the large membrane system, and it also underwent a 1 μs long production run.

2.7. Simulation of Water Systems

A system of neat Martini water was prepared for MD simulations in both the NVT and NPT ensembles. An initial volume of 6 nm × 6 nm × 6 nm contained 1530 Martini water particles. Similarly to the Martini membrane systems, the new-rf (11) simulation parameters with a 20 fs time step were used. Also, the ratio between water particles and antifreeze particles WF was set to 9:1. The system was equilibrated for 200 ns in both the NVT and NPT ensembles with rl = 1.422 nm. Detailed analyses were performed on a few μs long production runs. We note that Martini water is, in effect, a Lennard-Jones (LJ) fluid lacking long-range electrostatics.
The distortion of the cubic box containing Martini water was studied in the NPT ensemble to examine the broken spatial isotropy due to the use of an inadequate combination of rl and nstlist. A few varying conditions were examined, where the system was coupled to four different barostat types (semi-isotropically coupled PR, Berendsen and C-rescaling, (33) and anisotropically coupled PR) with consistent target pressures of 1 bar. To detect the effects of possible systematic asymmetries in the calculated pressures as biased box distortions, we used the semi-isotropic and anisotropic pressure coupling schemes, even though these are not normally used to simulate isotropic solvent systems. Before the production runs, the equilibrated systems were isotropically scaled by factors of 0.99, 1.00, and 1.01, respectively. This scaling was intended to mimic possible volume artifacts caused by the use of inadequate combinations of rl and nstlist. To account for possible anisotropy in the initial condition, the equilibrated systems were also rotated about the x, y, and z axes, respectively. This procedure was intended to eliminate any bias caused by the initial configuration of the system. For the same reason, the initial velocities of the particles were randomly generated, according to a Maxwell–Boltzmann distribution, for each of the 500 replicates, resulting in a sample size of 4500 runs for each barostat type. Then, the above procedures were performed using three different combinations: rl = 1.9 nm, nstlist = 1; rl = 1.9 nm, nstlist = 20; and rl = 1.28 nm, nstlist = 25. Hence, 54,000 solvent simulations entered the statistical analysis of possible anisotropy. Finally, hypothesizing the asymmetric implementation of the MxN algorithm to be the cause of the potential anisotropy, we enforced the 1 × 1 atom pair list by recompiling the GROMACS MD engine with -DGMX_SIMD = None. We then performed 4500 replicate simulations of the fully anisotropically coupled system with rl = 1.28 nm, nstlist = 25, and with 1 × 1 atom pair list setup.
Furthermore, the Martini water system was simulated in the NVT ensemble to illustrate that the various observed artifacts (except box rescaling) are not due to the barostat and occur independent of the ensemble type. As for the membrane systems, the respective production runs were 1 μs long. Except for the barostat settings, all input parameters were identical to those used for the NPT Martini water system.
Similarly, a cubic box containing pure TIP3P water (34) was generated via CHARMM-GUI in the NVT ensemble, (35) with dimensions of 5 nm × 5 nm × 5 nm. The system was equilibrated for 15 ns in both the NVT and NPT ensembles using rl = 1.422 nm, while the production runs were 200 ns long with a 2 fs time step. The nonbonded interaction cutoff was rc = 1.2 nm. The temperature was fixed at 310 K with the Nosé–Hoover thermostat. (36,37) For the calculation of the power spectral density of the pressure in TIP3P water, we used the v-rescale thermostat (32) to suppress the oscillatory contributions of the Nosé–Hoover thermostat. Finally, we used the PME algorithm with the default fourierspacing (0.12 nm).

2.8. Power Spectral Analysis

The power spectral density (PSD) of the scalar pressure was calculated using Welch’s method. (38) We used the time series of the pressure as input, which were calculated and saved at every time step (nstenergy = 1). The resulting PSD was plotted as a function of the frequency in units of 1/(nstlist × Δt). The visual inspection focused on peaks in the PSD as a means to identify characteristic time intervals of processes resulting in perturbations of the barostat action.

3. Results

ARTICLE SECTIONS
Jump To

3.1. Unphysical Distortion of the Large Martini Membrane

The temporal evolution of the large Martini membrane system simulated with default parameters and the PR barostat is shown in Figure 1A. Within the 350 ns of MD, the simulation box contracted in the xy membrane plane and expanded in the z direction. The box dimensions changed from 104.0 nm × 104.0 nm × 78.1 nm after equilibration to 97.1 nm × 97.1 nm × 89.7 nm at the end of the production run. This change in box shape left the overall volume approximately constant.
During the simulations, the large membrane buckled to form distinct folds in the x and y directions (Figure 1A). The potential energy and the enthalpy of the system increased by 13,200 and 24,600 kJ·mol–1, respectively, in 350 ns of MD. These increases are substantial also in relative terms, amounting to changes of ≈1.58 and ≈1.85%, respectively. Moreover, the membrane bending energy at the end of the MD simulation was estimated using eq 2 as Ebend ≈ 138 kBT ≈ 357 kJ·mol–1 for a bending rigidity of 25 kBT. (39,40) The observed large increases in the system energy, the enthalpy, and the membrane bending energy strongly indicate that the observed deformation is unphysical.
Pressure imbalances, not the barostats per se, appear to drive the box deformations. In MD simulations with semi-isotropically coupled PR (Figure 1A) and Berendsen barostats (Figure S1), we observed similar box deformations. Having thus ruled out an effect due to a specific barostat, we examined the components of the pressure tensor driving the barostat action. The running averages of the diagonal pressure tensor elements for the default setup (VBT = 0.005 kJ·mol–1·ps–1) are plotted as dashed lines in Figure 2A. Results are shown for the early phase of the simulation when the membrane is still flat (see Figure 1A). We observed that the three diagonal pressure tensor components deviate from the target pressure of 1 bar and each other.

Figure 2

Figure 2. Pressure tensor elements deviate from the target pressure in MD simulations with the default cutoff handling. Running averages of the diagonal elements of the pressure tensor are shown for (A) the large and (B) the small Martini membrane systems in NPT MD simulations and (C) for the Martini water system in an NVT simulation. The left column shows snapshots of the systems. The center and right columns show the running averages evaluated every nstenergy = 1 and 100 steps, respectively. Results obtained with the default simulation parameters for cutoff handling (VBT = 0.005 kJ·mol–1·ps–1) are shown as dashed lines (see the legend for color). The solid lines show results for a larger outer cutoff rl = 1.422 nm with nstlist = 20 fixed and dual pair list disabled. In (A, B), the target pressure of 1 bar in the NPT simulations is indicated by a dashed black line. For the NVT simulation in (C), the dashed black line indicates the consistent average obtained with rl = 1.422 nm and nstlist = 20.

High Resolution Image
Download MS PowerPoint Slide

3.2. Cutoff Handling Is Responsible for Membrane and Box Deformations

Differences in the apparent pressure average, as a function of the frequency of averaging, point to the underlying cause of the unphysical box distortions. The parameter nstenergy is the number of time steps between time points entering the pressure (and energy) averages. It should thus not affect the value of the average. However, for nstenergy = 1, we overestimated the pressure, and for nstenergy = 100, we underestimated it. These differences were significant and reproducible. Moreover, MD simulations of the small membrane system produced similar results (Figure 2B), pointing to the fact that we are dealing with a generic issue. As a possible explanation, we hypothesized that the pressure values calculated between neighbor list updates, which occur at intervals of nstlist = 20 or 25, differ from those right after neighbor list updates, with the former dominating the average for nstenergy = 1 and the latter for nstenergy = 100.
To test this hypothesis, we first performed MD simulations with nstlist = 1 to enforce neighbor list updates after every time step. As shown in Figure 1B, this eliminated the box and membrane deformations and made it possible to simulate a stable large membrane system without additional restraints. As a second test, we performed MD simulations in which we set a large outer cutoff distance of rl = 1.422 nm together with neighbor list updates every nstlist = 20 time steps. As shown by the solid lines in Figure 2A, the pressure components then converged consistently to the target pressure of 1 bar for nstenergy = 1 and 100, respectively. Further support came from MD simulations of the small membrane system, where we also observed convergence to the target pressures (Figure 2B). In addition, MD simulations of a box of Martini water in the NVT ensemble, i.e., without membrane and barostat, showed in essence the same effect of apparent deviations between the pressure averages with default cutoff settings and nstenergy = 1 and 100, and consistent averages with rl = 1.422 nm and nstlist = 20 (Figure 2C). Infrequent neighbor list updates thus emerged as a likely cause of pressure imbalances and associated box deformations.

3.3. Instantaneous Pressure Oscillates Between Neighbor List Updates

As a further test of the hypothesis that unresolved cutoff violations are at the heart of the observed problems, we examined the pressure as a function of the time between neighbor list updates (Figure 3). For the large and small membrane system in NPT ensembles and for the Martini water system in an NVT ensemble (top to bottom), we saved the time series of the pressure tensor components and averaged them over blocks of length nstlist starting immediately after a neighbor list update (time step 0). Results are shown for the default cutoff setting with nstlist = 25 (three left columns) and for the setting with rl = 1.422 nm and nstlist = 20 (right column). We averaged the in-plane pressure components P∥ = (Pxx + Pyy)/2 and show the normal pressure as P⊥ = Pzz.

Figure 3

Figure 3. Average pressure tensor components deviate from the target pressure between neighbor list updates. Averages were performed over blocks of nstlist time steps starting immediately after a neighbor list update (time step 0). Results are shown for the large (top) and small (center) membrane systems in NPT simulations and for the NVT Martini water system (bottom). Results include the default cutoff setting with nstlist = 25 (three left columns) and the setting with rl = 1.422 nm and nstlist = 20 (right column). The default rl values for the NPT large and small membranes and the NVT solvent were 1.269, 1.267, and 1.28 nm, respectively. For the runs in column 3, the dual pair list was disabled. We averaged the lateral pressure components P∥ = (Pxx + Pyy)/2 (orange curves) and showed the normal pressure as P⊥ = Pzz (blue curves). Dashed horizontal lines of matching colors denote the corresponding overall averages. Black horizontal dotted lines represent the reference pressure values. Cyan circles (left column) indicate deviations of the averaged pressures from the target. Green circles (column 3) indicate deviations just before the neighbor list update.

High Resolution Image
Download MS PowerPoint Slide
Consistent with our expectations, we found that with default cutoff settings, the average pressures immediately after a neighbor list update (i.e., at time step zero in Figure 3) are systematically lower than the pressures calculated between neighbor list updates (time steps > 0). However, we were initially puzzled by the observation that even at time step zero, the average pressure deviated from the target pressure (even though this is consistent with the findings in Figure 2 for nstenergy = 100). As a possible explanation, we considered that neighbor list updates came after nstlist = 25 time steps, yet thermostat and barostat actions after nsttcouple = nstpcouple = 20 time steps, and thus asynchronously (circles in the left column of Figure 3). Indeed, by setting nsttcouple = nstpcouple = nstlist = 25 time steps, the target and average pressures at time step 0 are consistent.
Counter to our expectations, we found the instantaneous pressure values to rise rapidly with time and to exhibit distinct oscillations. For missed interactions because of inadequate rl, we had expected a delayed and monotonic rise with time. As a source of this unexpected behavior, we identified the use of a dual pair list. When we disabled the dual pair-list evaluation by setting VBT = −1, both the rapid initial rise and the oscillations in the average pressures disappeared (column 3 in Figure 3). With the dual pair list enabled, the outer and the inner neighbor lists are maintained with two distinct intervals of 25 and 4 integration time steps, respectively. The zigzag oscillations in Figure 3 appear to be a superposition of two curves with these two periods. This argument is supported by simulations of Martini water in the NVT ensemble using a different hardware architecture, where the outer and inner lists were updated every 25 and 5 time steps, respectively (Figure S2). With the update intervals of the outer and the inner lists being multiples of five, the dominant oscillation period is also five time steps.

3.4. Pressure Deviations Correlate with Missed Particle Interactions

In Figure 3, the average pressure values right after neighbor list updates are close to the target values. However, at the time step just before the neighbor list update (at time step 24 in Figure 3, column 3), we noticed significant positive deviations of ΔP from the pressure at time step 0, as indicated by green circles. In Figure 4, we plot these deviations for lateral (ΔP∥) and normal pressures (ΔP⊥) as functions of the outer cutoff rl with the dual pair list disabled. Results are shown for Martini and TIP3P water. For reference, we also show nmissed for Martini water and nH–H, nO–H and nO–O for TIP3P water: nmissed is the expected mean number of unique cutoff violations of the water particles (eq S8); nH–H, nO–H and nO–O denote the expected number of unique cutoff violations of the hydrogen–hydrogen, oxygen–hydrogen, and oxygen–oxygen atom pairs, respectively, using eqs S8 and S11. Except for the shortest rl ≈ rc, we find that the pressure errors just before neighbor list updates follow the trend of missed interactions. We note further that the ΔP errors are positive for Martini water, consistent with missed attractive Lennard-Jones interactions, as these lead to an underestimation of the cohesiveness. We thus conclude that missed interactions due to overly short outer cutoffs are the primary contributor to deviations ΔP in the pressure.

Figure 4

Figure 4. Differences ΔP in the pressures calculated just before and right after neighbor list updates follow the trend of missed interactions. ΔP (left scale) for lateral (orange) and normal pressures (blue) are shown as functions of rl (A) for Martini water and (B) for TIP3P water, both simulated in NVT conditions with nstlist = 20 and dual pair list disabled. Standard deviations are shown as error bars. Also shown (right scale) is the expected number of missed interactions (black dashed lines). Considering the Martini water particles as point particles, the number of missed interactions (nmissed) is evaluated using eq S8. For TIP3P water, the expected number of the missed hydrogen–hydrogen (nH–H), oxygen–hydrogen (nO–H), and oxygen–oxygen (nO–O) electrostatic interactions are evaluated using eqs S8 and S11.

High Resolution Image
Download MS PowerPoint Slide
For TIP3P water, we found that the difference in the apparent pressure at time points just before and right after neighbor list updates tends to be negative, ΔP < 0 (Figures 4B and S3). The negative sign indicates that for TIP3P missed repulsive interactions dominate. From the rigid-rotor model described in the SI text, we indeed expect that at short times, the fast-moving hydrogen atoms with their low mass will result in missed repulsive hydrogen–hydrogen real-space electrostatic interactions. This is in contrast to Martini water, where ΔP > 0 results from missed attractive LJ interactions (Figure 4A). The nonmonotonic dependence of ΔP on rl for TIP3P water (rl ≲ 1.23 nm in Figure S3A and rl ≲ 1.28 nm in Figure S3B) may be caused by a partial cancellation of contributions from missed attractive (opposite-charge and LJ) and repulsive (like-charge) interactions. Moreover, for small buffers rl – rc and long time intervals between the neighbor list updates nstlist × Δt, the assumptions leading to eq S8 may no longer hold. A refined model could, for instance, use the estimated distribution of distances for missed pair interactions to estimate the impact on the virial.

3.5. Anisotropic Errors in Pressure Tensor Deform Box Shape

The errors ΔP in the pressure shown in Figure 4 tend to be somewhat anisotropic, even though the boxes had fixed cubic shape and volume. The lateral errors tend to be somewhat smaller than the normal errors, ΔP∥ < ΔP⊥. Therefore, we hypothesized that in MD simulations with both semi-isotropic and anisotropic barostats, we should see a tendency for the boxes to grow in the z direction. We tested this hypothesis by running repeated MD simulations of Martini water systems with semi-isotropic and fully anisotropic barostats of PR, Berendsen, and C-rescale type.
The results of these simulations confirm a tendency for the simulation box to expand preferentially along z by the action of the barostat (Table 1). This tendency is significant (p-value < 0.01 calculated from Pearson’s chi-squared test) for fully anisotropic pressure coupling with the PR barostat and two relatively poorer cutoff settings (rl = 1.9 nm, nstlist = 20 and rl = 1.28 nm, nstlist = 25). The observed tendency of the anisotropically coupled box to expand along z (Table 1) is consistent with the consistently larger error in the pressure along z seen in Figure 4. We note, however, that the direction of the box expansion is biased toward z but not deterministic. Spatial isotropy was restored when rl = 1.9 nm and nstlist = 1 were used (p-value ≈ 0.89). For this setting, all interactions should be counted at all time steps, and the dual pair list is turned off.
Table 1. Spatial Isotropy Can be Compromised by the Cutoff Treatment in NPT MD Simulations of Martini Watera
  semi-isotropic anisotropic
  Parrinello-Rahman Berendsen C-rescale Parrinello-Rahman
  contraction: elongation X: Y: Z
rl (nstlist) p-value p-value
  1992:2508 2274:2226 2232:2268 1486:1489:1525
1.9 nm (1) ≪0.001 0.484 0.709 0.890
  1310:3190 1997:2503 2239:2261 1387:1536:1577
1.9 nm (20) ≪0.001 ≪0.001 0.634 0.004
  1650:2850 1977:2523 2154:2346 1540:1322:1638
1.28 nm (25) ≪0.001 ≪0.001 ≪0.001 ≪0.001
1.28 nm (25)   1496:1531:1473
1 × 1 pair list   0.768
a

Four pressure coupling schemes were examined: semi-isotropically coupled PR, Berendsen, and C-rescale barostats, and an anisotropically coupled PR barostat. Four different combinations of rl and nstlist were tested, including a 1 × 1 atom pair list (column 1). Columns 2–5 list the number of times the system elongated along a specific principal axis. P-values were calculated under the null hypothesis that the probabilities are equal to 1/2 for contraction and elongation along z in the semi-isotropic case, and equal to 1/3 for expansions along x, y, and z in the fully anisotropic case

Interestingly, we observed a statistically significant asymmetry in the box shape changes even for rl = 1.9 nm and nstlist = 1 for MD simulations with a semi-isotropically coupled PR barostat (Table 1). Therefore, we checked alternative barostats with the same setting. For simulations with semi-isotropically coupled Berendsen and C-rescale barostats, the setting rl = 1.9 nm and nstlist = 1 did not result in any significant bias in our tests (Table 1). Therefore, we suspect that the asymmetry in box shape changes with semi-isotropically coupled PR barostat may be a result of the specifics of the barostat implementation. For less stringent settings, rl = 1.9 nm and nstlist = 20, we found that also the Berendsen barostat induced significant shape asymmetries but not the C-rescale barostat (Table 1). This further hints at the possibility that subtle differences in the barostat algorithm or implementation make the MD simulations susceptible to anisotropies in the pressure tensor. In practice, semi-isotropic coupling is typically used for systems with mechanical resistance, such as a lipid bilayer spanning the xy plane. For such systems, the consequences of the asymmetries detected here should be negligible.
A likely source of anisotropies is the asymmetric implementation of the MxN algorithm. (2,8) The algorithm constructs the neighbor lists by gridding the xy plane and binning the particles along the z axis. Clusters with insufficient numbers of particles are filled with dummy particles, which have zero contribution on the cluster volume. (8) Hence, the average cluster dimensions along the z axis would be smaller than those along the lateral axes. This effectively results in an anisotropic buffer thickness. Consequently, we expect that a larger fraction of cohesive interactions is missed along the z axis, which would explain the observed spatial anisotropies. To test this hypothesis, we enforced a 1 × 1 atom pair list (one particle per cluster) by recompiling the GROMACS MD engine with -DGMX_SIMD = None. We then performed 4500 replicate MD simulations of the fully anisotropically coupled system containing pure Martini water with rl = 1.28 nm and nstlist = 25, as before. Consistent with our hypothesis, we found that by enforcing a 1 × 1 atom pair list, the asymmetry in the box shape changes disappeared (last row in Table 1).

4. Discussion

ARTICLE SECTIONS
Jump To
Incomplete neighbor lists resulting in missed pair interactions emerge as the primary cause of the various artifacts observed. Fast-moving particles occasionally cross from outside the outer cutoff, r > rl, to distances within the cutoff sphere, r < rc, of interaction partners between updates of the neighbor list. As a result, some nonbonded interactions are missed in the force evaluations at proceeding time steps up to the next neighbor list update. Affected are all interactions that rely on the neighbor list. This includes the real-space part of Coulomb interactions evaluated with the PME algorithm, (24) as shown for TIP3P water. It should also affect the real-space part of other power–law interactions evaluated with the PME algorithm. (24,25)
As a partial fix to stabilize the simulations in cases where the outer cutoff rl is too short, one can apply barostats (and thermostats) immediately after neighbor list updates and thus with all interactions accounted for. To leave room for efficiency optimizations for specific MD runs, the actual combination of rl and nstlist in GROMACS is determined by the Verlet-buffer-tolerance VBT. In addition, GROMACS does not ensure that nsttcouple and nstpcouple are identical to nstlist by default. All barostat types in GROMACS act with a fixed frequency (nstpcouple) according to the instantaneous pressure tensor. If nstpcouple is not an integer multiple of nstlist and if the outer cutoff rl is too short, then the barostat will at regular intervals operate with instantaneous pressures calculated with an incomplete list of pair interactions. In the case of Martini water, the missed interactions are attractive. Consequently, the barostat then acts to reduce the artificially high apparent pressure by increasing the system volume. However, such an increase in volume is artificial, the effect of which is only partially restored at those time points where barostat action immediately follows a neighbor list update. It is thus advantageous to set nsttcouple = nstpcouple = n × nstlist with integer n ≥ 1. Indeed, the unphysical distortion of the large Martini membrane system is greatly suppressed if nstpcouple = nstlist even with a relatively short outer cutoff of rl = 1.269 nm (Figure S4). However, the undulation of the membrane shown in Figure S4 is still noticeably more pronounced than that simulated with nstlist = 1, resulting in a slightly larger Lz = 78.2 nm, compared to Lz = 77.9 nm in Figure 1B. Hence, setting nsttcouple = nstpcouple = nstlist in combination with default rl alone may be insufficient.
Setting the outer cutoff rl to values somewhat larger than the default further stabilizes the simulations. In numerical tests, we found it sufficient to obtain the default values of rl and nstlist for the same system but with double the integration time step, and then to enforce these values manually together with nstpcouple = nstlist = nsttcouple. With VBT set to −1, this procedure also automatically disables dual pair list evaluation.

5. Recommendations for Practitioners

ARTICLE SECTIONS
Jump To

5.1. Diagnosis

Deviations between target and calculated pressures in NPT simulations serve as the primary indicator of possible issues with neighbor list construction and cutoff treatment. As shown in Figure 2, missed interactions as a result of an inadequate cutoff treatment tend to result in noticeable deviations between target and actual pressures and in small but again noticeable anisotropy as manifested by differences among the diagonal elements of the pressure tensor.
As a complication, the pressure and energy are calculated only at every nstcalcenergy time step, whereas neighbor lists are updated at every nstlist time step. If nstcalcenergy is an integer multiple of nstlist, the pressure tensor is always evaluated immediately after a neighbor list update. This synchrony can mask deviations (Figure 2B right). More frequent or asynchronous pressure calculation (e.g., by setting nstenergy = 1) can then reveal cutoff issues, as seen by comparing Figure 2B center and right. However, variations of nstlist along atomistic MD trajectories can further complicate the analysis.
Trial trajectories with pressures calculated every time step (nstenergy = 1) provide the basis for a more detailed analysis. Running averages (Figure 2) and averages over blocks of nstlist time steps (Figure 3) help to pinpoint problems with missed interactions and the resulting pressure imbalances.
Oscillations in the scalar pressure as a result of neighbor list updates and barostat action can be revealed by a power spectral analysis. For an inadequate outer cutoff, the PSD of the scalar pressure shows distinct spikes at frequencies that are integer multiples of 1/(nstlist × Δt), as shown in Figure 5. Their amplitude is modulated by the update frequency of the inner neighbor list. By enforcing a larger outer cutoff rl with VBT set to −1 and without dual cutoff, these oscillations and the resulting features in the power spectral density disappear (lower curve in Figure 5). For the pressure in MD simulations of TIP3P water, we similarly observed distinct spikes in the PSD at integer multiples of 1/(nstlist × Δt), which disappeared for large values of rl and with the MxN algorithm disabled (Figure S5). Note that we used the v-rescale thermostat to avoid the oscillatory contributions to the pressure PSD in simulations with the Nosé–Hoover thermostat. (36,37)

Figure 5

Figure 5. Power spectral density (PSD) of the scalar pressure calculated for Martini water in the NVT ensemble. Results are shown for rl = 1.27 nm and nstlist = 25 (top curves; scaled by factor 100) with dual pair list and inner neighbor list updates every four (orange) and five (blue line) integration steps, respectively. Except for the oscillations at frequency multiples of 1/(nstlist × Δt), the two curves nearly superimpose. Also shown (fluorescent green, bottom curve) is the PSD for a simulation with rl = 1.422 nm and nstlist = 25 without dual cutoff, which does not show oscillations.

High Resolution Image
Download MS PowerPoint Slide

5.2. Recommendations

For the existing software, pending updates, we recommend to manually define the outer cutoff rl and the neighbor list update frequency nstlist as the default values obtained for the identical system with twice the time step. This requires setting verlet-buffer-tolerance = −1. For Martini systems using the new-rf parameters with nstlist = 20, rl should be at least 1.35 nm. This can be achieved by setting VBT = ∼0.0002 kJ·mol–1·ps–1 (Figure S6). For atomistic systems with nstlist = 20, rl should be at least 1.3 nm. Finally, nsttcouple and nstpcouple should be exact multiples of nstlist. However, these measures may significantly reduce the computational performance, especially for large systems. Depending on the hardware specifications, we observed a decrease in performance of up to ∼30%.
In the future, to minimize the impact on computational cost, a spatially isotropic neighbor search algorithm is desired also for the SIMD architecture. Spatial isotropy can for instance be restored to a significant degree with minimum computational overhead if the axis of the one-dimensional search is cycled through x, y, and z instead of keeping it fixed at z. More conservative default choices of rl and nstlist will also help to ensure stable MD simulations. The choice of rl and nstlist can be based on simple expressions for the probability of missed interactions such as eqs S8 and S11, including for atomistic MD simulations. Finally, all nstX variables (where X is energy, tcouple, pcouple, etc.) should be exact multiples of nstlist. In this way, energy output as well as thermostat and barostat actions benefit from freshly updated neighbor lists. If the dual pair list is enabled, the interval between updates of inner and outer neighbor lists should be set similarly to ensure that barostats, in particular, act with all interactions considered.

6. Concluding Remarks

ARTICLE SECTIONS
Jump To
The efficient calculation of pair interactions is at the heart of modern MD simulation codes. The construction of neighbor lists is a critical factor to reduce computational cost and to achieve near-linear scaling with system size, as well as efficient parallelization. We found that even a small number of missed nonbonded interactions can impact the accuracy and qualitative behavior of MD simulations. For large membrane systems, we observed that a slight but systematic error in the pressure with default settings for cutoff distance and neighbor list updates resulted in artificial membrane buckling caused by the counteraction of the barostat. We suspect that similar behavior motivated the imposition of restraints on lipid motion normal to the membrane in earlier studies. (15−20) Slight but systematic anisotropies in the errors of the pressure result in anisotropic deformations of the box shape. We also observed distinct beating effects in the pressure time series. As underlying causes of these different but related problems, we identified (i) missed pair interactions as a result of too infrequent neighbor list updates, (ii) slight anisotropies in the neighbor list construction, and (iii) neighbor list and barostat updates at incommensurable time intervals. In most current MD simulations, in particular, at all-atom resolution, we expect the slight pressure imbalances to have minimal impact. However, as MD simulations are used to study ever-larger systems, such as the coarse-grained membrane systems in Figure 1, the issues will have to be addressed. Of particular concern are large systems with inherent anisotropies such as membrane systems or systems containing quasi-infinite molecules and molecular assemblies, including DNA and protein fibers spanning the periodic box, as these could become deformed by the action of the barostat in response to small but systematic errors in the pressure tensor components. Particular care should be taken in calculations of quantities related to stress or pressure differences such as interfacial tension and line tension. Whereas the immediate fixes of longer outer cutoff lengths rl and disabled dual pair lists (VBT = −1) are associated with increased computational cost, we are confident that with the root causes identified, adjustments in algorithms and code can be made that will resolve the issues without major computational overhead.

Supporting Information

ARTICLE SECTIONS
Jump To

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.3c00777.

  • Detailed derivation of number of missed interactions; the rigid-rotor model for missed interactions in systems with rigid or near-rigid molecules; power spectral density of pressure fluctuations in TIP3P water; the shape of the average pressure tensor between neighbor list updates depends on the update frequencies of the inner and outer list; difference ΔP in the pressure just before and right after a neighbor list update in the NVT TIP3P solvent system; in a simulation of the large Martini membrane system with nstpcouple = nstlist = 25, the unphysical box distortion is greatly suppressed; power spectral density (PSD) of the pressure fluctuations in MD simulations of TIP3P water in the NVT ensemble with nstlist = 20; probability pmissed of missed interactions (gray line) as function of the VBT parameter in GROMACS for the NVT Martini water system (PDF)

Terms & Conditions

Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

Author Information

ARTICLE SECTIONS
Jump To
  • Corresponding Author
    • Gerhard Hummer - Department of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany;  Institute of Biophysics, Goethe University Frankfurt, Max-von-Laue-Straße 1, 60438 Frankfurt am Main, Germany;  Orcidhttps://orcid.org/0000-0001-7768-746X;  Email: [email protected]
  • Authors
    • Hyuntae Kim - Department of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany;  International Max Planck Research School on Cellular Biophysics, Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany;  Orcidhttps://orcid.org/0009-0009-1446-1354
    • Balázs Fábián - Department of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany;  Orcidhttps://orcid.org/0000-0002-6881-716X
  • Funding

    Open access funded by Max Planck Society.

  • Notes
    The authors declare no competing financial interest.

Acknowledgments

ARTICLE SECTIONS
Jump To

The authors thank the International Max Planck Research School (IMPRS) on Cellular Biophysics (H.K.), the Max Planck Society (H.K., B.F., G.H.), and the Alexander von Humboldt-Foundation (B.F.) for their support. The authors thank Sebastian Kehl (MPDCF) for discussions about the gmx source code. The authors thank the GROMACS developers for extensive and insightful discussions about the potential updates in upcoming GROMACS versions and are glad to acknowledge that the issues listed in this manuscript will be reflected upon in GROMACS 2024.

References

ARTICLE SECTIONS
Jump To

This article references 40 other publications.

  1. 1
    Gelpi, J.; Hospital, A.; Goñi, R.; Orozco, M. Molecular dynamics simulations: Advances and applications. Adv. Appl. Bioinf. Chem. 2015, 8, 37– 47,  DOI: 10.2147/AABC.S70333
    Google Scholar
    1
    Molecular dynamics simulations: advances and applications
    Hospital Adam; Goni Josep Ramon; Orozco Modesto; Gelpi Josep L
    Advances and applications in bioinformatics and chemistry : AABC (2015), 8 (), 37-47 ISSN:1178-6949.
    Molecular dynamics simulations have evolved into a mature technique that can be used effectively to understand macromolecular structure-to-function relationships. Present simulation times are close to biologically relevant ones. Information gathered about the dynamic properties of macromolecules is rich enough to shift the usual paradigm of structural bioinformatics from studying single structures to analyze conformational ensembles. Here, we describe the foundations of molecular dynamics and the improvements made in the direction of getting such ensemble. Specific application of the technique to three main issues (allosteric regulation, docking, and structure refinement) is discussed.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BC28rgs1ShtQ%253D%253D&md5=1b8c8639bd3356dfb4d700a3ec325729
  2. 2
    Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1–2, 19– 25,  DOI: 10.1016/j.softx.2015.06.001
    Google Scholar
    There is no corresponding record for this reference.
  3. 3
    Lundborg, M.; Lindahl, E. Automatic GROMACS Topology Generation and Comparisons of Force Fields for Solvation Free Energy Calculations. J. Phys. Chem. B 2015, 119, 810– 823,  DOI: 10.1021/jp505332p
    Google Scholar
    3
    Automatic GROMACS Topology Generation and Comparisons of Force Fields for Solvation Free Energy Calculations
    Lundborg, Magnus; Lindahl, Erik
    Journal of Physical Chemistry B (2015), 119 (3), 810-823CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
    Free energy calcn. has long been an important goal for mol. dynamics simulation and force field development, but historically it has been challenged by limited performance, accuracy, and creation of topologies for arbitrary small mols. This has made it difficult to systematically compare different sets of parameters to improve existing force fields, but in the past few years several authors have developed increasingly automated procedures to generate parameters for force fields such as Amber, CHARMM, and OPLS. The authors present a new framework that enables fully automated generation of GROMACS topologies for any of these force fields and an automated setup for parallel adaptive optimization of high-throughput free energy calcn. by adjusting lambda point placement on the fly. As a small example of this automated pipeline, the authors have calcd. solvation free energies of 50 different small mols. using the GAFF, OPLS-AA, and CGenFF force fields and four different water models, and by including the often neglected polarization costs, the authors show that the common charge models are somewhat under-polarized.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2cXhvVWns7zF&md5=3be571dedfcae41667937c7108217be4
  4. 4
    Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812– 7824,  DOI: 10.1021/jp071097f
    Google Scholar
    4
    The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations
    Marrink, Siewert J.; Risselada, H. Jelger; Yefimov, Serge; Tieleman, D. Peter; De Vries, Alex H.
    Journal of Physical Chemistry B (2007), 111 (27), 7812-7824CODEN: JPCBFK; ISSN:1520-6106. (American Chemical Society)
    We present an improved and extended version of our coarse grained lipid model. The new version, coined the MARTINI force field, is parametrized in a systematic way, based on the reprodn. of partitioning free energies between polar and apolar phases of a large no. of chem. compds. To reproduce the free energies of these chem. building blocks, the no. of possible interaction levels of the coarse-grained sites has increased compared to those of the previous model. Application of the new model to lipid bilayers shows an improved behavior in terms of the stress profile across the bilayer and the tendency to form pores. An extension of the force field now also allows the simulation of planar (ring) compds., including sterols. Application to a bilayer/cholesterol system at various concns. shows the typical cholesterol condensation effect similar to that obsd. in all atom representations.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD2sXmsVKmsLc%253D&md5=428d5750b94652e4917d905a30658235
  5. 5
    Gecht, M.; Siggel, M.; Linke, M.; Hummer, G.; Köfinger, J. MDBenchmark: A toolkit to optimize the performance of molecular dynamics simulations. J. Chem. Phys. 2020, 153, 144105  DOI: 10.1063/5.0019045
    Google Scholar
    5
    MDBenchmark: A toolkit to optimize the performance of molecular dynamics simulations
    Gecht, Michael; Siggel, Marc; Linke, Max; Hummer, Gerhard; Koefinger, Juergen
    Journal of Chemical Physics (2020), 153 (14), 144105CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
    Despite the impending flattening of Moores law, the system size, complexity, and length of mol. dynamics (MD) simulations keep on increasing, thanks to effective code parallelization and optimization combined with algorithmic developments. Going forward, exascale computing poses new challenges to the efficient execution and management of MD simulations. The diversity and rapid developments of hardware architectures, software environments, and MD engines make it necessary that users can easily run benchmarks to optimally set up simulations, both with respect to time-to-soln. and overall efficiency. To this end, we have developed the software MDBenchmark to streamline the setup, submission, and anal. of simulation benchmarks and scaling studies. The software design is open and as such not restricted to any specific MD engine or job queuing system. To illustrate the necessity and benefits of running benchmarks and the capabilities of MDBenchmark, we measure the performance of a diverse set of 23 MD simulation systems using GROMACS 2018. We compare the scaling of simulations with the no. of nodes for central processing unit (CPU)-only and mixed CPU-graphics processing unit (GPU) nodes and study the performance that can be achieved when running multiple simulations on a single node. In all these cases, we optimize the nos. of message passing interface (MPI) ranks and open multi-processing (OpenMP) threads, which is crucial to maximizing performance. Our results demonstrate the importance of benchmarking for finding the optimal system and hardware specific simulation parameters. Running MD simulations with optimized settings leads to a significant performance increase that reduces the monetary, energetic, and environmental costs of MD simulations. (c) 2020 American Institute of Physics.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3cXitVCmurvF&md5=d8cd0cbf9ed72f74d83f9170dacd1a11
  6. 6
    Kutzner, C.; Páll, S.; Fechner, M.; Esztermann, A.; de Groot, B. L.; Grubmüller, H. Best bang for your buck: GPU nodes for GROMACS biomolecular simulations. J. Comput. Chem. 2015, 36, 1990– 2008,  DOI: 10.1002/jcc.24030
    Google Scholar
    6
    Best bang for your buck: GPU nodes for GROMACS biomolecular simulations
    Kutzner, Carsten; Pall, Szilard; Fechner, Martin; Esztermann, Ansgar; de Groot, Bert L.; Grubmueller, Helmut
    Journal of Computational Chemistry (2015), 36 (26), 1990-2008CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
    The mol. dynamics simulation package GROMACS runs efficiently on a wide variety of hardware from commodity workstations to high performance computing clusters. Hardware features are well-exploited with a combination of single instruction multiple data, multithreading, and message passing interface (MPI)-based single program multiple data/multiple program multiple data parallelism while graphics processing units (GPUs) can be used as accelerators to compute interactions off-loaded from the CPU. Here, we evaluate which hardware produces trajectories with GROMACS 4.6 or 5.0 in the most economical way. We have assembled and benchmarked compute nodes with various CPU/GPU combinations to identify optimal compns. in terms of raw trajectory prodn. rate, performance-to-price ratio, energy efficiency, and several other criteria. Although hardware prices are naturally subject to trends and fluctuations, general tendencies are clearly visible. Adding any type of GPU significantly boosts a node's simulation performance. For inexpensive consumer-class GPUs this improvement equally reflects in the performance-to-price ratio. Although memory issues in consumer-class GPUs could pass unnoticed as these cards do not support error checking and correction memory, unreliable GPUs can be sorted out with memory checking tools. Apart from the obvious determinants for cost-efficiency like hardware expenses and raw performance, the energy consumption of a node is a major cost factor. Over the typical hardware lifetime until replacement of a few years, the costs for elec. power and cooling can become larger than the costs of the hardware itself. Taking that into account, nodes with a well-balanced ratio of CPU and consumer-class GPU resources produce the max. amt. of GROMACS trajectory over their lifetime. © 2015 The Authors. Journal of Computational Chem. Published by Wiley Periodicals, Inc.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2MXht1Krur%252FM&md5=2effcd08bab28e1c572dc1dac25e9c01
  7. 7
    Páll, S.; Abraham, M. J.; Kutzner, C.; Hess, B.; Lindahl, E. Tackling Exascale Software Challenges in Molecular Dynamics Simulations with GROMACS. In Lecture Notes in Computer Science; Springer, 2015; pp 3– 27.
    Google Scholar
    There is no corresponding record for this reference.
  8. 8
    Páll, S.; Hess, B. A flexible algorithm for calculating pair interactions on SIMD architectures. Comput. Phys. Commun. 2013, 184, 2641– 2650,  DOI: 10.1016/j.cpc.2013.06.003
    Google Scholar
    8
    A flexible algorithm for calculating pair interactions on SIMD architectures
    Pall, Szilard; Hess, Berk
    Computer Physics Communications (2013), 184 (12), 2641-2650CODEN: CPHCBZ; ISSN:0010-4655. (Elsevier B.V.)
    Calcg. interactions or correlations between pairs of particles is typically the most time-consuming task in particle simulation or correlation anal. Straightforward implementations using a double loop over particle pairs have traditionally worked well, esp. since compilers usually do a good job of unrolling the inner loop. In order to reach high performance on modern CPU and accelerator architectures, single-instruction multiple-data (SIMD) parallelization has become essential. Avoiding memory bottlenecks is also increasingly important and requires reducing the ratio of memory to arithmetic operations. Moreover, when pairs only interact within a certain cut-off distance, good SIMD utilization can only be achieved by reordering input and output data, which quickly becomes a limiting factor. Here we present an algorithm for SIMD parallelization based on grouping a fixed no. of particles, e.g. 2, 4, or 8, into spatial clusters. Calcg. all interactions between particles in a pair of such clusters improves data reuse compared to the traditional scheme and results in a more efficient SIMD parallelization. Adjusting the cluster size allows the algorithm to map to SIMD units of various widths. This flexibility not only enables fast and efficient implementation on current CPUs and accelerator architectures like GPUs or Intel MIC, but it also makes the algorithm future-proof. We present the algorithm with an application to mol. dynamics simulations, where we can also make use of the effective buffering the method introduces.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC3sXhtFCht7zM&md5=6ec6327a4dc3fbc551350132aef15485
  9. 9
    Buyya, R.; Vecchiola, C.; Selvi, S. T. Principles of Parallel and Distributed Computing. In Mastering Cloud Computing; Morgan Kaufmann: Boston, 2013; Chapter 2, pp 29– 70.
    Google Scholar
    There is no corresponding record for this reference.
  10. 10
    Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435– 447,  DOI: 10.1021/ct700301q
    Google Scholar
    10
    GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation
    Hess, Berk; Kutzner, Carsten; van der Spoel, David; Lindahl, Erik
    Journal of Chemical Theory and Computation (2008), 4 (3), 435-447CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
    Mol. simulation is an extremely useful, but computationally very expensive tool for studies of chem. and biomol. systems. Here, we present a new implementation of our mol. simulation toolkit GROMACS which now both achieves extremely high performance on single processors from algorithmic optimizations and hand-coded routines and simultaneously scales very well on parallel machines. The code encompasses a minimal-communication domain decompn. algorithm, full dynamic load balancing, a state-of-the-art parallel constraint solver, and efficient virtual site algorithms that allow removal of hydrogen atom degrees of freedom to enable integration time steps up to 5 fs for atomistic simulations also in parallel. To improve the scaling properties of the common particle mesh Ewald electrostatics algorithms, we have in addn. used a Multiple-Program, Multiple-Data approach, with sep. node domains responsible for direct and reciprocal space interactions. Not only does this combination of algorithms enable extremely long simulations of large systems but also it provides that simulation performance on quite modest nos. of std. cluster nodes.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD1cXhsVSqurc%253D&md5=d53c94901386260221792ea30f151c5f
  11. 11
    de Jong, D. H.; Baoukina, S.; Ingólfsson, H. I.; Marrink, S. J. Martini straight: Boosting performance using a shorter cutoff and GPUs. Comput. Phys. Commun. 2016, 199, 1– 7,  DOI: 10.1016/j.cpc.2015.09.014
    Google Scholar
    11
    Martini straight: Boosting performance using a shorter cutoff and GPUs
    de Jong, Djurre H.; Baoukina, Svetlana; Ingolfsson, Helgi I.; Marrink, Siewert J.
    Computer Physics Communications (2016), 199 (), 1-7CODEN: CPHCBZ; ISSN:0010-4655. (Elsevier B.V.)
    In mol. dynamics simulations, sufficient sampling is of key importance and a continuous challenge in the field. The coarse grain Martini force field has been widely used to enhance sampling. In its original implementation, this force field applied a shifted Lennard-Jones potential for the non-bonded van der Waals interactions, to avoid problems related to a relatively short cutoff. Here we investigate the use of a straight cutoff Lennard-Jones potential with potential modifiers. Together with a Verlet neighbor search algorithm, the modified potential allows the use of GPUs to accelerate the computations in Gromacs. We find that this alternative potential has little influence on most of the properties studied, including partitioning free energies, bulk liq. properties and bilayer properties. At the same time, energy conservation is kept within reasonable bounds. We conclude that the newly proposed straight cutoff approach is a viable alternative to the std. shifted potentials used in Martini, offering significant speedup even in the absence of GPUs.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2MXhs1Krtb%252FJ&md5=9bb06db6dbe3171eb07bf7509105fc67
  12. 12
    Abraham, M.; Alekseenko, A.; Bergh, C.; Blau, C.; Briand, E.; Doijade, M.; Fleischmann, S.; Gapsys, V.; Garg, G.; Gorelov, S.; Gouaillardet, G.; Gray, A.; Irrgang, M. E.; Jalalypour, F.; Jordan, J.; Junghans, C.; Kanduri, P.; Keller, S.; Kutzner, C.; Lemkul, J. A.; Lundborg, M.; Merz, P.; Miletic, V.; Morozov, D.; Páll, S.; Schulz, R.; Shirts, M.; Shvetsov, A.; Soproni, B.; van der Spoel, D.; Turner, P.; Uphoff, C.; Villa, A.; Wingbermühle, S.; Zhmurov, A.; Bauer, P.; Hess, B.; Lindahl, E. GROMACS 2023.1 Manual , 2023 DOI: 10.5281/zenodo.7852189 .
    Google Scholar
    There is no corresponding record for this reference.
  13. 13
    Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684– 3690,  DOI: 10.1063/1.448118
    Google Scholar
    13
    Molecular dynamics with coupling to an external bath
    Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; DiNola, A.; Haak, J. R.
    Journal of Chemical Physics (1984), 81 (8), 3684-90CODEN: JCPSA6; ISSN:0021-9606.
    In mol. dynamics (MD) simulations, the need often arises to maintain such parameters as temp. or pressure rather than energy and vol., or to impose gradients for studying transport properties in nonequil. MD. A method is described to realize coupling to an external bath with const. temp. or pressure with adjustable time consts. for the coupling. The method is easily extendable to other variables and to gradients, and can be applied also to polyat. mols. involving internal constraints. The influence of coupling time consts. on dynamical variables is evaluated. A leap-frog algorithm is presented for the general case involving constraints with coupling to both a const. temp. and a const. pressure bath.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaL2cXmtlGksbY%253D&md5=5510dc00297d63b91ee3a7a4ae5aacb1
  14. 14
    Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182– 7190,  DOI: 10.1063/1.328693
    Google Scholar
    14
    Polymorphic transitions in single crystals: A new molecular dynamics method
    Parrinello, M.; Rahman, A.
    Journal of Applied Physics (1981), 52 (12), 7182-90CODEN: JAPIAU; ISSN:0021-8979.
    A Lagrangian formulation is introduced; it can be used to make mol. dynamics (MD) calcns. on systems under the most general, externally applied, conditions of stress. In this formulation the MD cell shape and size can change according to dynamic equations given by this Lagrangian. This MD technique was used to the study of structural transitions of a Ni single crystal under uniform uniaxial compressive and tensile loads. Some results regarding the stress-strain relation obtained by static calcns. are invalid at finite temp. Under compressive loading, the model of Ni shows a bifurcation in its stress-strain relation; this bifurcation provides a link in configuration space between cubic and hexagonal close packing. Such a transition could perhaps be obsd. exptl. under extreme conditions of shock.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaL38XislSnuw%253D%253D&md5=a0a5617389f6cabbf2a405c649aadf03
  15. 15
    Ingólfsson, H. I.; Melo, M. N.; van Eerden, F. J.; Arnarez, C.; Lopez, C. A.; Wassenaar, T. A.; Periole, X.; de Vries, A. H.; Tieleman, D. P.; Marrink, S. J. Lipid Organization of the Plasma Membrane. J. Am. Chem. Soc. 2014, 136, 14554– 14559,  DOI: 10.1021/ja507832e
    Google Scholar
    15
    Lipid Organization of the Plasma Membrane
    Ingolfsson, Helgi I.; Melo, Manuel N.; van Eerden, Floris J.; Arnarez, Clement; Lopez, Cesar A.; Wassenaar, Tsjerk A.; Periole, Xavier; de Vries, Alex H.; Tieleman, D. Peter; Marrink, Siewert J.
    Journal of the American Chemical Society (2014), 136 (41), 14554-14559CODEN: JACSAT; ISSN:0002-7863. (American Chemical Society)
    The detailed organization of cellular membranes remains rather elusive. Based on large-scale mol. dynamics simulations, we provide a high-resoln. view of the lipid organization of a plasma membrane at an unprecedented level of complexity. Our plasma membrane model consists of 63 different lipid species, combining 14 types of headgroups and 11 types of tails asym. distributed across the two leaflets, closely mimicking an idealized mammalian plasma membrane. We observe an enrichment of cholesterol in the outer leaflet and a general non-ideal lateral mixing of the different lipid species. Transient domains with liq.-ordered character form and disappear on the microsecond time scale. These domains are coupled across the two membrane leaflets. In the outer leaflet, distinct nanodomains consisting of gangliosides are obsd. Phosphoinositides show preferential clustering in the inner leaflet. Our data provide a key view on the lateral organization of lipids in one of life's fundamental structures, the cell membrane.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2cXhsFKnsbjI&md5=3fcecc9c57a02ac5a43b899401e22467
  16. 16
    Larsen, A. H. Molecular Dynamics Simulations of Curved Lipid Membranes. Int. J. Mol. Sci. 2022, 23, 8098  DOI: 10.3390/ijms23158098
    Google Scholar
    16
    Molecular Dynamics Simulations of Curved Lipid Membranes
    Larsen, Andreas Haahr
    International Journal of Molecular Sciences (2022), 23 (15), 8098CODEN: IJMCFK; ISSN:1422-0067. (MDPI AG)
    A review. Eukaryotic cells contain membranes with various curvatures, from the near-plane plasma membrane to the highly curved membranes of organelles, vesicles, and membrane protrusions. These curvatures are generated and sustained by curvature-inducing proteins, peptides, and lipids, and describing these mechanisms is an important scientific challenge. In addn. to that, some mols. can sense membrane curvature and thereby be trafficked to specific locations. The description of curvature sensing is another fundamental challenge. Curved lipid membranes and their interplay with membrane-assocd. proteins can be investigated with mol. dynamics (MD) simulations. Various methods for simulating curved membranes with MD are discussed here, including tools for setting up simulation of vesicles and methods for sustaining membrane curvature. The latter are divided into methods that exploit scaffolding virtual beads, methods that use curvature-inducing mols., and methods applying virtual forces. The variety of simulation tools allow researcher to closely match the conditions of exptl. studies of membrane curvatures.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB38XitFCgsLfN&md5=172716d0c5a1b96eb89185053301b6f4
  17. 17
    Vögele, M.; Köfinger, J.; Hummer, G. Hydrodynamics of Diffusion in Lipid Membrane Simulations. Phys. Rev. Lett. 2018, 120, 268104  DOI: 10.1103/PhysRevLett.120.268104
    Google Scholar
    17
    Hydrodynamics of Diffusion in Lipid Membrane Simulations
    Vogele Martin; Kofinger Jurgen; Hummer Gerhard; Hummer Gerhard
    Physical review letters (2018), 120 (26), 268104 ISSN:.
    By performing molecular dynamics simulations with up to 132 million coarse-grained particles in half-micron sized boxes, we show that hydrodynamics quantitatively explains the finite-size effects on diffusion of lipids, proteins, and carbon nanotubes in membranes. The resulting Oseen correction allows us to extract infinite-system diffusion coefficients and membrane surface viscosities from membrane simulations despite the logarithmic divergence of apparent diffusivities with increasing box width. The hydrodynamic theory of diffusion applies also to membranes with asymmetric leaflets and embedded proteins, and to a complex plasma-membrane mimetic.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BB3c%252FnvVKitA%253D%253D&md5=7d6c5e2e9f9f4b20057c2463cac64f8d
  18. 18
    Duboué-Dijon, E.; Hénin, J. Building intuition for binding free energy calculations: Bound state definition, restraints, and symmetry. J. Chem. Phys. 2021, 154, 204101  DOI: 10.1063/5.0057845
    Google Scholar
    18
    Building intuition for binding free energy calculations: Bound state definition, restraints, and symmetry
    Duboue-Dijon, E.; Henin, J.
    Journal of Chemical Physics (2021), 154 (20), 204101CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
    The theory behind computation of abs. binding free energies using explicit-solvent mol. simulations is well-established, yet somewhat complex, with counter-intuitive aspects. This leads to frequent frustration, common misconceptions, and sometimes erroneous numerical treatment. To improve this, we present the main practically relevant segments of the theory with const. ref. to phys. intuition. We pinpoint the role of the implicit or explicit definition of the bound state (or the binding site) to make a robust link between an exptl. measurement and a computational result. We clarify the role of symmetry and discuss cases where symmetry no. corrections have been misinterpreted. In particular, we argue that symmetry corrections as classically presented are a source of confusion and could be advantageously replaced by restraint free energy contributions. We establish that contrary to a common intuition, partial or missing sampling of some modes of sym. bound states does not affect the calcd. decoupling free energies. Finally, we review these questions and pitfalls in the context of a few common practical situations: binding to a sym. receptor (equiv. binding sites), binding of a sym. ligand (equiv. poses), and formation of a sym. complex, in the case of homodimerization. (c) 2021 American Institute of Physics.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3MXhtFGltb%252FN&md5=33d019450c8d2d0e5ad33f16477a7d4a
  19. 19
    Mori, T.; Miyashita, N.; Im, W.; Feig, M.; Sugita, Y. Molecular dynamics simulations of biological membranes and membrane proteins using enhanced conformational sampling algorithms. Biochim. Biophys. Acta., Biomembr. 2016, 1858, 1635– 1651,  DOI: 10.1016/j.bbamem.2015.12.032
    Google Scholar
    19
    Molecular dynamics simulations of biological membranes and membrane proteins using enhanced conformational sampling algorithms
    Mori, Takaharu; Miyashita, Naoyuki; Im, Wonpil; Feig, Michael; Sugita, Yuji
    Biochimica et Biophysica Acta, Biomembranes (2016), 1858 (7_Part_B), 1635-1651CODEN: BBBMBS; ISSN:0005-2736. (Elsevier B.V.)
    A review. This paper reviews various enhanced conformational sampling methods and explicit/implicit solvent/membrane models, as well as their recent applications to the exploration of the structure and dynamics of membranes and membrane proteins. Mol. dynamics simulations have become an essential tool to investigate biol. problems, and their success relies on proper mol. models together with efficient conformational sampling methods. The implicit representation of solvent/membrane environments is reasonable approxn. to the explicit all-atom models, considering the balance between computational cost and simulation accuracy. Implicit models can be easily combined with replica-exchange mol. dynamics methods to explore a wider conformational space of a protein. Other mol. models and enhanced conformational sampling methods are also briefly discussed. As application examples, we introduce recent simulation studies of glycophorin A, phospholamban, amyloid precursor protein, and mixed lipid bilayers and discuss the accuracy and efficiency of each simulation model and method. This article is part of a Special Issue entitled: Membrane Proteins. Guest Editors: J.C. Gumbart and Sergei Noskov.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC28XjvFygtg%253D%253D&md5=02a175fae7d42d40bc672a1e5509b55b
  20. 20
    Kolossváry, I.; Sherman, W. Comprehensive Approach to Simulating Large Scale Conformational Changes in Biological Systems Utilizing a Path Collective Variable and New Barrier Restraint. J. Phys. Chem. B 2023, 127, 5214– 5229,  DOI: 10.1021/acs.jpcb.3c02028
    Google Scholar
    20
    Comprehensive Approach to Simulating Large Scale Conformational Changes in Biological Systems Utilizing a Path Collective Variable and New Barrier Restraint
    Kolossvary, Istvan; Sherman, Woody
    Journal of Physical Chemistry B (2023), 127 (23), 5214-5229CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
    Conformational sampling of complex biomols. is an emerging frontier in drug discovery. Advances in lab-based structural biol. and related computational approaches like AlphaFold have made great strides in obtaining static protein structures for biol. relevant targets. However, biol. is in const. motion, and many important biol. processes rely on conformationally driven events. Conventional mol. dynamics (MD) simulations run on std. hardware are impractical for many drug design projects, where conformationally driven biol. events can take microseconds to milliseconds or longer. An alternative approach is to focus the search on a limited region of conformational space defined by a putative reaction coordinate (i.e., path collective variable). The search space is typically limited by applying restraints, which can be guided by insights about the underlying biol. process of interest. The challenge is striking a balance between the degree to which the system is constrained and still allowing for natural motions along the path. A plethora of restraints exist to limit the size of conformational search space, although each has drawbacks when simulating complex biol. motions. In this work, we present a three-stage procedure to construct realistic path collective variables (PCVs) and introduce a new kind of barrier restraint that is particularly well suited for complex conformationally driven biol. events, such as allosteric modulations and conformational signaling. The PCV presented here is all-atom (as opposed to C-alpha or backbone only) and is derived from all-atom MD trajectory frames. The new restraint relies on a barrier function (specifically, the scaled reciprocal function), which we show is particularly beneficial in the context of mol. dynamics, where near-hard-wall restraints are needed with zero tolerance to restraint violation. We have implemented our PCV and barrier restraint within a hybrid sampling framework that combines well-tempered metadynamics and extended-Lagrangian adaptive biasing force (meta-eABF). We use three particular examples of high pharmaceutical interest to demonstrate the value of this approach: (1) sampling the distance from ubiquitin to a protein of interest within the supramol. cullin-RING ligase complex, (2) stabilizing the wild-type conformation of the oncogenic mutant JAK2-V617F pseudokinase domain, and (3) inducing an activated state of the stimulator of interferon genes (STING) protein obsd. upon ligand binding. For examples 2 and 3, we present statistical anal. of meta-eABF free energy ests. and, for each case, code for reproducing this work.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3sXhtFKjtLzJ&md5=c4e43cfe37b8575041cd6e06467f2c05
  21. 21
    Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604– 613,  DOI: 10.1016/j.cpc.2013.09.018
    Google Scholar
    21
    PLUMED 2: New feathers for an old bird
    Tribello, Gareth A.; Bonomi, Massimiliano; Branduardi, Davide; Camilloni, Carlo; Bussi, Giovanni
    Computer Physics Communications (2014), 185 (2), 604-613CODEN: CPHCBZ; ISSN:0010-4655. (Elsevier B.V.)
    Enhancing sampling and analyzing simulations are central issues in mol. simulation. Recently, we introduced PLUMED, an open-source plug-in that provides some of the most popular mol. dynamics (MD) codes with implementations of a variety of different enhanced sampling algorithms and collective variables (CVs). The rapid changes in this field, in particular new directions in enhanced sampling and dimensionality redn. together with new hardware, require a code that is more flexible and more efficient. We therefore present PLUMED 2 here-a complete rewrite of the code in an object-oriented programming language (C++). This new version introduces greater flexibility and greater modularity, which both extends its core capabilities and makes it far easier to add new methods and CVs. It also has a simpler interface with the MD engines and provides a single software library contg. both tools and core facilities. Ultimately, the new code better serves the ever-growing community of users and contributors in coping with the new challenges arising in the field.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC3sXhs1yqs7fJ&md5=292009aab558d0ef1108bb9a5f036c40
  22. 22
    Frenkel, D.; Smit, B. Introduction. In Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, 2002; Vol. 1.
    Google Scholar
    There is no corresponding record for this reference.
  23. 23
    Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids, 2nd ed.; Oxford University Press, 2017.
    Google Scholar
    There is no corresponding record for this reference.
  24. 24
    Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577– 8593,  DOI: 10.1063/1.470117
    Google Scholar
    24
    A smooth particle mesh Ewald method
    Essmann, Ulrich; Perera, Lalith; Berkowitz, Max L.; Darden, Tom; Lee, Hsing; Pedersen, Lee G.
    Journal of Chemical Physics (1995), 103 (19), 8577-93CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
    The previously developed particle mesh Ewald method is reformulated in terms of efficient B-spline interpolation of the structure factors. This reformulation allows a natural extension of the method to potentials of the form 1/rp with p ≥ 1. Furthermore, efficient calcn. of the virial tensor follows. Use of B-splines in the place of Lagrange interpolation leads to analytic gradients as well as a significant improvement in the accuracy. The authors demonstrate that arbitrary accuracy can be achieved, independent of system size N, at a cost that scales as N log(N). For biomol. systems with many thousands of atoms and this method permits the use of Ewald summation at a computational cost comparable to that of a simple truncation method of 10 Å or less.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK2MXptlehtrw%253D&md5=092a679dd3bee08da28df41e302383a7
  25. 25
    Sega, M.; Dellago, C. Long-range dispersion effects on the water/vapor interface simulated using the most common models. J. Phys. Chem. B 2017, 121, 3798– 3803,  DOI: 10.1021/acs.jpcb.6b12437
    Google Scholar
    25
    Long-Range Dispersion Effects on the Water/Vapor Interface Simulated Using the Most Common Models
    Sega, Marcello; Dellago, Christoph
    Journal of Physical Chemistry B (2017), 121 (15), 3798-3803CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
    The long-range contribution to dispersion forces is known to have a major impact on the properties of inhomogeneous fluids, and its correct treatment is increasingly recognized as being a necessary requirement to avoid cutoff-related artifacts. Although anal. corrections for quantities like the surface tension are known, these cannot take into account the structural changes induced by the long-range contributions. The authors analyze the interfacial properties of seven popular water models, comparing the results with the cutoff version of the dispersion potential. The differences in surface tension ests. are in all cases found to be less than 2 mN/m.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2sXjtVait7o%253D&md5=a32c3a9e14c07687bb492635c4b5c74d
  26. 26
    Panigrahy, R. An Improved Algorithm Finding Nearest Neighbor Using Kd-trees. In Lecture Notes in Computer Science; Springer: Berlin Heidelberg pp 387– 398.
    Google Scholar
    There is no corresponding record for this reference.
  27. 27
    Páll, S.; Zhmurov, A.; Bauer, P.; Abraham, M.; Lundborg, M.; Gray, A.; Hess, B.; Lindahl, E. Heterogeneous parallelization and acceleration of molecular dynamics simulations in GROMACS. J. Chem. Phys. 2020, 153, 134110  DOI: 10.1063/5.0018516
    Google Scholar
    27
    Heterogeneous parallelization and acceleration of molecular dynamics simulations in GROMACS
    Pall, Szilard; Zhmurov, Artem; Bauer, Paul; Abraham, Mark; Lundborg, Magnus; Gray, Alan; Hess, Berk; Lindahl, Erik
    Journal of Chemical Physics (2020), 153 (13), 134110CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
    The introduction of accelerator devices such as graphics processing units (GPUs) has had profound impact on mol. dynamics simulations and has enabled order-of-magnitude performance advances using commodity hardware. To fully reap these benefits, it has been necessary to reformulate some of the most fundamental algorithms, including the Verlet list, pair searching, and cutoffs. Here, we present the heterogeneous parallelization and acceleration design of mol. dynamics implemented in the GROMACS codebase over the last decade. The setup involves a general cluster-based approach to pair lists and non-bonded pair interactions that utilizes both GPU and central processing unit (CPU) single instruction, multiple data acceleration efficiently, including the ability to load-balance tasks between CPUs and GPUs. The algorithm work efficiency is tuned for each type of hardware, and to use accelerators more efficiently, we introduce dual pair lists with rolling pruning updates. Combined with new direct GPU-GPU communication and GPU integration, this enables excellent performance from single GPU simulations through strong scaling across multiple GPUs and efficient multi-node parallelization. (c) 2020 American Institute of Physics.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3cXitVWrtLnE&md5=6997d6a941338d496269c295828c22f9
  28. 28
    Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; in ’t Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D.; Shan, R.; Stevens, M. J.; Tranchida, J.; Trott, C.; Plimpton, S. J. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comput. Phys. Commun. 2022, 271, 108171  DOI: 10.1016/j.cpc.2021.108171
    Google Scholar
    28
    LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
    Thompson, Aidan P.; Aktulga, H. Metin; Berger, Richard; Bolintineanu, Dan S.; Brown, W. Michael; Crozier, Paul S.; in 't Veld, Pieter J.; Kohlmeyer, Axel; Moore, Stan G.; Nguyen, Trung Dac; Shan, Ray; Stevens, Mark J.; Tranchida, Julien; Trott, Christian; Plimpton, Steven J.
    Computer Physics Communications (2022), 271 (), 108171CODEN: CPHCBZ; ISSN:0010-4655. (Elsevier B.V.)
    Since the classical mol. dynamics simulator LAMMPS was released as an open source code in 2004, it has become a widely-used tool for particle-based modeling of materials at length scales ranging from at. to mesoscale to continuum. Reasons for its popularity are that it provides a wide variety of particle interaction models for different materials, that it runs on any platform from a single CPU core to the largest supercomputers with accelerators, and that it gives users control over simulation details, either via the input script or by adding code for new interat. potentials, constraints, diagnostics, or other features needed for their models. As a result, hundreds of people have contributed new capabilities to LAMMPS and it has grown from fifty thousand lines of code in 2004 to a million lines today. In this paper several of the fundamental algorithms used in LAMMPS are described along with the design strategies which have made it flexible for both users and developers. We also highlight some capabilities recently added to the code which were enabled by this flexibility, including dynamic load balancing, on-the-fly visualization, magnetic spin dynamics models, and quantum-accuracy machine learning interat. potentials.Program Title: Large-scale Atomic/Mol. Massively Parallel Simulator (LAMMPS)CPC Library link to program files:https://doi.org/10.17632/cxbxs9btsv.1Developer's repository link:https://github.com/lammps/lammpsLicensing provisions: GPLv2Programming language: C++, Python, C, FortranSupplementary material:https://www.lammps.orgNature of problem: Many science applications in physics, chem., materials science, and related fields require parallel, scalable, and efficient generation of long, stable classical particle dynamics trajectories. Within this common problem definition, there lies a great diversity of use cases, distinguished by different particle interaction models, external constraints, as well as timescales and lengthscales ranging from at. to mesoscale to macroscopic.Soln. method: The LAMMPS code uses parallel spatial decompn., distributed neighbor lists, and parallel FFTs for long-range Coulombic interactions [1]. The time integration algorithm is based on the Stormer-Verlet symplectic integrator [2], which provides better stability than higher-order non-symplectic methods. In addn., LAMMPS supports a wide range of interat. potentials, constraints, diagnostics, software interfaces, and pre- and post-processing features.Addnl. comments including restrictions and unusual features: This paper serves as the definitive ref. for the LAMMPS code.S. Plimpton, Fast parallel algorithms for short-range mol. dynamics. Phys. 117 (1995) 1-19.L. Verlet, Computer expts. on classical fluids: I. Thermodynamical properties of Lennard-Jones mols., Phys. Rev. 159 (1967) 98-103.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3MXitlSrsb7O&md5=cd0bfd050820e97c11779003add20ed3
  29. 29
    Helfrich, W. Elastic Properties of Lipid Bilayers: Theory and Possible Experiments. Z. Naturforsch. C 1973, 28, 693– 703,  DOI: 10.1515/znc-1973-11-1209
    Google Scholar
    29
    Elastic properties of lipid bilayers. Theory and possible experiments
    Helfrich, W.
    Zeitschrift fuer Naturforschung, Teil C: Biochemie, Biophysik, Biologie, Virologie (1973), 28 (11-12), 693-703CODEN: ZNFCAP; ISSN:0341-0471.
    A theory for the elastic behavior of lipid bilayers is presented involving 3 types of strains, i.e. stretching, tilt, and curvature. Possible expts. to det. the elastic properties of lipid bilayers by magnetic and pressure deformation measurements are discussed.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaE2cXnvFersA%253D%253D&md5=c64c8c2c43d3defe6b78d2c06abe0630
  30. 30
    Bhaskara, R. M.; Grumati, P.; Garcia-Pardo, J.; Kalayil, S.; Covarrubias-Pinto, A.; Chen, W.; Kudryashev, M.; Dikic, I.; Hummer, G. Curvature induction and membrane remodeling by FAM134B reticulon homology domain assist selective ER-phagy. Nat. Commun. 2019, 10, 2370  DOI: 10.1038/s41467-019-10345-3
    Google Scholar
    30
    Curvature induction and membrane remodeling by FAM134B reticulon homology domain assist selective ER-phagy
    Bhaskara Ramachandra M; Hummer Gerhard; Grumati Paolo; Covarrubias-Pinto Adriana; Dikic Ivan; Garcia-Pardo Javier; Kalayil Sissy; Chen Wenbo; Kudryashev Mikhail; Dikic Ivan; Garcia-Pardo Javier; Chen Wenbo; Kudryashev Mikhail; Hummer Gerhard
    Nature communications (2019), 10 (1), 2370 ISSN:.
    FAM134B/RETREG1 is a selective ER-phagy receptor that regulates the size and shape of the endoplasmic reticulum. The structure of its reticulon-homology domain (RHD), an element shared with other ER-shaping proteins, and the mechanism of membrane shaping remain poorly understood. Using molecular modeling and molecular dynamics (MD) simulations, we assemble a structural model for the RHD of FAM134B. Through MD simulations of FAM134B in flat and curved membranes, we relate the dynamic RHD structure with its two wedge-shaped transmembrane helical hairpins and two amphipathic helices to FAM134B functions in membrane-curvature induction and curvature-mediated protein sorting. FAM134B clustering, as expected to occur in autophagic puncta, amplifies the membrane-shaping effects. Electron microscopy of in vitro liposome remodeling experiments support the membrane remodeling functions of the different RHD structural elements. Disruption of the RHD structure affects selective autophagy flux and leads to disease states.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BB3M3hs1ahug%253D%253D&md5=10b720cc239064d028e0436dc9e1ae4d
  31. 31
    Wassenaar, T. A.; Ingólfsson, H. I.; Böckmann, R. A.; Tieleman, D. P.; Marrink, S. J. Computational Lipidomics with insane: A Versatile Tool for Generating Custom Membranes for Molecular Simulations. J. Chem. Theory Comput. 2015, 11, 2144– 2155,  DOI: 10.1021/acs.jctc.5b00209
    Google Scholar
    31
    Computational Lipidomics with insane: A Versatile Tool for Generating Custom Membranes for Molecular Simulations
    Wassenaar, Tsjerk A.; Ingolfsson, Helgi I.; Boeckmann, Rainer A.; Tieleman, D. Peter; Marrink, Siewert J.
    Journal of Chemical Theory and Computation (2015), 11 (5), 2144-2155CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
    For simulations of membranes and membrane proteins, the generation of the lipid bilayer is a crit. step in the setup of the system. Membranes comprising multiple components pose a particular challenge, because the relative abundances need to be controlled and the equilibration of the system may take several microseconds. Here the authors present a comprehensive method for building membrane contg. systems, characterized by simplicity and versatility. The program uses preset, coarse-grain lipid templates to build the membrane, and also allows on-the-fly generation of simple lipid types by specifying the headgroup, linker, and lipid tails on the command line. The resulting models can be equilibrated, after which a relaxed atomistic model can be obtained by reverse transformation. For multicomponent membranes, this provides an efficient means for generating equilibrated atomistic models. The method is called insane, an acronym for INSert membrANE. The program has been made available, together with the complementary method for reverse transformation, at http://cgmartini.nl/. This work highlights the key features of insane and presents a survey of properties for a large range of lipids as a start of a computational lipidomics project.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2MXmtlamsbo%253D&md5=8a754c43db80752045ed3bfc3e64b5db
  32. 32
    Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101  DOI: 10.1063/1.2408420
    Google Scholar
    32
    Canonical sampling through velocity rescaling
    Bussi, Giovanni; Donadio, Davide; Parrinello, Michele
    Journal of Chemical Physics (2007), 126 (1), 014101/1-014101/7CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
    The authors present a new mol. dynamics algorithm for sampling the canonical distribution. In this approach the velocities of all the particles are rescaled by a properly chosen random factor. The algorithm is formally justified and it is shown that, in spite of its stochastic nature, a quantity can still be defined that remains const. during the evolution. In numerical applications this quantity can be used to measure the accuracy of the sampling. The authors illustrate the properties of this new method on Lennard-Jones and TIP4P water models in the solid and liq. phases. Its performance is excellent and largely independent of the thermostat parameter also with regard to the dynamic properties.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD2sXosVCltg%253D%253D&md5=9c182b57bfc65bca6be23c8c76b4be77
  33. 33
    Bernetti, M.; Bussi, G. Pressure control using stochastic cell rescaling. J. Chem. Phys. 2020, 153, 114107  DOI: 10.1063/5.0020514
    Google Scholar
    33
    Pressure control using stochastic cell rescaling
    Bernetti, Mattia; Bussi, Giovanni
    Journal of Chemical Physics (2020), 153 (11), 114107CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
    Mol. dynamics simulations require barostats to be performed at a const. pressure. The usual recipe is to employ the Berendsen barostat first, which displays a first-order vol. relaxation efficient in equilibration but results in incorrect vol. fluctuations, followed by a second-order or a Monte Carlo barostat for prodn. runs. In this paper, we introduce stochastic cell rescaling, a first-order barostat that samples the correct vol. fluctuations by including a suitable noise term. The algorithm is shown to report vol. fluctuations compatible with the isobaric ensemble and its anisotropic variant is tested on a membrane simulation. Stochastic cell rescaling can be straightforwardly implemented in the existing codes and can be used effectively in both equilibration and prodn. phases. (c) 2020 American Institute of Physics.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3cXhvVOnt7nN&md5=0021e74a81e9956f6b64873d904a66f2
  34. 34
    Mark, P.; Nilsson, L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A 2001, 105, 9954– 9960,  DOI: 10.1021/jp003020w
    Google Scholar
    34
    Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K
    Mark, Pekka; Nilsson, Lennart
    Journal of Physical Chemistry A (2001), 105 (43), 9954-9960CODEN: JPCAFH; ISSN:1089-5639. (American Chemical Society)
    Mol. dynamics simulations of five water models, the TIP3P (original and modified), SPC (original and refined), and SPC/E (original), were performed using the CHARMM mol. mechanics program. All simulations were carried out in the microcanonical NVE ensemble, using 901 water mols. in a cubic simulation cell furnished with periodic boundary conditions at 298 K. The SHAKE algorithm was used to keep water mols. rigid. Nanosecond trajectories were calcd. with all water models for high statistical accuracy. The characteristic self-diffusion coeffs. D and radial distribution functions, gOO, gOH, and gHH for all five water models were detd. and compared to exptl. data. The effects of velocity rescaling on the self-diffusion coeff. D were examd. All these empirical water models used in this study are similar by having three interaction sites, but the small differences in their pair potentials composed of Lennard-Jones (LJ) and Coulombic terms give significant differences in the calcd. self-diffusion coeffs., and in the height of the second peak of the radial distribution function gOO.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD3MXntlWrurs%253D&md5=fecd3db40210b04e8b2a933ea07b131e
  35. 35
    Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: A web-based graphical user interface for CHARMM. J. Comput. Chem. 2008, 29, 1859– 1865,  DOI: 10.1002/jcc.20945
    Google Scholar
    35
    CHARMM-GUI: a web-based graphical user interface for CHARMM
    Jo, Sunhwan; Kim, Taehoon; Iyer, Vidyashankara G.; Im, Wonpil
    Journal of Computational Chemistry (2008), 29 (11), 1859-1865CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
    CHARMM is an academic research program used widely for macromol. mechanics and dynamics with versatile anal. and manipulation tools of at. coordinates and dynamics trajectories. CHARMM-GUI, http://www.charmm-gui.org, has been developed to provide a web-based graphical user interface to generate various input files and mol. systems to facilitate and standardize the usage of common and advanced simulation techniques in CHARMM. The web environment provides an ideal platform to build and validate a mol. model system in an interactive fashion such that, if a problem is found through visual inspection, one can go back to the previous setup and regenerate the whole system again. In this article, we describe the currently available functional modules of CHARMM-GUI Input Generator that form a basis for the advanced simulation techniques. Future directions of the CHARMM-GUI development project are also discussed briefly together with other features in the CHARMM-GUI website, such as Archive and Movie Gallery.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD1cXosVKksbc%253D&md5=112a3dd61d792b040f9f716b32220d7e
  36. 36
    Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985, 31, 1695– 1697,  DOI: 10.1103/PhysRevA.31.1695
    Google Scholar
    36
    Canonical dynamics: Equilibrium phase-space distributions
    Hoover
    Physical review. A, General physics (1985), 31 (3), 1695-1697 ISSN:0556-2791.
    There is no expanded citation for this reference.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BC2sjotlWltA%253D%253D&md5=99a2477835b37592226a5d18a760685c
  37. 37
    Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511– 519,  DOI: 10.1063/1.447334
    Google Scholar
    37
    A unified formulation of the constant-temperature molecular-dynamics methods
    Nose, Shuichi
    Journal of Chemical Physics (1984), 81 (1), 511-19CODEN: JCPSA6; ISSN:0021-9606.
    Three recently proposed const. temp. mol. dynamics methods [N., (1984) (1); W. G. Hoover et al., (1982) (2); D. J. Evans and G. P. Morris, (1983) (2); and J. M. Haile and S. Gupta, 1983) (3)] are examd. anal. via calcg. the equil. distribution functions and comparing them with that of the canonical ensemble. Except for effects due to momentum and angular momentum conservation, method (1) yields the rigorous canonical distribution in both momentum and coordinate space. Method (2) can be made rigorous in coordinate space, and can be derived from method (1) by imposing a specific constraint. Method (3) is not rigorous and gives a deviation of order N-1/2 from the canonical distribution (N the no. of particles). The results for the const. temp.-const. pressure ensemble are similar to the canonical ensemble case.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaL2cXkvFOrs7k%253D&md5=2974515ec89e5601868e35871c0f19c2
  38. 38
    Welch, P. The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Trans. Audio Electroacoust. 1967, 15, 70– 73,  DOI: 10.1109/TAU.1967.1161901
    Google Scholar
    There is no corresponding record for this reference.
  39. 39
    Hernández-Muñoz, J.; Bresme, F.; Tarazona, P.; Chacón, E. Bending Modulus of Lipid Membranes from Density Correlation Functions. J. Chem. Theory Comput. 2022, 18, 3151– 3163,  DOI: 10.1021/acs.jctc.2c00099
    Google Scholar
    39
    Bending Modulus of Lipid Membranes from Density Correlation Functions
    Hernandez-Munoz, Jose; Bresme, Fernando; Tarazona, Pedro; Chacon, Enrique
    Journal of Chemical Theory and Computation (2022), 18 (5), 3151-3163CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
    The bending modulus κ quantifies the elasticity of biol. membranes in terms of the free energy cost of increasing the membrane corrugation. Mol. dynamics (MD) simulations provide a powerful approach to quantify κ by analyzing the thermal fluctuations of the lipid bilayer. However, existing methods require the identification and filtering of non-mesoscopic fluctuation modes. State of the art methods rely on identifying a smooth surface to describe the membrane shape. These methods introduce uncertainties in calcg. κ since they rely on different criteria to select the relevant fluctuation modes. Here, we present a method to compute κ using mol. simulations. Our approach circumvents the need to define a mesoscopic surface or an orientation field for the lipid tails explicitly. The bending and tilt moduli can be extd. from the anal. of the d. correlation function (DCF). The method introduced here builds on the Bedeaux and Weeks (BW) theory for the DCF of fluctuating interfaces and on the coupled undulatory (CU) mode introduced by us in previous work. We test the BW-DCF method by computing the elastic properties of lipid membranes with different system sizes (from 500 to 6000 lipid mols.) and using coarse-grained (for POPC and DPPC lipids) and fully atomistic models (for DPPC). Further, we quantify the impact of cholesterol on the bending modulus of DPPC bilayers. We compare our results with bending moduli obtained with X-ray diffraction data and different computer simulation methods.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB38XptVOitLs%253D&md5=1ce4b80044739c9c7b78430e2d1aa421
  40. 40
    Eid, J.; Razmazma, H.; Jraij, A.; Ebrahimi, A.; Monticelli, L. On Calculating the Bending Modulus of Lipid Bilayer Membranes from Buckling Simulations. J. Phys. Chem. B 2020, 124, 6299– 6311,  DOI: 10.1021/acs.jpcb.0c04253
    Google Scholar
    40
    On calculating the bending modulus of lipid bilayer membranes from buckling simulations
    Eid, Jad; Razmazma, Hafez; Jraij, Alia; Ebrahimi, Ali; Monticelli, Luca
    Journal of Physical Chemistry B (2020), 124 (29), 6299-6311CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
    The bending modulus is an important phys. const. characterizing lipid membranes. Different methods have been devised for calcg. the bending modulus from simulations, and one of them, named the buckling method, is nowadays widely used due to its simplicity and numerical stability. However, questions remain on the reproducibility, finite size effects, and interpretation of results on lipid mixts. Here we explore the dependence of simulation results on the system size and the strain. We find that the dimensions of the box have a negligible impact on the results when the system size is beyond a certain threshold. We then calc. the bending rigidity for of a series of common single-component lipid bilayers (PC, PS, PE, PG, and SM), as well as a no. of binary and ternary lipid mixts. We find that bending moduli of lipid mixts. can be predicted from the weighted av. of the moduli of the individual components, as long as the mixt. is homogeneous. For phase-sepd. mixts., the apparent elastic modulus is closer to the value of the softer component. Predictions of the bending modulus based on the area compressibility modulus are found to be generally unreliable.
    >> More from SciFinder ®
    https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3cXht1yntbnL&md5=3aa840e0da955031d874d2b047514f6f

Cited By

ARTICLE SECTIONS
Jump To

This article is cited by 1 publications.

  1. Denys Biriukov, Robert Vácha. Pathways to a Shiny Future: Building the Foundation for Computational Physical Chemistry and Biophysics in 2050. ACS Physical Chemistry Au 2024, Article ASAP.
Download PDF

  • Abstract

    High Resolution Image
    Download MS PowerPoint Slide

    Figure 1

    Figure 1. Large Martini POPC bilayer crumples in MD simulations with default simulation parameters, yet stays flat with frequent neighbor list updates. (A) Snapshots of the membrane (phosphate groups in gold) in MD simulation with default parameters. The Verlet-buffer-tolerance, which denotes the maximally allowed energy drift per particle between neighbor list updates due to missed nonbonded interactions, was set to VBT = 0.005 kJ·mol–1·ps–1; the outer cutoff was rl = 1.269 nm; the number of time steps between neighbor list updates was nstlist = 25; and dual pair list was enabled. (B) Snapshots in an MD simulation with neighbor list updates enforced at every time step (nstlist = 1). Snapshots are at time points 50, 150, 250, and 350 ns (left to right). Simulation boxes are indicated as blue lines, and box heights Lz are listed.

    High Resolution Image
    Download MS PowerPoint Slide

    Figure 2

    Figure 2. Pressure tensor elements deviate from the target pressure in MD simulations with the default cutoff handling. Running averages of the diagonal elements of the pressure tensor are shown for (A) the large and (B) the small Martini membrane systems in NPT MD simulations and (C) for the Martini water system in an NVT simulation. The left column shows snapshots of the systems. The center and right columns show the running averages evaluated every nstenergy = 1 and 100 steps, respectively. Results obtained with the default simulation parameters for cutoff handling (VBT = 0.005 kJ·mol–1·ps–1) are shown as dashed lines (see the legend for color). The solid lines show results for a larger outer cutoff rl = 1.422 nm with nstlist = 20 fixed and dual pair list disabled. In (A, B), the target pressure of 1 bar in the NPT simulations is indicated by a dashed black line. For the NVT simulation in (C), the dashed black line indicates the consistent average obtained with rl = 1.422 nm and nstlist = 20.

    High Resolution Image
    Download MS PowerPoint Slide

    Figure 3

    Figure 3. Average pressure tensor components deviate from the target pressure between neighbor list updates. Averages were performed over blocks of nstlist time steps starting immediately after a neighbor list update (time step 0). Results are shown for the large (top) and small (center) membrane systems in NPT simulations and for the NVT Martini water system (bottom). Results include the default cutoff setting with nstlist = 25 (three left columns) and the setting with rl = 1.422 nm and nstlist = 20 (right column). The default rl values for the NPT large and small membranes and the NVT solvent were 1.269, 1.267, and 1.28 nm, respectively. For the runs in column 3, the dual pair list was disabled. We averaged the lateral pressure components P∥ = (Pxx + Pyy)/2 (orange curves) and showed the normal pressure as P⊥ = Pzz (blue curves). Dashed horizontal lines of matching colors denote the corresponding overall averages. Black horizontal dotted lines represent the reference pressure values. Cyan circles (left column) indicate deviations of the averaged pressures from the target. Green circles (column 3) indicate deviations just before the neighbor list update.

    High Resolution Image
    Download MS PowerPoint Slide

    Figure 4

    Figure 4. Differences ΔP in the pressures calculated just before and right after neighbor list updates follow the trend of missed interactions. ΔP (left scale) for lateral (orange) and normal pressures (blue) are shown as functions of rl (A) for Martini water and (B) for TIP3P water, both simulated in NVT conditions with nstlist = 20 and dual pair list disabled. Standard deviations are shown as error bars. Also shown (right scale) is the expected number of missed interactions (black dashed lines). Considering the Martini water particles as point particles, the number of missed interactions (nmissed) is evaluated using eq S8. For TIP3P water, the expected number of the missed hydrogen–hydrogen (nH–H), oxygen–hydrogen (nO–H), and oxygen–oxygen (nO–O) electrostatic interactions are evaluated using eqs S8 and S11.

    High Resolution Image
    Download MS PowerPoint Slide

    Figure 5

    Figure 5. Power spectral density (PSD) of the scalar pressure calculated for Martini water in the NVT ensemble. Results are shown for rl = 1.27 nm and nstlist = 25 (top curves; scaled by factor 100) with dual pair list and inner neighbor list updates every four (orange) and five (blue line) integration steps, respectively. Except for the oscillations at frequency multiples of 1/(nstlist × Δt), the two curves nearly superimpose. Also shown (fluorescent green, bottom curve) is the PSD for a simulation with rl = 1.422 nm and nstlist = 25 without dual cutoff, which does not show oscillations.

    High Resolution Image
    Download MS PowerPoint Slide

  • References

    ARTICLE SECTIONS
    Jump To

    This article references 40 other publications.

    1. 1
      Gelpi, J.; Hospital, A.; Goñi, R.; Orozco, M. Molecular dynamics simulations: Advances and applications. Adv. Appl. Bioinf. Chem. 2015, 8, 37– 47,  DOI: 10.2147/AABC.S70333
      1
      Molecular dynamics simulations: advances and applications
      Hospital Adam; Goni Josep Ramon; Orozco Modesto; Gelpi Josep L
      Advances and applications in bioinformatics and chemistry : AABC (2015), 8 (), 37-47 ISSN:1178-6949.
      Molecular dynamics simulations have evolved into a mature technique that can be used effectively to understand macromolecular structure-to-function relationships. Present simulation times are close to biologically relevant ones. Information gathered about the dynamic properties of macromolecules is rich enough to shift the usual paradigm of structural bioinformatics from studying single structures to analyze conformational ensembles. Here, we describe the foundations of molecular dynamics and the improvements made in the direction of getting such ensemble. Specific application of the technique to three main issues (allosteric regulation, docking, and structure refinement) is discussed.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BC28rgs1ShtQ%253D%253D&md5=1b8c8639bd3356dfb4d700a3ec325729
    2. 2
      Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1–2, 19– 25,  DOI: 10.1016/j.softx.2015.06.001
      There is no corresponding record for this reference.
    3. 3
      Lundborg, M.; Lindahl, E. Automatic GROMACS Topology Generation and Comparisons of Force Fields for Solvation Free Energy Calculations. J. Phys. Chem. B 2015, 119, 810– 823,  DOI: 10.1021/jp505332p
      3
      Automatic GROMACS Topology Generation and Comparisons of Force Fields for Solvation Free Energy Calculations
      Lundborg, Magnus; Lindahl, Erik
      Journal of Physical Chemistry B (2015), 119 (3), 810-823CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
      Free energy calcn. has long been an important goal for mol. dynamics simulation and force field development, but historically it has been challenged by limited performance, accuracy, and creation of topologies for arbitrary small mols. This has made it difficult to systematically compare different sets of parameters to improve existing force fields, but in the past few years several authors have developed increasingly automated procedures to generate parameters for force fields such as Amber, CHARMM, and OPLS. The authors present a new framework that enables fully automated generation of GROMACS topologies for any of these force fields and an automated setup for parallel adaptive optimization of high-throughput free energy calcn. by adjusting lambda point placement on the fly. As a small example of this automated pipeline, the authors have calcd. solvation free energies of 50 different small mols. using the GAFF, OPLS-AA, and CGenFF force fields and four different water models, and by including the often neglected polarization costs, the authors show that the common charge models are somewhat under-polarized.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2cXhvVWns7zF&md5=3be571dedfcae41667937c7108217be4
    4. 4
      Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812– 7824,  DOI: 10.1021/jp071097f
      4
      The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations
      Marrink, Siewert J.; Risselada, H. Jelger; Yefimov, Serge; Tieleman, D. Peter; De Vries, Alex H.
      Journal of Physical Chemistry B (2007), 111 (27), 7812-7824CODEN: JPCBFK; ISSN:1520-6106. (American Chemical Society)
      We present an improved and extended version of our coarse grained lipid model. The new version, coined the MARTINI force field, is parametrized in a systematic way, based on the reprodn. of partitioning free energies between polar and apolar phases of a large no. of chem. compds. To reproduce the free energies of these chem. building blocks, the no. of possible interaction levels of the coarse-grained sites has increased compared to those of the previous model. Application of the new model to lipid bilayers shows an improved behavior in terms of the stress profile across the bilayer and the tendency to form pores. An extension of the force field now also allows the simulation of planar (ring) compds., including sterols. Application to a bilayer/cholesterol system at various concns. shows the typical cholesterol condensation effect similar to that obsd. in all atom representations.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD2sXmsVKmsLc%253D&md5=428d5750b94652e4917d905a30658235
    5. 5
      Gecht, M.; Siggel, M.; Linke, M.; Hummer, G.; Köfinger, J. MDBenchmark: A toolkit to optimize the performance of molecular dynamics simulations. J. Chem. Phys. 2020, 153, 144105  DOI: 10.1063/5.0019045
      5
      MDBenchmark: A toolkit to optimize the performance of molecular dynamics simulations
      Gecht, Michael; Siggel, Marc; Linke, Max; Hummer, Gerhard; Koefinger, Juergen
      Journal of Chemical Physics (2020), 153 (14), 144105CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
      Despite the impending flattening of Moores law, the system size, complexity, and length of mol. dynamics (MD) simulations keep on increasing, thanks to effective code parallelization and optimization combined with algorithmic developments. Going forward, exascale computing poses new challenges to the efficient execution and management of MD simulations. The diversity and rapid developments of hardware architectures, software environments, and MD engines make it necessary that users can easily run benchmarks to optimally set up simulations, both with respect to time-to-soln. and overall efficiency. To this end, we have developed the software MDBenchmark to streamline the setup, submission, and anal. of simulation benchmarks and scaling studies. The software design is open and as such not restricted to any specific MD engine or job queuing system. To illustrate the necessity and benefits of running benchmarks and the capabilities of MDBenchmark, we measure the performance of a diverse set of 23 MD simulation systems using GROMACS 2018. We compare the scaling of simulations with the no. of nodes for central processing unit (CPU)-only and mixed CPU-graphics processing unit (GPU) nodes and study the performance that can be achieved when running multiple simulations on a single node. In all these cases, we optimize the nos. of message passing interface (MPI) ranks and open multi-processing (OpenMP) threads, which is crucial to maximizing performance. Our results demonstrate the importance of benchmarking for finding the optimal system and hardware specific simulation parameters. Running MD simulations with optimized settings leads to a significant performance increase that reduces the monetary, energetic, and environmental costs of MD simulations. (c) 2020 American Institute of Physics.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3cXitVCmurvF&md5=d8cd0cbf9ed72f74d83f9170dacd1a11
    6. 6
      Kutzner, C.; Páll, S.; Fechner, M.; Esztermann, A.; de Groot, B. L.; Grubmüller, H. Best bang for your buck: GPU nodes for GROMACS biomolecular simulations. J. Comput. Chem. 2015, 36, 1990– 2008,  DOI: 10.1002/jcc.24030
      6
      Best bang for your buck: GPU nodes for GROMACS biomolecular simulations
      Kutzner, Carsten; Pall, Szilard; Fechner, Martin; Esztermann, Ansgar; de Groot, Bert L.; Grubmueller, Helmut
      Journal of Computational Chemistry (2015), 36 (26), 1990-2008CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
      The mol. dynamics simulation package GROMACS runs efficiently on a wide variety of hardware from commodity workstations to high performance computing clusters. Hardware features are well-exploited with a combination of single instruction multiple data, multithreading, and message passing interface (MPI)-based single program multiple data/multiple program multiple data parallelism while graphics processing units (GPUs) can be used as accelerators to compute interactions off-loaded from the CPU. Here, we evaluate which hardware produces trajectories with GROMACS 4.6 or 5.0 in the most economical way. We have assembled and benchmarked compute nodes with various CPU/GPU combinations to identify optimal compns. in terms of raw trajectory prodn. rate, performance-to-price ratio, energy efficiency, and several other criteria. Although hardware prices are naturally subject to trends and fluctuations, general tendencies are clearly visible. Adding any type of GPU significantly boosts a node's simulation performance. For inexpensive consumer-class GPUs this improvement equally reflects in the performance-to-price ratio. Although memory issues in consumer-class GPUs could pass unnoticed as these cards do not support error checking and correction memory, unreliable GPUs can be sorted out with memory checking tools. Apart from the obvious determinants for cost-efficiency like hardware expenses and raw performance, the energy consumption of a node is a major cost factor. Over the typical hardware lifetime until replacement of a few years, the costs for elec. power and cooling can become larger than the costs of the hardware itself. Taking that into account, nodes with a well-balanced ratio of CPU and consumer-class GPU resources produce the max. amt. of GROMACS trajectory over their lifetime. © 2015 The Authors. Journal of Computational Chem. Published by Wiley Periodicals, Inc.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2MXht1Krur%252FM&md5=2effcd08bab28e1c572dc1dac25e9c01
    7. 7
      Páll, S.; Abraham, M. J.; Kutzner, C.; Hess, B.; Lindahl, E. Tackling Exascale Software Challenges in Molecular Dynamics Simulations with GROMACS. In Lecture Notes in Computer Science; Springer, 2015; pp 3– 27.
      There is no corresponding record for this reference.
    8. 8
      Páll, S.; Hess, B. A flexible algorithm for calculating pair interactions on SIMD architectures. Comput. Phys. Commun. 2013, 184, 2641– 2650,  DOI: 10.1016/j.cpc.2013.06.003
      8
      A flexible algorithm for calculating pair interactions on SIMD architectures
      Pall, Szilard; Hess, Berk
      Computer Physics Communications (2013), 184 (12), 2641-2650CODEN: CPHCBZ; ISSN:0010-4655. (Elsevier B.V.)
      Calcg. interactions or correlations between pairs of particles is typically the most time-consuming task in particle simulation or correlation anal. Straightforward implementations using a double loop over particle pairs have traditionally worked well, esp. since compilers usually do a good job of unrolling the inner loop. In order to reach high performance on modern CPU and accelerator architectures, single-instruction multiple-data (SIMD) parallelization has become essential. Avoiding memory bottlenecks is also increasingly important and requires reducing the ratio of memory to arithmetic operations. Moreover, when pairs only interact within a certain cut-off distance, good SIMD utilization can only be achieved by reordering input and output data, which quickly becomes a limiting factor. Here we present an algorithm for SIMD parallelization based on grouping a fixed no. of particles, e.g. 2, 4, or 8, into spatial clusters. Calcg. all interactions between particles in a pair of such clusters improves data reuse compared to the traditional scheme and results in a more efficient SIMD parallelization. Adjusting the cluster size allows the algorithm to map to SIMD units of various widths. This flexibility not only enables fast and efficient implementation on current CPUs and accelerator architectures like GPUs or Intel MIC, but it also makes the algorithm future-proof. We present the algorithm with an application to mol. dynamics simulations, where we can also make use of the effective buffering the method introduces.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC3sXhtFCht7zM&md5=6ec6327a4dc3fbc551350132aef15485
    9. 9
      Buyya, R.; Vecchiola, C.; Selvi, S. T. Principles of Parallel and Distributed Computing. In Mastering Cloud Computing; Morgan Kaufmann: Boston, 2013; Chapter 2, pp 29– 70.
      There is no corresponding record for this reference.
    10. 10
      Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435– 447,  DOI: 10.1021/ct700301q
      10
      GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation
      Hess, Berk; Kutzner, Carsten; van der Spoel, David; Lindahl, Erik
      Journal of Chemical Theory and Computation (2008), 4 (3), 435-447CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
      Mol. simulation is an extremely useful, but computationally very expensive tool for studies of chem. and biomol. systems. Here, we present a new implementation of our mol. simulation toolkit GROMACS which now both achieves extremely high performance on single processors from algorithmic optimizations and hand-coded routines and simultaneously scales very well on parallel machines. The code encompasses a minimal-communication domain decompn. algorithm, full dynamic load balancing, a state-of-the-art parallel constraint solver, and efficient virtual site algorithms that allow removal of hydrogen atom degrees of freedom to enable integration time steps up to 5 fs for atomistic simulations also in parallel. To improve the scaling properties of the common particle mesh Ewald electrostatics algorithms, we have in addn. used a Multiple-Program, Multiple-Data approach, with sep. node domains responsible for direct and reciprocal space interactions. Not only does this combination of algorithms enable extremely long simulations of large systems but also it provides that simulation performance on quite modest nos. of std. cluster nodes.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD1cXhsVSqurc%253D&md5=d53c94901386260221792ea30f151c5f
    11. 11
      de Jong, D. H.; Baoukina, S.; Ingólfsson, H. I.; Marrink, S. J. Martini straight: Boosting performance using a shorter cutoff and GPUs. Comput. Phys. Commun. 2016, 199, 1– 7,  DOI: 10.1016/j.cpc.2015.09.014
      11
      Martini straight: Boosting performance using a shorter cutoff and GPUs
      de Jong, Djurre H.; Baoukina, Svetlana; Ingolfsson, Helgi I.; Marrink, Siewert J.
      Computer Physics Communications (2016), 199 (), 1-7CODEN: CPHCBZ; ISSN:0010-4655. (Elsevier B.V.)
      In mol. dynamics simulations, sufficient sampling is of key importance and a continuous challenge in the field. The coarse grain Martini force field has been widely used to enhance sampling. In its original implementation, this force field applied a shifted Lennard-Jones potential for the non-bonded van der Waals interactions, to avoid problems related to a relatively short cutoff. Here we investigate the use of a straight cutoff Lennard-Jones potential with potential modifiers. Together with a Verlet neighbor search algorithm, the modified potential allows the use of GPUs to accelerate the computations in Gromacs. We find that this alternative potential has little influence on most of the properties studied, including partitioning free energies, bulk liq. properties and bilayer properties. At the same time, energy conservation is kept within reasonable bounds. We conclude that the newly proposed straight cutoff approach is a viable alternative to the std. shifted potentials used in Martini, offering significant speedup even in the absence of GPUs.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2MXhs1Krtb%252FJ&md5=9bb06db6dbe3171eb07bf7509105fc67
    12. 12
      Abraham, M.; Alekseenko, A.; Bergh, C.; Blau, C.; Briand, E.; Doijade, M.; Fleischmann, S.; Gapsys, V.; Garg, G.; Gorelov, S.; Gouaillardet, G.; Gray, A.; Irrgang, M. E.; Jalalypour, F.; Jordan, J.; Junghans, C.; Kanduri, P.; Keller, S.; Kutzner, C.; Lemkul, J. A.; Lundborg, M.; Merz, P.; Miletic, V.; Morozov, D.; Páll, S.; Schulz, R.; Shirts, M.; Shvetsov, A.; Soproni, B.; van der Spoel, D.; Turner, P.; Uphoff, C.; Villa, A.; Wingbermühle, S.; Zhmurov, A.; Bauer, P.; Hess, B.; Lindahl, E. GROMACS 2023.1 Manual , 2023 DOI: 10.5281/zenodo.7852189 .
      There is no corresponding record for this reference.
    13. 13
      Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684– 3690,  DOI: 10.1063/1.448118
      13
      Molecular dynamics with coupling to an external bath
      Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; DiNola, A.; Haak, J. R.
      Journal of Chemical Physics (1984), 81 (8), 3684-90CODEN: JCPSA6; ISSN:0021-9606.
      In mol. dynamics (MD) simulations, the need often arises to maintain such parameters as temp. or pressure rather than energy and vol., or to impose gradients for studying transport properties in nonequil. MD. A method is described to realize coupling to an external bath with const. temp. or pressure with adjustable time consts. for the coupling. The method is easily extendable to other variables and to gradients, and can be applied also to polyat. mols. involving internal constraints. The influence of coupling time consts. on dynamical variables is evaluated. A leap-frog algorithm is presented for the general case involving constraints with coupling to both a const. temp. and a const. pressure bath.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaL2cXmtlGksbY%253D&md5=5510dc00297d63b91ee3a7a4ae5aacb1
    14. 14
      Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182– 7190,  DOI: 10.1063/1.328693
      14
      Polymorphic transitions in single crystals: A new molecular dynamics method
      Parrinello, M.; Rahman, A.
      Journal of Applied Physics (1981), 52 (12), 7182-90CODEN: JAPIAU; ISSN:0021-8979.
      A Lagrangian formulation is introduced; it can be used to make mol. dynamics (MD) calcns. on systems under the most general, externally applied, conditions of stress. In this formulation the MD cell shape and size can change according to dynamic equations given by this Lagrangian. This MD technique was used to the study of structural transitions of a Ni single crystal under uniform uniaxial compressive and tensile loads. Some results regarding the stress-strain relation obtained by static calcns. are invalid at finite temp. Under compressive loading, the model of Ni shows a bifurcation in its stress-strain relation; this bifurcation provides a link in configuration space between cubic and hexagonal close packing. Such a transition could perhaps be obsd. exptl. under extreme conditions of shock.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaL38XislSnuw%253D%253D&md5=a0a5617389f6cabbf2a405c649aadf03
    15. 15
      Ingólfsson, H. I.; Melo, M. N.; van Eerden, F. J.; Arnarez, C.; Lopez, C. A.; Wassenaar, T. A.; Periole, X.; de Vries, A. H.; Tieleman, D. P.; Marrink, S. J. Lipid Organization of the Plasma Membrane. J. Am. Chem. Soc. 2014, 136, 14554– 14559,  DOI: 10.1021/ja507832e
      15
      Lipid Organization of the Plasma Membrane
      Ingolfsson, Helgi I.; Melo, Manuel N.; van Eerden, Floris J.; Arnarez, Clement; Lopez, Cesar A.; Wassenaar, Tsjerk A.; Periole, Xavier; de Vries, Alex H.; Tieleman, D. Peter; Marrink, Siewert J.
      Journal of the American Chemical Society (2014), 136 (41), 14554-14559CODEN: JACSAT; ISSN:0002-7863. (American Chemical Society)
      The detailed organization of cellular membranes remains rather elusive. Based on large-scale mol. dynamics simulations, we provide a high-resoln. view of the lipid organization of a plasma membrane at an unprecedented level of complexity. Our plasma membrane model consists of 63 different lipid species, combining 14 types of headgroups and 11 types of tails asym. distributed across the two leaflets, closely mimicking an idealized mammalian plasma membrane. We observe an enrichment of cholesterol in the outer leaflet and a general non-ideal lateral mixing of the different lipid species. Transient domains with liq.-ordered character form and disappear on the microsecond time scale. These domains are coupled across the two membrane leaflets. In the outer leaflet, distinct nanodomains consisting of gangliosides are obsd. Phosphoinositides show preferential clustering in the inner leaflet. Our data provide a key view on the lateral organization of lipids in one of life's fundamental structures, the cell membrane.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2cXhsFKnsbjI&md5=3fcecc9c57a02ac5a43b899401e22467
    16. 16
      Larsen, A. H. Molecular Dynamics Simulations of Curved Lipid Membranes. Int. J. Mol. Sci. 2022, 23, 8098  DOI: 10.3390/ijms23158098
      16
      Molecular Dynamics Simulations of Curved Lipid Membranes
      Larsen, Andreas Haahr
      International Journal of Molecular Sciences (2022), 23 (15), 8098CODEN: IJMCFK; ISSN:1422-0067. (MDPI AG)
      A review. Eukaryotic cells contain membranes with various curvatures, from the near-plane plasma membrane to the highly curved membranes of organelles, vesicles, and membrane protrusions. These curvatures are generated and sustained by curvature-inducing proteins, peptides, and lipids, and describing these mechanisms is an important scientific challenge. In addn. to that, some mols. can sense membrane curvature and thereby be trafficked to specific locations. The description of curvature sensing is another fundamental challenge. Curved lipid membranes and their interplay with membrane-assocd. proteins can be investigated with mol. dynamics (MD) simulations. Various methods for simulating curved membranes with MD are discussed here, including tools for setting up simulation of vesicles and methods for sustaining membrane curvature. The latter are divided into methods that exploit scaffolding virtual beads, methods that use curvature-inducing mols., and methods applying virtual forces. The variety of simulation tools allow researcher to closely match the conditions of exptl. studies of membrane curvatures.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB38XitFCgsLfN&md5=172716d0c5a1b96eb89185053301b6f4
    17. 17
      Vögele, M.; Köfinger, J.; Hummer, G. Hydrodynamics of Diffusion in Lipid Membrane Simulations. Phys. Rev. Lett. 2018, 120, 268104  DOI: 10.1103/PhysRevLett.120.268104
      17
      Hydrodynamics of Diffusion in Lipid Membrane Simulations
      Vogele Martin; Kofinger Jurgen; Hummer Gerhard; Hummer Gerhard
      Physical review letters (2018), 120 (26), 268104 ISSN:.
      By performing molecular dynamics simulations with up to 132 million coarse-grained particles in half-micron sized boxes, we show that hydrodynamics quantitatively explains the finite-size effects on diffusion of lipids, proteins, and carbon nanotubes in membranes. The resulting Oseen correction allows us to extract infinite-system diffusion coefficients and membrane surface viscosities from membrane simulations despite the logarithmic divergence of apparent diffusivities with increasing box width. The hydrodynamic theory of diffusion applies also to membranes with asymmetric leaflets and embedded proteins, and to a complex plasma-membrane mimetic.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BB3c%252FnvVKitA%253D%253D&md5=7d6c5e2e9f9f4b20057c2463cac64f8d
    18. 18
      Duboué-Dijon, E.; Hénin, J. Building intuition for binding free energy calculations: Bound state definition, restraints, and symmetry. J. Chem. Phys. 2021, 154, 204101  DOI: 10.1063/5.0057845
      18
      Building intuition for binding free energy calculations: Bound state definition, restraints, and symmetry
      Duboue-Dijon, E.; Henin, J.
      Journal of Chemical Physics (2021), 154 (20), 204101CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
      The theory behind computation of abs. binding free energies using explicit-solvent mol. simulations is well-established, yet somewhat complex, with counter-intuitive aspects. This leads to frequent frustration, common misconceptions, and sometimes erroneous numerical treatment. To improve this, we present the main practically relevant segments of the theory with const. ref. to phys. intuition. We pinpoint the role of the implicit or explicit definition of the bound state (or the binding site) to make a robust link between an exptl. measurement and a computational result. We clarify the role of symmetry and discuss cases where symmetry no. corrections have been misinterpreted. In particular, we argue that symmetry corrections as classically presented are a source of confusion and could be advantageously replaced by restraint free energy contributions. We establish that contrary to a common intuition, partial or missing sampling of some modes of sym. bound states does not affect the calcd. decoupling free energies. Finally, we review these questions and pitfalls in the context of a few common practical situations: binding to a sym. receptor (equiv. binding sites), binding of a sym. ligand (equiv. poses), and formation of a sym. complex, in the case of homodimerization. (c) 2021 American Institute of Physics.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3MXhtFGltb%252FN&md5=33d019450c8d2d0e5ad33f16477a7d4a
    19. 19
      Mori, T.; Miyashita, N.; Im, W.; Feig, M.; Sugita, Y. Molecular dynamics simulations of biological membranes and membrane proteins using enhanced conformational sampling algorithms. Biochim. Biophys. Acta., Biomembr. 2016, 1858, 1635– 1651,  DOI: 10.1016/j.bbamem.2015.12.032
      19
      Molecular dynamics simulations of biological membranes and membrane proteins using enhanced conformational sampling algorithms
      Mori, Takaharu; Miyashita, Naoyuki; Im, Wonpil; Feig, Michael; Sugita, Yuji
      Biochimica et Biophysica Acta, Biomembranes (2016), 1858 (7_Part_B), 1635-1651CODEN: BBBMBS; ISSN:0005-2736. (Elsevier B.V.)
      A review. This paper reviews various enhanced conformational sampling methods and explicit/implicit solvent/membrane models, as well as their recent applications to the exploration of the structure and dynamics of membranes and membrane proteins. Mol. dynamics simulations have become an essential tool to investigate biol. problems, and their success relies on proper mol. models together with efficient conformational sampling methods. The implicit representation of solvent/membrane environments is reasonable approxn. to the explicit all-atom models, considering the balance between computational cost and simulation accuracy. Implicit models can be easily combined with replica-exchange mol. dynamics methods to explore a wider conformational space of a protein. Other mol. models and enhanced conformational sampling methods are also briefly discussed. As application examples, we introduce recent simulation studies of glycophorin A, phospholamban, amyloid precursor protein, and mixed lipid bilayers and discuss the accuracy and efficiency of each simulation model and method. This article is part of a Special Issue entitled: Membrane Proteins. Guest Editors: J.C. Gumbart and Sergei Noskov.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC28XjvFygtg%253D%253D&md5=02a175fae7d42d40bc672a1e5509b55b
    20. 20
      Kolossváry, I.; Sherman, W. Comprehensive Approach to Simulating Large Scale Conformational Changes in Biological Systems Utilizing a Path Collective Variable and New Barrier Restraint. J. Phys. Chem. B 2023, 127, 5214– 5229,  DOI: 10.1021/acs.jpcb.3c02028
      20
      Comprehensive Approach to Simulating Large Scale Conformational Changes in Biological Systems Utilizing a Path Collective Variable and New Barrier Restraint
      Kolossvary, Istvan; Sherman, Woody
      Journal of Physical Chemistry B (2023), 127 (23), 5214-5229CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
      Conformational sampling of complex biomols. is an emerging frontier in drug discovery. Advances in lab-based structural biol. and related computational approaches like AlphaFold have made great strides in obtaining static protein structures for biol. relevant targets. However, biol. is in const. motion, and many important biol. processes rely on conformationally driven events. Conventional mol. dynamics (MD) simulations run on std. hardware are impractical for many drug design projects, where conformationally driven biol. events can take microseconds to milliseconds or longer. An alternative approach is to focus the search on a limited region of conformational space defined by a putative reaction coordinate (i.e., path collective variable). The search space is typically limited by applying restraints, which can be guided by insights about the underlying biol. process of interest. The challenge is striking a balance between the degree to which the system is constrained and still allowing for natural motions along the path. A plethora of restraints exist to limit the size of conformational search space, although each has drawbacks when simulating complex biol. motions. In this work, we present a three-stage procedure to construct realistic path collective variables (PCVs) and introduce a new kind of barrier restraint that is particularly well suited for complex conformationally driven biol. events, such as allosteric modulations and conformational signaling. The PCV presented here is all-atom (as opposed to C-alpha or backbone only) and is derived from all-atom MD trajectory frames. The new restraint relies on a barrier function (specifically, the scaled reciprocal function), which we show is particularly beneficial in the context of mol. dynamics, where near-hard-wall restraints are needed with zero tolerance to restraint violation. We have implemented our PCV and barrier restraint within a hybrid sampling framework that combines well-tempered metadynamics and extended-Lagrangian adaptive biasing force (meta-eABF). We use three particular examples of high pharmaceutical interest to demonstrate the value of this approach: (1) sampling the distance from ubiquitin to a protein of interest within the supramol. cullin-RING ligase complex, (2) stabilizing the wild-type conformation of the oncogenic mutant JAK2-V617F pseudokinase domain, and (3) inducing an activated state of the stimulator of interferon genes (STING) protein obsd. upon ligand binding. For examples 2 and 3, we present statistical anal. of meta-eABF free energy ests. and, for each case, code for reproducing this work.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3sXhtFKjtLzJ&md5=c4e43cfe37b8575041cd6e06467f2c05
    21. 21
      Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604– 613,  DOI: 10.1016/j.cpc.2013.09.018
      21
      PLUMED 2: New feathers for an old bird
      Tribello, Gareth A.; Bonomi, Massimiliano; Branduardi, Davide; Camilloni, Carlo; Bussi, Giovanni
      Computer Physics Communications (2014), 185 (2), 604-613CODEN: CPHCBZ; ISSN:0010-4655. (Elsevier B.V.)
      Enhancing sampling and analyzing simulations are central issues in mol. simulation. Recently, we introduced PLUMED, an open-source plug-in that provides some of the most popular mol. dynamics (MD) codes with implementations of a variety of different enhanced sampling algorithms and collective variables (CVs). The rapid changes in this field, in particular new directions in enhanced sampling and dimensionality redn. together with new hardware, require a code that is more flexible and more efficient. We therefore present PLUMED 2 here-a complete rewrite of the code in an object-oriented programming language (C++). This new version introduces greater flexibility and greater modularity, which both extends its core capabilities and makes it far easier to add new methods and CVs. It also has a simpler interface with the MD engines and provides a single software library contg. both tools and core facilities. Ultimately, the new code better serves the ever-growing community of users and contributors in coping with the new challenges arising in the field.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC3sXhs1yqs7fJ&md5=292009aab558d0ef1108bb9a5f036c40
    22. 22
      Frenkel, D.; Smit, B. Introduction. In Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, 2002; Vol. 1.
      There is no corresponding record for this reference.
    23. 23
      Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids, 2nd ed.; Oxford University Press, 2017.
      There is no corresponding record for this reference.
    24. 24
      Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577– 8593,  DOI: 10.1063/1.470117
      24
      A smooth particle mesh Ewald method
      Essmann, Ulrich; Perera, Lalith; Berkowitz, Max L.; Darden, Tom; Lee, Hsing; Pedersen, Lee G.
      Journal of Chemical Physics (1995), 103 (19), 8577-93CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
      The previously developed particle mesh Ewald method is reformulated in terms of efficient B-spline interpolation of the structure factors. This reformulation allows a natural extension of the method to potentials of the form 1/rp with p ≥ 1. Furthermore, efficient calcn. of the virial tensor follows. Use of B-splines in the place of Lagrange interpolation leads to analytic gradients as well as a significant improvement in the accuracy. The authors demonstrate that arbitrary accuracy can be achieved, independent of system size N, at a cost that scales as N log(N). For biomol. systems with many thousands of atoms and this method permits the use of Ewald summation at a computational cost comparable to that of a simple truncation method of 10 Å or less.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK2MXptlehtrw%253D&md5=092a679dd3bee08da28df41e302383a7
    25. 25
      Sega, M.; Dellago, C. Long-range dispersion effects on the water/vapor interface simulated using the most common models. J. Phys. Chem. B 2017, 121, 3798– 3803,  DOI: 10.1021/acs.jpcb.6b12437
      25
      Long-Range Dispersion Effects on the Water/Vapor Interface Simulated Using the Most Common Models
      Sega, Marcello; Dellago, Christoph
      Journal of Physical Chemistry B (2017), 121 (15), 3798-3803CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
      The long-range contribution to dispersion forces is known to have a major impact on the properties of inhomogeneous fluids, and its correct treatment is increasingly recognized as being a necessary requirement to avoid cutoff-related artifacts. Although anal. corrections for quantities like the surface tension are known, these cannot take into account the structural changes induced by the long-range contributions. The authors analyze the interfacial properties of seven popular water models, comparing the results with the cutoff version of the dispersion potential. The differences in surface tension ests. are in all cases found to be less than 2 mN/m.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2sXjtVait7o%253D&md5=a32c3a9e14c07687bb492635c4b5c74d
    26. 26
      Panigrahy, R. An Improved Algorithm Finding Nearest Neighbor Using Kd-trees. In Lecture Notes in Computer Science; Springer: Berlin Heidelberg pp 387– 398.
      There is no corresponding record for this reference.
    27. 27
      Páll, S.; Zhmurov, A.; Bauer, P.; Abraham, M.; Lundborg, M.; Gray, A.; Hess, B.; Lindahl, E. Heterogeneous parallelization and acceleration of molecular dynamics simulations in GROMACS. J. Chem. Phys. 2020, 153, 134110  DOI: 10.1063/5.0018516
      27
      Heterogeneous parallelization and acceleration of molecular dynamics simulations in GROMACS
      Pall, Szilard; Zhmurov, Artem; Bauer, Paul; Abraham, Mark; Lundborg, Magnus; Gray, Alan; Hess, Berk; Lindahl, Erik
      Journal of Chemical Physics (2020), 153 (13), 134110CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
      The introduction of accelerator devices such as graphics processing units (GPUs) has had profound impact on mol. dynamics simulations and has enabled order-of-magnitude performance advances using commodity hardware. To fully reap these benefits, it has been necessary to reformulate some of the most fundamental algorithms, including the Verlet list, pair searching, and cutoffs. Here, we present the heterogeneous parallelization and acceleration design of mol. dynamics implemented in the GROMACS codebase over the last decade. The setup involves a general cluster-based approach to pair lists and non-bonded pair interactions that utilizes both GPU and central processing unit (CPU) single instruction, multiple data acceleration efficiently, including the ability to load-balance tasks between CPUs and GPUs. The algorithm work efficiency is tuned for each type of hardware, and to use accelerators more efficiently, we introduce dual pair lists with rolling pruning updates. Combined with new direct GPU-GPU communication and GPU integration, this enables excellent performance from single GPU simulations through strong scaling across multiple GPUs and efficient multi-node parallelization. (c) 2020 American Institute of Physics.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3cXitVWrtLnE&md5=6997d6a941338d496269c295828c22f9
    28. 28
      Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; in ’t Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D.; Shan, R.; Stevens, M. J.; Tranchida, J.; Trott, C.; Plimpton, S. J. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comput. Phys. Commun. 2022, 271, 108171  DOI: 10.1016/j.cpc.2021.108171
      28
      LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
      Thompson, Aidan P.; Aktulga, H. Metin; Berger, Richard; Bolintineanu, Dan S.; Brown, W. Michael; Crozier, Paul S.; in 't Veld, Pieter J.; Kohlmeyer, Axel; Moore, Stan G.; Nguyen, Trung Dac; Shan, Ray; Stevens, Mark J.; Tranchida, Julien; Trott, Christian; Plimpton, Steven J.
      Computer Physics Communications (2022), 271 (), 108171CODEN: CPHCBZ; ISSN:0010-4655. (Elsevier B.V.)
      Since the classical mol. dynamics simulator LAMMPS was released as an open source code in 2004, it has become a widely-used tool for particle-based modeling of materials at length scales ranging from at. to mesoscale to continuum. Reasons for its popularity are that it provides a wide variety of particle interaction models for different materials, that it runs on any platform from a single CPU core to the largest supercomputers with accelerators, and that it gives users control over simulation details, either via the input script or by adding code for new interat. potentials, constraints, diagnostics, or other features needed for their models. As a result, hundreds of people have contributed new capabilities to LAMMPS and it has grown from fifty thousand lines of code in 2004 to a million lines today. In this paper several of the fundamental algorithms used in LAMMPS are described along with the design strategies which have made it flexible for both users and developers. We also highlight some capabilities recently added to the code which were enabled by this flexibility, including dynamic load balancing, on-the-fly visualization, magnetic spin dynamics models, and quantum-accuracy machine learning interat. potentials.Program Title: Large-scale Atomic/Mol. Massively Parallel Simulator (LAMMPS)CPC Library link to program files:https://doi.org/10.17632/cxbxs9btsv.1Developer's repository link:https://github.com/lammps/lammpsLicensing provisions: GPLv2Programming language: C++, Python, C, FortranSupplementary material:https://www.lammps.orgNature of problem: Many science applications in physics, chem., materials science, and related fields require parallel, scalable, and efficient generation of long, stable classical particle dynamics trajectories. Within this common problem definition, there lies a great diversity of use cases, distinguished by different particle interaction models, external constraints, as well as timescales and lengthscales ranging from at. to mesoscale to macroscopic.Soln. method: The LAMMPS code uses parallel spatial decompn., distributed neighbor lists, and parallel FFTs for long-range Coulombic interactions [1]. The time integration algorithm is based on the Stormer-Verlet symplectic integrator [2], which provides better stability than higher-order non-symplectic methods. In addn., LAMMPS supports a wide range of interat. potentials, constraints, diagnostics, software interfaces, and pre- and post-processing features.Addnl. comments including restrictions and unusual features: This paper serves as the definitive ref. for the LAMMPS code.S. Plimpton, Fast parallel algorithms for short-range mol. dynamics. Phys. 117 (1995) 1-19.L. Verlet, Computer expts. on classical fluids: I. Thermodynamical properties of Lennard-Jones mols., Phys. Rev. 159 (1967) 98-103.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3MXitlSrsb7O&md5=cd0bfd050820e97c11779003add20ed3
    29. 29
      Helfrich, W. Elastic Properties of Lipid Bilayers: Theory and Possible Experiments. Z. Naturforsch. C 1973, 28, 693– 703,  DOI: 10.1515/znc-1973-11-1209
      29
      Elastic properties of lipid bilayers. Theory and possible experiments
      Helfrich, W.
      Zeitschrift fuer Naturforschung, Teil C: Biochemie, Biophysik, Biologie, Virologie (1973), 28 (11-12), 693-703CODEN: ZNFCAP; ISSN:0341-0471.
      A theory for the elastic behavior of lipid bilayers is presented involving 3 types of strains, i.e. stretching, tilt, and curvature. Possible expts. to det. the elastic properties of lipid bilayers by magnetic and pressure deformation measurements are discussed.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaE2cXnvFersA%253D%253D&md5=c64c8c2c43d3defe6b78d2c06abe0630
    30. 30
      Bhaskara, R. M.; Grumati, P.; Garcia-Pardo, J.; Kalayil, S.; Covarrubias-Pinto, A.; Chen, W.; Kudryashev, M.; Dikic, I.; Hummer, G. Curvature induction and membrane remodeling by FAM134B reticulon homology domain assist selective ER-phagy. Nat. Commun. 2019, 10, 2370  DOI: 10.1038/s41467-019-10345-3
      30
      Curvature induction and membrane remodeling by FAM134B reticulon homology domain assist selective ER-phagy
      Bhaskara Ramachandra M; Hummer Gerhard; Grumati Paolo; Covarrubias-Pinto Adriana; Dikic Ivan; Garcia-Pardo Javier; Kalayil Sissy; Chen Wenbo; Kudryashev Mikhail; Dikic Ivan; Garcia-Pardo Javier; Chen Wenbo; Kudryashev Mikhail; Hummer Gerhard
      Nature communications (2019), 10 (1), 2370 ISSN:.
      FAM134B/RETREG1 is a selective ER-phagy receptor that regulates the size and shape of the endoplasmic reticulum. The structure of its reticulon-homology domain (RHD), an element shared with other ER-shaping proteins, and the mechanism of membrane shaping remain poorly understood. Using molecular modeling and molecular dynamics (MD) simulations, we assemble a structural model for the RHD of FAM134B. Through MD simulations of FAM134B in flat and curved membranes, we relate the dynamic RHD structure with its two wedge-shaped transmembrane helical hairpins and two amphipathic helices to FAM134B functions in membrane-curvature induction and curvature-mediated protein sorting. FAM134B clustering, as expected to occur in autophagic puncta, amplifies the membrane-shaping effects. Electron microscopy of in vitro liposome remodeling experiments support the membrane remodeling functions of the different RHD structural elements. Disruption of the RHD structure affects selective autophagy flux and leads to disease states.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BB3M3hs1ahug%253D%253D&md5=10b720cc239064d028e0436dc9e1ae4d
    31. 31
      Wassenaar, T. A.; Ingólfsson, H. I.; Böckmann, R. A.; Tieleman, D. P.; Marrink, S. J. Computational Lipidomics with insane: A Versatile Tool for Generating Custom Membranes for Molecular Simulations. J. Chem. Theory Comput. 2015, 11, 2144– 2155,  DOI: 10.1021/acs.jctc.5b00209
      31
      Computational Lipidomics with insane: A Versatile Tool for Generating Custom Membranes for Molecular Simulations
      Wassenaar, Tsjerk A.; Ingolfsson, Helgi I.; Boeckmann, Rainer A.; Tieleman, D. Peter; Marrink, Siewert J.
      Journal of Chemical Theory and Computation (2015), 11 (5), 2144-2155CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
      For simulations of membranes and membrane proteins, the generation of the lipid bilayer is a crit. step in the setup of the system. Membranes comprising multiple components pose a particular challenge, because the relative abundances need to be controlled and the equilibration of the system may take several microseconds. Here the authors present a comprehensive method for building membrane contg. systems, characterized by simplicity and versatility. The program uses preset, coarse-grain lipid templates to build the membrane, and also allows on-the-fly generation of simple lipid types by specifying the headgroup, linker, and lipid tails on the command line. The resulting models can be equilibrated, after which a relaxed atomistic model can be obtained by reverse transformation. For multicomponent membranes, this provides an efficient means for generating equilibrated atomistic models. The method is called insane, an acronym for INSert membrANE. The program has been made available, together with the complementary method for reverse transformation, at http://cgmartini.nl/. This work highlights the key features of insane and presents a survey of properties for a large range of lipids as a start of a computational lipidomics project.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2MXmtlamsbo%253D&md5=8a754c43db80752045ed3bfc3e64b5db
    32. 32
      Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101  DOI: 10.1063/1.2408420
      32
      Canonical sampling through velocity rescaling
      Bussi, Giovanni; Donadio, Davide; Parrinello, Michele
      Journal of Chemical Physics (2007), 126 (1), 014101/1-014101/7CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
      The authors present a new mol. dynamics algorithm for sampling the canonical distribution. In this approach the velocities of all the particles are rescaled by a properly chosen random factor. The algorithm is formally justified and it is shown that, in spite of its stochastic nature, a quantity can still be defined that remains const. during the evolution. In numerical applications this quantity can be used to measure the accuracy of the sampling. The authors illustrate the properties of this new method on Lennard-Jones and TIP4P water models in the solid and liq. phases. Its performance is excellent and largely independent of the thermostat parameter also with regard to the dynamic properties.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD2sXosVCltg%253D%253D&md5=9c182b57bfc65bca6be23c8c76b4be77
    33. 33
      Bernetti, M.; Bussi, G. Pressure control using stochastic cell rescaling. J. Chem. Phys. 2020, 153, 114107  DOI: 10.1063/5.0020514
      33
      Pressure control using stochastic cell rescaling
      Bernetti, Mattia; Bussi, Giovanni
      Journal of Chemical Physics (2020), 153 (11), 114107CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
      Mol. dynamics simulations require barostats to be performed at a const. pressure. The usual recipe is to employ the Berendsen barostat first, which displays a first-order vol. relaxation efficient in equilibration but results in incorrect vol. fluctuations, followed by a second-order or a Monte Carlo barostat for prodn. runs. In this paper, we introduce stochastic cell rescaling, a first-order barostat that samples the correct vol. fluctuations by including a suitable noise term. The algorithm is shown to report vol. fluctuations compatible with the isobaric ensemble and its anisotropic variant is tested on a membrane simulation. Stochastic cell rescaling can be straightforwardly implemented in the existing codes and can be used effectively in both equilibration and prodn. phases. (c) 2020 American Institute of Physics.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3cXhvVOnt7nN&md5=0021e74a81e9956f6b64873d904a66f2
    34. 34
      Mark, P.; Nilsson, L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A 2001, 105, 9954– 9960,  DOI: 10.1021/jp003020w
      34
      Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K
      Mark, Pekka; Nilsson, Lennart
      Journal of Physical Chemistry A (2001), 105 (43), 9954-9960CODEN: JPCAFH; ISSN:1089-5639. (American Chemical Society)
      Mol. dynamics simulations of five water models, the TIP3P (original and modified), SPC (original and refined), and SPC/E (original), were performed using the CHARMM mol. mechanics program. All simulations were carried out in the microcanonical NVE ensemble, using 901 water mols. in a cubic simulation cell furnished with periodic boundary conditions at 298 K. The SHAKE algorithm was used to keep water mols. rigid. Nanosecond trajectories were calcd. with all water models for high statistical accuracy. The characteristic self-diffusion coeffs. D and radial distribution functions, gOO, gOH, and gHH for all five water models were detd. and compared to exptl. data. The effects of velocity rescaling on the self-diffusion coeff. D were examd. All these empirical water models used in this study are similar by having three interaction sites, but the small differences in their pair potentials composed of Lennard-Jones (LJ) and Coulombic terms give significant differences in the calcd. self-diffusion coeffs., and in the height of the second peak of the radial distribution function gOO.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD3MXntlWrurs%253D&md5=fecd3db40210b04e8b2a933ea07b131e
    35. 35
      Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: A web-based graphical user interface for CHARMM. J. Comput. Chem. 2008, 29, 1859– 1865,  DOI: 10.1002/jcc.20945
      35
      CHARMM-GUI: a web-based graphical user interface for CHARMM
      Jo, Sunhwan; Kim, Taehoon; Iyer, Vidyashankara G.; Im, Wonpil
      Journal of Computational Chemistry (2008), 29 (11), 1859-1865CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
      CHARMM is an academic research program used widely for macromol. mechanics and dynamics with versatile anal. and manipulation tools of at. coordinates and dynamics trajectories. CHARMM-GUI, http://www.charmm-gui.org, has been developed to provide a web-based graphical user interface to generate various input files and mol. systems to facilitate and standardize the usage of common and advanced simulation techniques in CHARMM. The web environment provides an ideal platform to build and validate a mol. model system in an interactive fashion such that, if a problem is found through visual inspection, one can go back to the previous setup and regenerate the whole system again. In this article, we describe the currently available functional modules of CHARMM-GUI Input Generator that form a basis for the advanced simulation techniques. Future directions of the CHARMM-GUI development project are also discussed briefly together with other features in the CHARMM-GUI website, such as Archive and Movie Gallery.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD1cXosVKksbc%253D&md5=112a3dd61d792b040f9f716b32220d7e
    36. 36
      Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985, 31, 1695– 1697,  DOI: 10.1103/PhysRevA.31.1695
      36
      Canonical dynamics: Equilibrium phase-space distributions
      Hoover
      Physical review. A, General physics (1985), 31 (3), 1695-1697 ISSN:0556-2791.
      There is no expanded citation for this reference.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BC2sjotlWltA%253D%253D&md5=99a2477835b37592226a5d18a760685c
    37. 37
      Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511– 519,  DOI: 10.1063/1.447334
      37
      A unified formulation of the constant-temperature molecular-dynamics methods
      Nose, Shuichi
      Journal of Chemical Physics (1984), 81 (1), 511-19CODEN: JCPSA6; ISSN:0021-9606.
      Three recently proposed const. temp. mol. dynamics methods [N., (1984) (1); W. G. Hoover et al., (1982) (2); D. J. Evans and G. P. Morris, (1983) (2); and J. M. Haile and S. Gupta, 1983) (3)] are examd. anal. via calcg. the equil. distribution functions and comparing them with that of the canonical ensemble. Except for effects due to momentum and angular momentum conservation, method (1) yields the rigorous canonical distribution in both momentum and coordinate space. Method (2) can be made rigorous in coordinate space, and can be derived from method (1) by imposing a specific constraint. Method (3) is not rigorous and gives a deviation of order N-1/2 from the canonical distribution (N the no. of particles). The results for the const. temp.-const. pressure ensemble are similar to the canonical ensemble case.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaL2cXkvFOrs7k%253D&md5=2974515ec89e5601868e35871c0f19c2
    38. 38
      Welch, P. The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Trans. Audio Electroacoust. 1967, 15, 70– 73,  DOI: 10.1109/TAU.1967.1161901
      There is no corresponding record for this reference.
    39. 39
      Hernández-Muñoz, J.; Bresme, F.; Tarazona, P.; Chacón, E. Bending Modulus of Lipid Membranes from Density Correlation Functions. J. Chem. Theory Comput. 2022, 18, 3151– 3163,  DOI: 10.1021/acs.jctc.2c00099
      39
      Bending Modulus of Lipid Membranes from Density Correlation Functions
      Hernandez-Munoz, Jose; Bresme, Fernando; Tarazona, Pedro; Chacon, Enrique
      Journal of Chemical Theory and Computation (2022), 18 (5), 3151-3163CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
      The bending modulus κ quantifies the elasticity of biol. membranes in terms of the free energy cost of increasing the membrane corrugation. Mol. dynamics (MD) simulations provide a powerful approach to quantify κ by analyzing the thermal fluctuations of the lipid bilayer. However, existing methods require the identification and filtering of non-mesoscopic fluctuation modes. State of the art methods rely on identifying a smooth surface to describe the membrane shape. These methods introduce uncertainties in calcg. κ since they rely on different criteria to select the relevant fluctuation modes. Here, we present a method to compute κ using mol. simulations. Our approach circumvents the need to define a mesoscopic surface or an orientation field for the lipid tails explicitly. The bending and tilt moduli can be extd. from the anal. of the d. correlation function (DCF). The method introduced here builds on the Bedeaux and Weeks (BW) theory for the DCF of fluctuating interfaces and on the coupled undulatory (CU) mode introduced by us in previous work. We test the BW-DCF method by computing the elastic properties of lipid membranes with different system sizes (from 500 to 6000 lipid mols.) and using coarse-grained (for POPC and DPPC lipids) and fully atomistic models (for DPPC). Further, we quantify the impact of cholesterol on the bending modulus of DPPC bilayers. We compare our results with bending moduli obtained with X-ray diffraction data and different computer simulation methods.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB38XptVOitLs%253D&md5=1ce4b80044739c9c7b78430e2d1aa421
    40. 40
      Eid, J.; Razmazma, H.; Jraij, A.; Ebrahimi, A.; Monticelli, L. On Calculating the Bending Modulus of Lipid Bilayer Membranes from Buckling Simulations. J. Phys. Chem. B 2020, 124, 6299– 6311,  DOI: 10.1021/acs.jpcb.0c04253
      40
      On calculating the bending modulus of lipid bilayer membranes from buckling simulations
      Eid, Jad; Razmazma, Hafez; Jraij, Alia; Ebrahimi, Ali; Monticelli, Luca
      Journal of Physical Chemistry B (2020), 124 (29), 6299-6311CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
      The bending modulus is an important phys. const. characterizing lipid membranes. Different methods have been devised for calcg. the bending modulus from simulations, and one of them, named the buckling method, is nowadays widely used due to its simplicity and numerical stability. However, questions remain on the reproducibility, finite size effects, and interpretation of results on lipid mixts. Here we explore the dependence of simulation results on the system size and the strain. We find that the dimensions of the box have a negligible impact on the results when the system size is beyond a certain threshold. We then calc. the bending rigidity for of a series of common single-component lipid bilayers (PC, PS, PE, PG, and SM), as well as a no. of binary and ternary lipid mixts. We find that bending moduli of lipid mixts. can be predicted from the weighted av. of the moduli of the individual components, as long as the mixt. is homogeneous. For phase-sepd. mixts., the apparent elastic modulus is closer to the value of the softer component. Predictions of the bending modulus based on the area compressibility modulus are found to be generally unreliable.
      >> More from SciFinder ®
      https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3cXht1yntbnL&md5=3aa840e0da955031d874d2b047514f6f

  • Supporting Information

    Supporting Information

    ARTICLE SECTIONS
    Jump To

    The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.3c00777.

    • Detailed derivation of number of missed interactions; the rigid-rotor model for missed interactions in systems with rigid or near-rigid molecules; power spectral density of pressure fluctuations in TIP3P water; the shape of the average pressure tensor between neighbor list updates depends on the update frequencies of the inner and outer list; difference ΔP in the pressure just before and right after a neighbor list update in the NVT TIP3P solvent system; in a simulation of the large Martini membrane system with nstpcouple = nstlist = 25, the unphysical box distortion is greatly suppressed; power spectral density (PSD) of the pressure fluctuations in MD simulations of TIP3P water in the NVT ensemble with nstlist = 20; probability pmissed of missed interactions (gray line) as function of the VBT parameter in GROMACS for the NVT Martini water system (PDF)


    Terms & Conditions

    Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

PDF [ 10MB]

Pair your accounts.

Export articles to Mendeley

Get article recommendations from ACS based on references in your Mendeley library.

Pair your accounts.

Export articles to Mendeley

Get article recommendations from ACS based on references in your Mendeley library.

You’ve supercharged your research process with ACS and Mendeley!

STEP 1:
Click to create an ACS ID

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

MENDELEY PAIRING EXPIRED
Your Mendeley pairing has expired. Please reconnect