# Small ionic radii limit time step in Martini 3 molecular dynamics
simulations

Balázs Fábián<EMAIL_ADDRESS>[  Sebastian Thallmair  [  Gerhard
Hummer  [

###### Abstract

Among other improvements, the Martini 3 coarse-grained force field provides a
more accurate description of the solvation of protein pockets and channels
through the consistent use of various bead types and sizes. Here, we show that
the representation of Na+ and Cl- ions as “tiny” (TQ5) beads limits the
accessible time step to 25 fs. By contrast, with Martini 2, time steps of
30-40 fs were possible for lipid bilayer systems without proteins. This
limitation is relevant for, e.g., phase separating lipid mixtures that require
long equilibration times. We derive a quantitative kinetic model of time-
integration instabilities in molecular dynamics (MD) as a function of time
step, ion concentration and mass, system size, and simulation time. With this
model, we demonstrate that ion-water interactions are the main source of
instability at physiological conditions, followed closely by ion-ion
interactions. We show that increasing the ionic masses makes it possible to
use time steps up to 40 fs with minimal impact on static equilibrium
properties and on dynamical quantities such as lipid and ion diffusion
coefficients. Increasing the size of the bead representing the ions (and thus
changing their hydration) also permits longer time steps. The use of larger
time steps in Martini 3 simulations results in a more efficient exploration of
configuration space. The kinetic model of MD simulation crashes can be used to
determine the maximum allowed time step whenever sampling efficiency is
critical.

MPI-BP]Department of Theoretical Biophysics, Max Planck Institute of
Biophysics, Max-von-Laue Straße 3, 60438, Frankfurt am Main, Germany
FIAS]Frankfurt Institute for Advanced Studies, Ruth-Moufang-Straße 1, 60438
Frankfurt am Main, Germany MPI-BP]Department of Theoretical Biophysics, Max
Planck Institute of Biophysics, Max-von-Laue Straße 3, 60438, Frankfurt am
Main, Germany

##  1 Introduction

Compared to its predecessors1, 2, the recent Martini 3 force field3
constitutes a significant advance in biomolecular simulations. Martini 3
offers consistently parametrized coarse-grained (CG) interaction sites of
different size and a better coverage of the chemical space by a whole range of
new bead types. This results in an improved representation of molecular shape,
packing, and interactions in general. The diversity of the ion models was
strongly increased, now featuring five bead types Q1-Q5 in three resolutions
ranging from “tiny” (T, 2-to-1 mapping) via “small” (S, 3-to-1 mapping) to
“regular” (R, 4-to-1 mapping) beads. The parametrization of the beads features
a Martini Hofmeister series from hard, more inorganic ions (Q5) to soft, more
organic ions (Q1). In addition, specific interactions (e.g., of cation-$\pi$
type) can be incorporated using labels which were previously only available in
the polarizable version. 4 The new features of Martini 3 resulted in a change
of the resolution of hard ions such as Na+ and Cl- which are routinely added
to biomolecular simulation systems to establish physiological conditions in
terms of ionic strength. The size of both Na+ and Cl- ions changed from R to
T, which means that in Martini 3 the ions are modeled without a hydration
shell. This is in contrast to Martini 2, where a hydration shell was
considered to be included in the larger R-bead type used for Na+ and Cl-.

While Martini 2 simulations were routinely run with time steps of $\Delta
t=30$ fs 5, 6, 7 according to the recommendation of the new-rf parameter set
8, the current recommended time step is $\Delta t=20$ fs, which was used for
all test systems during the parametrization of the force field. In several
previous studies employing Martini 2 even $\Delta t=40$ fs was used 2, 9, 10.
However, this was deemed excessive by the developers of the model 11, and the
use of such large time steps is generally discouraged.

Here, we show that the representation of Martini 3 ions as “tiny” charged
beads of type TQ5 requires integration time step no larger than 25 fs. For
$\Delta t>25$ fs the tendency of the time integration to become unstable and
crash the simulations increases significantly. We develop a quantitative
kinetic model of the rate of crashing as a function of integration time step,
ion concentration, and ionic mass. By fitting the three parameters of this
model globally to the observed crash statistics of the MD runs, we identify
ion-ion and ion-water interactions as the main culprits. We explore two
possible solutions that permit simulations to run with larger time steps: by
modifying the ionic mass and by changing the bead type of the ions. While the
former keeps the desired resolution with almost no perturbations to the
system, the latter is consistent with the multi-resolution approach of Martini
3. The effects and trade-offs made by these solutions are evaluated and
discussed in detail. Moreover, the statistical model employed here to estimate
the rate of crashing can be applied to determine an optimal time step for MD-
based high-throughput campaigns. The formalism used here to determine the
maximum allowed time step should prove useful whenever sampling efficiency is
critical.

##  2 Theory

###  2.1 Altering the mass or bead type

In classical molecular dynamics (MD) simulations, Newton’s equations of motion
(usually augmented with thermostats and barostats) are integrated numerically
with finite time steps. The integration time step $\Delta t$ is chosen in a
trade off between accuracy, e.g., to conserve the total energy, and the
computational efficiency of sampling configuration space. To maximize
efficiency, the time step is commonly chosen close to the stability limit of
time integration. 12 For classical non-polarizable models, this limit is
determined by the fastest molecular motions, in particular the rattling of
stiff covalent bonds and tight non-covalent interations, and the hard
collisions between molecules on highly anharmonic potential surfaces.

For NaCl solutions and lipid bilayers, we observed that Martini 3 MD
simulations required shorter time steps than in equivalent setups simulated
with the Martini 2 model. As the main culprit, we identified the newly
introduced ions modeled as “tiny” charged (TQ5) beads and their ion-ion
interactions. As possible remedies, we consider (i) changing the ionic masses
and (ii) changing the bead type from tiny to small (SQ5) or regular (RQ5).

The motivation behind changing the ionic mass is the invariance of Newton’s
equations of motion under a uniform scaling of all masses and of time by
$\alpha$ and $\sqrt{\alpha}$, respectively. Whereas scaling all particle
masses does not improve configuration space sampling, we expect that
increasing the ionic mass will permit a longer time step. The displacements of
heavier ions at each time step are smaller, which should reduce the risk of an
uncontrolled clash causing an instability in the numerical time integration.
If such ion displacements are limiting the time step, we expect that the
maximum allowed time step for stable time integration scales with the ionic
mass $m$ as

$\displaystyle\frac{\Delta t^{\prime}}{\Delta
t}\approx\sqrt{\frac{m^{\prime}}{m}}$ (1)

where we ignored the motions of the other particles. Then, starting from a
time step $\Delta t=20$ fs, doubling and quadrupling the ionic mass $m$ would
make it possible to use $\Delta t^{\prime}=30$ and 40 fs, respectively.

Importantly, changes in particle mass leave the equilibrium structural
properties unchanged because the classical mechanical partition function
separates into kinetic and configurational contributions. However, the
dynamics of the system is generally modified in a nontrivial manner.
Nevertheless, with a longer time step one can expect a faster sampling of the
relevant configuration space for the molecules with unmodified masses.

In Martini simulations, usually only relatively slow diffusive motions are of
interest, because the fast intramolecular and collision dynamics is impacted
significantly by coarse-graining. Translational and rotational diffusion are
dominated by hydrodynamic effects, as reflected in the Stokes–Einstein theory
of 3D diffusion 13 and in the Saffman–Delbrück theory of 2D diffusion in lipid
bilayers embedded in a 3D solvent 14. Both theories predict that the diffusion
coefficient of a dissolved molecular species does not depend directly on its
mass, only on its size and the viscosity of the surrounding media. Therefore,
we expect that ionic mass changes have minimal impact on the diffusive
dynamics.

Nevertheless, changes in the mass of some of the particles can potentially
change the viscosity of the aqueous solvent and the membrane. Here, we
consider increasing the mass of the ions described by TQ5 beads, leaving all
other masses untouched. Importantly, these particles constitute only a tiny
fraction of the molecules in the aqueous solvent and do not enter within the
lipid bilayers. Therefore, we expect that doubling their mass will have only a
very limited effect on the viscosity of solvent and membrane, and thus in turn
on the diffusive dynamics of all species, including the ions themselves. We
test this assumption by monitoring the diffusivity of solvent and lipid
species as a function of ionic mass. By keeping the box dimensions constant in
this comparison, we do not need to correct for the large finite size effects
in the computed diffusion coefficients 15, 16.

Modification of the bead type of the ions is based on the notion that in
Martini 3, a change in the bead size of the ions represents a change in their
hydration. Bead-size changes are thus consistent with the coarse-graining
philosophy of the force field (see Martini 3 Supplementary Information 3).
Moreover, the current resolution of lipid head groups such as the choline
moiety of phosphatidyl-choline (PC) lipids is already low due to the employed
6-1 mapping 2, 17, where 6 non-hydrogen atoms are grouped together into a
single interaction site. Therefore, adequate modeling of protein-less bilayers
currently does not require a high resolution of the ions. We thus expect that
an increase in the ion bead size will not significantly perturb the membrane
properties.

###  2.2 Crash probability as function of time step

In the following, we develop a kinetic model of the probability that a
simulation with a given time step will crash during a given total simulation
time. We assume that, rarely, the position updates in time integration place
particles within a strongly repulsive region of the potential surface, such
that the subsequent force evaluation leads to numerical instabilities either
immediately or after a few uncontrolled further integration steps.

Let $\Delta t$ be the integration time step. With $v$ the velocity in one
dimension, the corresponding displacement is then $\Delta x=v\Delta t$,
ignoring the higher-order acceleration terms. We now assume that displacements
$\Delta x>\Delta x_{\mathrm{crit}}$ are critical, and lead to crashes. The
velocities of a particle of mass $m$ at inverse temperature $\beta$ satisfy a
Maxwell-Boltzmann distribution, $p(v)dv=\exp(-\beta mv^{2}/2)\sqrt{\beta
m/(2\pi)}dv$. The probability of a critical displacement during a simulation
time step becomes

$\displaystyle P(\Delta x>\Delta x_{\mathrm{crit}})$
$\displaystyle=\int_{\Delta x_{\mathrm{crit}}/\Delta
t}^{\infty}dv\,\exp(-\beta mv^{2}/2)\sqrt{\beta m/(2\pi)}$
$\displaystyle\approx\frac{\gamma\Delta t\exp(-1/(2\gamma^{2}\Delta
t^{2}))}{\sqrt{2\pi}}$ (2)

with $1/\gamma^{2}=\beta m\Delta x_{\mathrm{crit}}^{2}$. Here we used the
asymptotic expansion of the complementary error function. In a simulation
system, we expect that a certain fraction $\sqrt{2\pi}f$ of the particle
displacements (here, those of the “tiny” ions and their interaction partners)
can result in such catastrophic clashes. The probability
$Q=\sqrt{2\pi}fP(\Delta x>\Delta x_{\mathrm{crit}})$ of a catastrophic
displacement in a given time step is approximately

$Q=f\gamma\Delta t\exp(-1/(2\gamma^{2}\Delta t^{2}))$ (3)

Assuming uncorrelated events and $0<Q\ll 1$, the probability that a simulation
run of $n$ time steps will not crash is

$P=(1-Q)^{n}\approx\exp\left[-tf\gamma\exp(-1/(2\gamma^{2}\Delta
t^{2}))\right]$ (4)

with $t=n\Delta t$ the total time. The crash times are thus predicted to be
distributed exponentially,
$p(t_{\text{crash}})dt_{\text{crash}}=k_{\mathrm{crash}}\exp(-k_{\mathrm{crash}}t_{\text{crash}})dt_{\text{crash}}$
with a rate that depends on the simulation time step as

$k_{\mathrm{crash}}(\Delta t)=f\gamma\exp(-1/(2\gamma^{2}\Delta
t^{2}))=k_{\mathrm{crash}}^{0}\exp(-1/(2\gamma^{2}\Delta t^{2}))$ (5)

The crash rate has two parameters, $f$ and $\gamma$, whose product defines
$k_{\mathrm{crash}}^{0}=f\gamma$. In our systems, we expect that the
dimensionless factor $f$ depends on the concentration $c_{\mathrm{ion}}$ of
ions as $f=(a{c_{\mathrm{ion}}}^{2}+bc_{\mathrm{ion}})V$ with coefficients $a$
and $b$, and $V$ the system volume. Based on the definition, and consistent
with Eq. (1), we expect $\gamma$ to depend on the mass of the ions
approximately as $\gamma\propto 1/\sqrt{m}$.

##  3 Methods

###  3.1 Simulation of NaCl solutions

To demonstrate the effect of changing ionic masses (or the lack thereof) on
the viscosity $\mu_{\mathrm{f}}$ of the fluid surrounding the bilayer, we
prepared water boxes with 0, 0.075, 0.15 and 0.3 M NaCl concentrations using
insane.py17 and the option -salt. All cubic boxes had an initial edge length
of 10 nm. We simulated these NaCl solutions with the standard Martini 3 CG
force field3 and the Gromacs package18 (version 2020.1). To eliminate
differences stemming from finite size effects between these systems15, we took
the smallest box as reference and removed the appropriate number of molecules
from all others. This procedure resulted in equilibrium box lengths between
8.93 and 8.97 nm for all NaCl solution systems. No corrections were applied to
account for finite size effects on diffusion. At every concentration, we
performed simulations with time steps of 20, 30 and 40 fs, and using the
original, doubled and quadrupled ionic masses (denoted as 1$\times$m,
2$\times$m, 4$\times$m, respectively) to systematically test the effects on
the diffusion coefficient, and therefore on the viscosity of the bulk medium.
Following a short energy minimization using a steepest descent algorithm,
every system was equilibrated during a 1 µs run with a 20 fs time step. This
was followed by 3 µs-long production runs using 20, 30 or 40 fs time steps
with writing the positions to disk every 0.5 ns. The temperature of the
systems was kept at 310 K using a velocity re-scale thermostat 19. The
pressure was maintained at 1 bar by an isotropic Parrinello–Rahman barostat 20
(p${}_{\mathrm{ref}}=1$ bar and $\tau_{\mathrm{p}}=12$ ps).

To gain insight into the origin of the simulation crashes, we performed a
large number of additional simulations using various combinations of
concentration (0.075 M and in the range of 0.1 to 0.53 M with steps of
approximately 0.05 M) and time step (every fs from 29 to 32 fs). In this set,
every simulation used the settings as described above, but with a fixed number
of desired steps ($n_{\mathrm{steps}}=10^{8}$) instead of fixed total
simulation time, and without trying to match the overall volumes. The majority
of these simulations could not achieve the required
$n_{\mathrm{steps}}=10^{8}$, and the only observable extracted from them was
the rate of crashing, as described below. To reduce the statistical
uncertainty in the evaluation of crash rate, every simulation was repeated 40
times.

Figure 1: Snapshots of the NaCl solution system (LEFT) and the DPPC lipid
bilayer system (RIGHT). Sodium and chloride ions are shown as yellow and green
spheres, and DPPC lipids as orange tubes. Water is represented as transparent
surface. The images were rendered using VMD 21.

###  3.2 Simulation of DPPC bilayers

We also performed MD simulations of fully solvated DPPC bilayers containing
Na+ and Cl- ions with the Martini 2 and Martini 3 CG force fields. The initial
systems were prepared using insane.py17. All bilayers were constructed in an
initial box of 18$\times$18$\times$18 nm3 and solvated by 14405 CG waters
(pure W for Martini 3, 10% WF antifreeze2 for Martini 2). The equilibrium
membrane area in the simulation boxes was about 17x17 nm2, larger than the
$\approx$8x8 nm2 used in the Martini 3 validation test. Every simulation
system contained 0.15 M NaCl, as added with insane.py using the flag -salt
0.15. We performed all simulations with Gromacs 2020.1, and used 20, 30 and 40
fs time steps. We simulated a Martini 2 membrane for reference, and Martini 3
membranes using TQ5 (standard Martini 3), SQ5, RQ5(=Q5) beads as monovalent
ions. Finally, we tested the effect of increasing the ionic mass by a factor
of 2. The doubled mass corresponds to using the interactions of TQ5 along with
the mass of RQ5. The abbreviated labels of the above lipid simulations are M2,
M3(-TQ5), M3-SQ5, M3-RQ5 and M3-2$\times$m, respectively. All membranes were
energy minimized using a steepest descent algorithm, and further minimized by
running 50 ns trajectories using 5 fs time steps. This was followed by a 5 µs
equilibration using a 20 fs time step. Then, production runs of 10 µs were
performed using 20, 30 and 40 fs time steps and writing the positions to disk
every 0.5 ns. The temperature of the systems was kept at 310 K by a velocity
re-scale thermostat19. The pressure of the systems was controlled in a semi-
isotropic manner with a compressibility of $3\times 10^{-4}$ bar-1. For the
initial energy minimization, we employed the Berendsen barostat22
(p${}_{\mathrm{ref}}=1$ bar and $\tau_{\mathrm{p}}=1$ ps), while the
equilibration and production proceeded with the Parrinello-Rahman barostat20
($p_{\mathrm{ref}}=1$ bar and $\tau_{\mathrm{p}}=12$ ps). Figure 1 shows
illustrative snapshots of the simulated NaCl solutions and lipid bilayers.

###  3.3 Analysis

In our numerical tests, we ran 40 simulations for each system. The time of
$t_{i}$ is the targeted run length, $t_{i}=\Delta t\times 10^{8}$, if run $i$
completed normally, or the time of the crash caused by an instability. For the
exact definition of a crash, see section 1.1 of SI. From the times $t_{i}$, we
estimate the rate of crashing (and thus the reciprocal of the mean time to a
crash) using a maximum-likelihood estimator for right-censored data and
Poisson statistics,

$k_{\mathrm{crash}}=\frac{1}{t_{\mathrm{crash}}}=\frac{n_{\mathrm{crash}}}{\sum_{i=1}^{40}t_{i}},$
(6)

where $0\leq n_{\mathrm{crash}}\leq 40$ is the number of crashes observed in
the 40 runs. We assessed the uncertainty of the estimator using the standard
error $\sigma=k_{\mathrm{crash}}/\sqrt{n_{\mathrm{crash}}}$ which corresponds
to the Cramér–Rao bound. This estimator gives us some guidance on the
molecular factors causing the crashes. For instance, if ion-water interactions
were dominant, we would expect $k_{\mathrm{crash}}\propto c_{\mathrm{ion}}$,
where $c_{\mathrm{ion}}$ is the ion concentration. By contrast, if ion-ion
interactions were dominant, we would expect
$k_{\mathrm{crash}}\propto(c_{\mathrm{ion}})^{2}$, as discussed in the Theory
section.

From the trajectories of the bulk NaCl solutions, we determined the diffusion
coefficients $D_{\mathrm{3D}}$ of water. We computed the diffusion
coefficients with an optimal Generalized Least Squares (GLS) estimator 23, as
implemented in the DiffusionGLS 24 package. The trajectories were unwrapped
using a scheme that correctly takes into account volume fluctuations in the
NPT ensemble 25, as implemented in qwarp 26.

We evaluated the properties of the membrane system using the same observables
as in the validation of the original Martini 3 force field: area-per-lipid
$A_{l}$, thickness $d$ (both of them obtained with FATSLiM27), area
compressibility modulus (computed from the projected area fluctuations
$K_{\mathrm{A}}$, without correction for finite size effect28), and $S_{n}$
order parameter (using do-order-gmx5.py, available on the Martini website
www.cgmartini.nl). In addition, we also calculated the density profile of the
ions along the membrane normal and the lateral diffusion coefficient
$D_{\mathrm{lat}}$ of the lipid centers of mass within the membrane.

##  4 Results and Discussion

###  4.1 Stability of Martini 3 NaCl solutions as a function of time step

We performed extensive MD simulations of the NaCl solutions at time steps
between 20 and 40 fs and a range of ion concentrations. By computing the rate
of crashing $k_{\mathrm{crash}}$ using Eq. (6), we found more frequent crashes
at longer time steps and with increasing number of ions, indicating an issue
with the stability of the numerical time integration associated with ions of
bead type TQ5. To describe the time step dependence of these crashes, we
fitted Eq. (5) to the $k_{\mathrm{crash}}$ values. As shown in Fig. 2, the
slope of the curves in the semi-logarithmic representation is constant.
Accordingly, a single $\gamma$ parameter accurately captures the time-step
dependence of the crash tendency independent of ion concentration.

Figure 2: Rate of crashing in the NaCl solutions as a function of the time
step size. Shown is $k_{\text{crash}}$ on a semi-logarithmic scale as a
function of $1/\Delta t^{2}$. The individual colors correspond to fixed
concentrations. The concentrations increase from the bottom curve (0.075 M) to
the top (0.5 M). The solid lines represent fits of Eq. (5) using a common
$\gamma$ value.

According to the theory, the prefactor $k_{\mathrm{crash}}^{0}=f\gamma$ in Eq.
(5) accounts for the ion concentration dependence. Indeed, the values of the
dimensionless concentration factor $f=k_{\mathrm{crash}}^{0}/\gamma$ obtained
from the intercept of the fits in Fig. 2 exhibit a linear-quadratic
concentration dependence on ion concentration, as shown in Fig. 3. We
incorporated this concentration dependence and the dependence on the overall
system volume $V$ into Eq. (5) in the form
$k_{\mathrm{crash}}^{0}\left(c_{\mathrm{ion}},V\right)=\gamma\left(ac_{\mathrm{ion}}^{2}+bc_{\mathrm{ion}}\right)V$.
The time step, concentration, and volume dependent crash rate then becomes

$k_{\mathrm{crash}}\left(c_{\mathrm{ion}},V,\Delta
t\right)=\gamma\left(ac_{\mathrm{ion}}^{2}+bc_{\mathrm{ion}}\right)V\;\exp{\left(-\frac{1}{2\gamma^{2}\Delta
t^{2}}\right)},$ (7)

where $a=7.5820\times 10^{-7}$ mM-2 nm-3 and $b=3.2043\times 10^{-4}$ mM-1
nm-3 and $\gamma=5.2568$ ps-1. Note that a rearrangement of the prefactor
gives $k_{\mathrm{crash}}^{0}\left(c_{\mathrm{ion}},V\right)=\gamma
n_{\mathrm{ions}}\left(ac_{\mathrm{ion}}+b\right)$, where
$n_{\mathrm{ions}}/2$ is the total number of positive and negative TQ5 ions in
the system.

Figure 3: Dependence of rate of crashing on the ion concentration. Shown is
the dimensionless concentration factor
$f=k_{\mathrm{crash}}^{0}/\gamma=\left(ac_{\mathrm{ion}}^{2}+bc_{\mathrm{ion}}\right)V$,
$c_{\mathrm{ion}}=2c_{\mathrm{NaCl}}$ extracted from the intercept of the
linear fits in Fig. 2 as a function of the ion concentration in NaCl water
solutions. The solid line shows a fit of a linear-quadratic dependence on
$c_{\mathrm{ion}}$ with a fixed value of $\gamma$.

Equation (5) with a linear-quadratic concentration dependence accurately
accounts for the observed crash rates $k_{\mathrm{crash}}$ across ion
concentrations and time steps. As shown in Fig. 4, we achieved an excellent
agreement between the numerical data for time steps between 29–32 fs and the
simple 3-parameter theory. As a final validation, we varied the system size.
As predicted by the theory, $k_{\mathrm{crash}}$ depends linearly on the
overall volume of the system (see Fig. LABEL:fig:si-crash-volume).
Incorporating the concentration and system size dependence into Eq. (5) gives
a general equation, Eq. (7), capable of predicting the expected number of
crashes. The analysis presented here required a careful control of certain
simulation parameters. In particular, we found that using a fixed value of
verlet-buffer-tolerance instead of a fixed neighborlist cut-off (rlist)
effectively masked the functional form of the concentration dependence of
$k_{\mathrm{crash}}$ (see Fig. LABEL:fig:si-crash-rlist).

Figure 4: Rate of crashing in the NaCl solutions as a function of NaCl
concentration analyzed for different time steps. The panels correspond to
$\Delta t=29,30,31$, and 32 fs. The solid black lines are fits of Eq. (6) with
three parameters $\gamma$, $a$ and $b$. The error bars represent the standard
error. Note that the $y$-scale varies between the panels and that
$c_{\mathrm{ion}}=2c_{\mathrm{NaCl}}$.

Besides estimating the crash rate, our model also allows to determine the
particle types involved in the crashes. At physiological ion concentrations,
$c_{\mathrm{NaCl}}=c_{\mathrm{ion}}/2\approx 150$ mM, crashes are more likely
to be caused by ion-water collisions than ion-ion collisions. According to the
theory, their ratio is about $ac_{\mathrm{ion}}/b\approx 0.71$. Supported by
the data, the theory also predicts an extremely steep increase in
$k_{\mathrm{crash}}$ as a function of the time step, which was also the main
reason for the use of such a seemingly narrow range of $\Delta t$ between
29–32 fs . With relatively short ($10^{8}$ steps) runs of a small ($\sim$6000
particles) test system, the theory predicts that one would have to perform
$\sim$10 000 seconds of simulation time at $\Delta t=20$ fs to reach a single
failure even in the most concentrated (500 mM) system. On the other extreme of
the scale, $\Delta t=40$ fs results in a failure after every $\sim$500 steps
(every 20 ps), on average. These numbers are supported by our observations in
simulations using the respective time steps. As a practical example, we
applied Eq. (7) to predict $k_{\mathrm{crash}}$ in a system of 150 mM NaCl
solution consisting of $10^{5}$ particles corresponding to a cubic box of edge
length approximately 23 nm. The estimated $k_{\mathrm{crash}}$ as a function
of the time step is presented in Fig. 5. Using $\Delta t=25$ fs, one can
expect one crash about every 366 µs of simulation time, while $\Delta t=26$ fs
already results in a crash every 41 µs. We expect that the addition of
proteins or a lipid bilayer to the systems does not fundamentally alter
$k_{\mathrm{crash}}$, provided they do not introduce additional numerical
instabilities. As a consequence, our scaling law limits the largest accessible
time step to $\Delta t\leq 25$ fs in any practical application.

Figure 5: Rate of crashing in a 150 mM NaCl solution consisting of $10^{5}$ CG
beads, as estimated using Eq. (7). Horizontal lines indicate
$k_{\text{crash}}$ for runs with time steps of $\Delta t=25$, 26, and 27 fs
(bottom to top). The lines are labeled by the expected mean time to a first
crash.

###  4.2 Increased ionic mass stabilizes time integration

In numerical tests of the NaCl solutions, we found that an increase in the
mass of the ions described by TQ5 beads resulted in more stable time
integration. According to Eq. (1), we expect that increasing the ionic masses
by a factor of $(m^{\prime}/m)^{2}=40^{2}/30^{2}\approx 1.7778$ will allow us
to use a time step longer by a factor $4/3$. In the more elaborate model Eq.
(5), the ionic mass $m$ enters through $\gamma$. For a modified ionic mass
$m^{\prime}$, we have $\gamma^{\prime}=\gamma\sqrt{m/m^{\prime}}$. In
practice, the mass $m^{\prime}$ may have to be adjusted to account for the
more complex coupled motions of the different particles. As shown in Fig.
LABEL:fig:si-crash-mass, we obtain excellent agreement of the crash rate
predicted by Eq. (5) with the data obtained for the heavier ions by using an
effective mass of $m^{\prime}=1.7188m$ instead of the nominal ionic mass of
$1.7778m$, which amounts to a difference of around 3.5%.

These findings further confirm that the ions of TQ5 type indeed determine the
stability of the NaCl solutions. Moreover, they also demonstrate that an ionic
mass increase improves the stability by allowing a time step scaled
approximately by the square root of ionic masses after and before mass
correction, in close agreement with the prediction of Eq. (1). Taking the 0.3
M test system as an example, doubling the masses at $\Delta t=30$ fs reduces
the estimated crash rate $k_{\mathrm{crash}}$ from $\approx 10/\mbox{\textmu
s}$ to $\approx 10^{-8}/\mbox{\textmu s}$, while quadrupling the masses at
$\Delta t=40$ fs reduces the estimated $k_{\mathrm{crash}}$ from $\approx
10^{5}/\mbox{\textmu s}$ to $\approx 10^{-10}/\mbox{\textmu s}$.

###  4.3 Increased ionic masses leave the structure of NaCl solutions
unchanged and minimally impact their dynamics

To investigate the perturbations caused by an increase of ionic masses, we
analyzed the structure and dynamics of the NaCl solutions. As expected, the
modification of the ionic masses does not impact the structural properties
(Fig. LABEL:fig:si-bulk-rdf). We found the ion-ion radial distribution
functions (RDFs) to be independent of the integration time step and the ionic
mass. The diffusion coefficients $D_{\mathrm{3D}}$ of CG water obtained from
these systems are collected in Table 4.3. The data confirm our expectations
based on hydrodynamic considerations 13, 14: changing the mass of the ions has
only small effects on the diffusion coefficient. At a NaCl concentration of
300 mM, a four-fold increase in ionic mass reduced the water diffusivity by
about 2%. We expect that ionic mass will impact the librational motion;
however, these non-diffusive motions occur on timescales below the frequency
of sample collection (Fig. LABEL:fig:si-solvent-gls). The observed decrease in
diffusivity with increasing NaCl concentration is in agreement with the
tendency observed in experiments29. However, we also noticed a mass-
independent small decrease in $D_{\mathrm{3D}}$ with increasing time step
$\Delta t$. When compared to the conventional approach of fitting a straight
line to an ad hoc interval of the MSD curve, the $D_{\mathrm{3D}}$ values
calculated with a GLS estimator show smaller errors, more consistent estimates
between replicas and clearer tendencies across the simulations (see Table
LABEL:tab:si-D3d along with Figs. LABEL:fig:si-solvent-gls and LABEL:fig:si-
dcoeff-estimator of the Supporting Information). Therefore, care must be taken
not only when one compares diffusion coefficients from simulations performed
using different time steps, but also when using different estimators 23.

Table 1: Estimated $D_{\mathrm{3D}}$ diffusion coefficient [nm2/ns] of the CG
water beads in Martini 3 NaCl solutions at various ion concentrations. Results
are listed for ionic masses multiplied by factors 1, 2, and 4; and for time
steps of 20, 30, and 40 fs. The errors are smaller than 0.01 nm2/ns. These
values were obtained from the correctly unwrapped trajectories 25 using
DiffusionGLS 23, 24. Entries are missing where simulations suffered from a
high rate of crashing $k_{\mathrm{crash}}$. Results are listed without
corrections for finite system size15.

4 | 0.0 M | 0.075 M | 0.15 M |  0.3 M
---|---|---|---|---
| 1$\times$m | 1$\times$m | 2$\times$m | 4$\times$m | 1$\times$m | 2$\times$m | 4$\times$m | 1$\times$m | 2$\times$m | 4$\times$m
4 20 fs | 2.43 | 2.26 | 2.25 | 2.25 | 2.11 | 2.10 | 2.08 | 1.81 | 1.80 | 1.77
30 fs | 2.28 | 2.15a | 2.14 | 2.14 | — | 2.00 | 2.00 | — | 1.73 | 1.71
40 fs | 2.10 | — | 1.99a | 1.98 | — | — | 1.86 | — | — | 1.62
4 |  |  |  |  |  |  |  |  |  |

a These simulations had $k_{\mathrm{crash}}$ around
$\mathcal{O}\left(1/\mathrm{simulation}\right)$

###  4.4 Ionic mass increase has no discernible effect on the structure and
dynamics in lipid membrane systems

Having established the lack of appreciable effects of ionic mass in bulk
aqueous solution, we turned our attention to the neat DPPC bilayers. As for
NaCl solutions, missing entries in all subsequent tables indicate that a high
crash rate made it impossible to obtain converged results. The average values
of area-per-lipid $A_{l}$ and membrane thickness $d$ are collected in Tables 2
and 3. DPPC lipids modelled with the standard Martini 3 model (denoted
M3(-TQ5)) have a smaller $A_{l}$ and correspondingly higher $d$ than in
Martini 2 (M2), consistent with its parametrization. Doubling the mass of TQ5
(system M3-2$\times$m) and subsequently increasing $\Delta t$ has no
detectable effect on either of these quantities.

Table 2: Average area-per-lipid $A_{l}$ [Å2] from simulations of DPPC bilayers. The error is less than $\pm 0.5$ $\textup{\AA}^{2}$. Entries are missing where simulations failed to run properly. M2 and M3 refer to Martini 2 and Martini 3, respectively. 4 | M2 | M3(-TQ5) | M3-SQ5 | M3-RQ5 | M3-2$\times$m
---|---|---|---|---|---
4 20 fs | 60.6 | 59.7 | 59.7 | 59.8 | 59.7
30 fs | 60.5 | — | 59.5 | 59.6 | 59.5
40 fs | 60.4 | — | — | 59.4 | —
4 |  |  |  |  |
Table 3: Membrane thickness $d$ [nm] from simulations of DPPC bilayers. The error is less than 0.025 nm. Entries are missing where simulations failed to run properly.. 4 | M2 | M3(-TQ5) | M3-SQ5 | M3-RQ5 | M3-2$\times$m
---|---|---|---|---|---
4 20 fs | 4.147 | 4.195 | 4.194 | 4.190 | 4.195
30 fs | 4.151 | — | 4.202 | 4.198 | 4.203
40 fs | 4.157 | — | — | 4.209 | —
4 |  |  |  |  |

The lipid models in Martini 3 are designed to be softer than those in Martini
2 3. Correspondingly, DPPC lipids in the Martini 3 validation exhibited
$K_{\mathrm{A}}$ values of 232.9$\pm$13.3 mN/m in excellent agreement with the
experimentally measured 231$\pm$20.0 mN/m30. The $K_{\mathrm{A}}$ values
calculated in this work are collected in Table 4. Our values are also in
general agreement with the experimental data, but somewhat below them. The
discrepancy is not surprising: $K_{\mathrm{A}}$ values are both sensitive to
the length of the time steps31 and to system size 31, 32, consistent with our
simulations. In particular, the dependence on system size originates from the
difference between real and projected membrane areas due to the undulations of
the membrane28. A word of caution is needed at this point: a common practice
when simulating large patches of planar membranes is to apply weak harmonic
restraints to a quarter of the lipids in one of the two leaflets33, 16. Such
“pinning” of a subset of lipids suppresses long wavelength undulations. The
suppression of undulations should be taken into account when comparing
$K_{\mathrm{A}}$ from simulations of different sizes. As we saw with $A_{l}$
and $d$, doubling the ionic masses does not influence $K_{\mathrm{A}}$.
However, doubling $\Delta t$ from 20 to 40 fs resulted in a noticeable
decrease (10-20%) of $K_{\mathrm{A}}$ for the M2 and M3-RQ5 systems, where
runs with $\Delta t=40$ fs were stable.

Table 4: Area compressibility modulus $K_{\mathrm{A}}$ [mN/m] from simulations of DPPC bilayers. The error is less than 2 mN/m. Entries are missing where simulations failed to run properly. Results are listed without corrections for finite size28. 4 | M2 | M3(-TQ5) | M3-SQ5 | M3-RQ5 | M3-2$\times$m
---|---|---|---|---|---
4 20 fs | 279 | 220 | 221 | 216 | 218
30 fs | 254 | — | 197 | 197 | 204
40 fs | 243 | — | — | 179 | —
4 |  |  |  |  |

The lipid tail order parameter values $S_{n}$ presented in Table 5 show a
higher degree of order in Martini 3 than in the previous version. Because
$S_{n}$ is computed as an average over the bonds involving lipid tail beads,
the values are unaffected by changes in the surrounding medium. An increase in
$\Delta t$ causes a minute increase in $S_{n}$, which matches our expectation
of straighter lipid chains based on the computed membrane thickness values.

Table 5: Average order parameter $S_{n}$ of the carbon chains from simulations of DPPC bilayers. The error is less than 0.001. Entries are missing where simulations failed to run properly. 4 | M2 | M3(-TQ5) | M3-SQ5 | M3-RQ5 | M3-2$\times$m
---|---|---|---|---|---
4 20 fs | 0.455 | 0.481 | 0.481 | 0.479 | 0.481
30 fs | 0.457 | — | 0.485 | 0.483 | 0.485
40 fs | 0.460 | — | — | 0.488 | —
4 |  |  |  |  |
Figure 6: Number density profile of ions along the direction normal to the
membrane. The mid-plane of the bilayer is located at $z=0$ nm. Results are
shown for a time step $\Delta t=20$ fs, and are representative of all time
steps.

In addition to the properties used in the optimization of the Martini 3 force
field, we calculated the density profile of the ions along the membrane normal
and lateral diffusion coefficient $D_{\mathrm{lat}}$ of the lipid centers of
mass. Similar to the radial distribution profiles of the NaCl solutions (Fig.
LABEL:fig:si-bulk-rdf), Fig. 6 demonstrates that the obtained density profiles
are insensitive to changes in the ionic mass, which is also true for the
separate density profiles of Na+ and Cl- (see Fig. LABEL:fig:si-membrane-
densprof-separate).

Compared to Martini 2, the lateral diffusion coefficient $D_{\mathrm{lat}}$ of
DPPC lipids is decreased in Martini 3 consistent with the higher order of the
lipid chains, as seen from Table 6. The GLS estimator again provided diffusion
coefficients with smaller statistical errors and clearer tendencies (see Table
LABEL:tab:si-Dlat and Fig. LABEL:fig:si-m3-gls of the Supporting Information).
Note, however, that these diffusion coefficients were not corrected for
finite-size effects 16. From the computed $D_{\mathrm{lat}}$ values it is
clear that simply changing the mass of the TQ5 bead without increasing $\Delta
t$ does not alter the diffusion coefficient of the lipids, which is in accord
with our results for the NaCl solutions.

Table 6: Lateral diffusion coefficients $D_{\mathrm{lat}}$ [10-2 nm2/ns] from simulations of DPPC bilayers. The error is smaller than 0.1$\times$10-2 nm2/ns. Entries are missing where simulations failed to run properly. 4 | M2 | M3(-TQ5) | M3-SQ5 | M3-RQ5 | M3-2$\times$m
---|---|---|---|---|---
4 20 fs | 5.80 | 5.58 | 5.64 | 5.64 | 5.54
30 fs | 5.69 | — | 5.54 | 5.49 | 5.46
40 fs | 5.51 | — | — | 5.31 | —
4 |  |  |  |  |

###  4.5 Bead type changes have minimal impact on the structure and dynamics
of the membrane

As an alternative solution to changing the masses, one can modify the bead
type used to represent ions in Martini 3. The use of SQ5 beads allowed us to
increase $\Delta t$ to 30 fs, meaning not a single crash was detected in the
40 replicas of the NaCl solutions. However, the combination of SQ5 beads and
$\Delta t=40$ fs resulted in an average of 5 crashes per simulation. Further
increasing the bead size to RQ5 made it possible to run simulations at $\Delta
t=40$ fs. The procedure of increasing the bead size is fully consistent with
the Martini 3 philosophy, and it represents the inclusion of hydrating waters
in the ionic bead. Just as in the case of mass changes, we assessed the
differences that occur upon changing the bead type from TQ5 to SQ5 or RQ5.

The average values of $A_{l}$, $d$, and $K_{\mathrm{A}}$ are listed in Tables
2, 3, and 4, respectively, in columns M3-SQ5 and M3-RQ5. Similarly to changing
the ionic mass, modifying the bead type while keeping $\Delta t$ constant has
virtually no effect on these values. Whereas increasing $\Delta t$ does not
affect $A_{l}$ and $d$ values, we found the $K_{\mathrm{A}}$ values to
decrease noticeably, by 10-20%. Nevertheless, the changes introduced between
Martini 2 and Martini 3 are largely preserved independent of the time step.
Moreover, the lipid tail order parameters $S_{n}$ in Table 5 are insensitive
to the ionic bead type, while an increase in $\Delta t$ produces an increase
comparable to that observed in the system with doubled masses (also in Table
5).

The modification of the bead type altered the distribution of ions around the
membrane. However, as Fig. 6 shows, these changes are minor compared to those
between the Martini 2 and 3 versions and between different atomistic force
fields 34. The alteration of the bead type also produced a detectable increase
in the lipid diffusivity, indicating a somewhat more fluid liquid phase (Table
6). However, the effects of changing the time step $\Delta t$ on lipid
diffusion are bigger than the effects of changing the ionic mass in the ranges
considered.

Finally, to assess the impact of changing the ionic mass or the bead type in
systems containing charged lipids, we also tested DPPC bilayers containing 5%
PI-(4,5)P2 (phosphatidylinositol 4,5-bisphosphate) 35 using the Martini 3
force field. The results are presented in Tables LABEL:tab:si-pip2-apl to
LABEL:tab:si-pip2-dlat and Figs. LABEL:fig:si-pip2-densprof and LABEL:fig:si-
pip2-densprof-separate. The changes in the computed quantities were
commensurate with those presented here for neat DPPC membranes. However, the
PI-(4,5)P2-containing simulations were limited to a maximum 30 fs time step.

##  5 Conclusions

The recently developed Martini 3 force field has successfully addressed
several fundamental issues raised by Alessandri et al.36 Most importantly, the
introduction of specific cross-interaction terms between particles of
different sizes created a consistent framework for the use of different bead
sizes.

Here, we showed that the introduction of the “tiny” Martini 3 bead type for
representing unhydrated Na+ and Cl- ions limits the accessible time step to
below 25 fs. To achieve this, we performed extensive statistical analyses of
the rate of crashing and developed a quantitative model of the crash rate
$k_{\mathrm{crash}}\left(c_{\mathrm{ion}},\Delta t\right)$. The model revealed
the role of ion-ion and ion-water collisions in bringing about the crashes.
The insight into the factors limiting the stability of MD time integration for
neat NaCl aqueous solutions could be transferred to larger and more complex
systems containing lipid bilayers. The knowledge of such stability limits
greatly facilitates the rational design of computational experiments and the
optimal use of available resources. A particular example is the free energy
method of Lechner et al. relying on fast-switching trajectories, where optimum
efficiency of the algorithm was achieved using $\Delta t$ just short of the
stability limit. 12, 37 Another example is the setup of in silico high
throughput campaigns such as the optimization of compositions of ionic
liquids38 or deep eutectic solvents39 for liquid-liquid extractions, where the
statistical model can be applied to determine an optimal time step for a given
set of molecular systems.

Moreover, our simulations of NaCl solutions demonstrated that increasing the
ionic mass has no significant effects on the structure or dynamics of the
system (beyond the timescale of librational motions). Doubling (quadrupling)
the masses allowed us to perform simulations at $\Delta t=30$ fs (40 fs), but
this resulted in somewhat slower diffusion. The properties of membrane systems
containing DPPC lipids were also insensitive to doubling the mass. The larger
time step had a minute effect on the structure of the bilayer, and the change
in $D_{\mathrm{lat}}$ was even smaller than in case of the NaCl solutions.

As an alternative consistent with the philosophy of Martini 3, we explored the
impact of changing the size of the charged bead used to represent ions. Again,
we noticed only a minor influence on the dynamical and structural properties
of protein-less bilayers, while making possible the use of an increased time
step.

Increased ionic masses or altered bead types allowed us to use time steps of
$\Delta t=30$ fs (with SQ5) or $40$ fs (with RQ5) to model protein-less
bilayers. Increased time steps can provide crucial speed-up, e.g., for
simulating phase separating lipid mixtures 40, 5. Although simulations of neat
DPPC bilayers were running stably even at 40 fs, we do not recommend going
beyond 30 fs. Major reasons to avoid time steps as long as $\Delta t=40$ fs,
besides the general consensus11, are the following: (i) the structural and
dynamical properties of the system undergo relatively larger changes between
30 and 40 fs than between 20 and 30 fs as seen, e.g., in the diffusion
coefficients of water in NaCl solutions; (ii) in the simulations of neat DPPC
bilayers, only the M3-RQ5 system could tolerate such a high time step; (iii)
the introduction of more finely mapped lipids such as phosphatidylinositol
4,5-bisphosphate35 precludes the use of 40 fs (see Supporting Information).
Our analysis indicated that to reach $\Delta t=30$ fs it suffices to either
double the ionic mass or change the bead type to SQ5, which both have only a
mild effect on the behavior of the system.

Another known issue that limits the time step in three-component phase
separating systems is the presence of insufficiently converged constraints 41.
We plan to address this problem in a future publication.

B.F. is grateful for Dr. Jakob Tómas Bullerjahn for fruitful discussions
regarding statistical data analysis. This research was supported by the Max
Planck Society (B.F. and G.H). B.F., S.T. and G.H. thank the Center for
Multiscale Modelling in Life Sciences (CMMS) sponsored by the Hessian Ministry
of Science and Art for funding. S.T. thanks the Alfons und Gertrud Kassel
foundation and the Dr. Rolf M. Schwiete foundation.

See the supplementary material for the definition of crash in simulations, the
influence of rlist, doubling of masses, and changing the volume on
$k_{\mathrm{crash}}$, the diffusion coefficients obtained with the
conventional estimator and its comparison to the GLS estimator, the ion-ion
radial distribution functions in the NaCl solutions, the density profiles of
ions around the DPPC membrane, and the results concerning the PI-(4,5)P2
containing membranes.

## References

  * Marrink et al. 2004 Marrink, S. J.; De Vries, A. H.; Mark, A. E. Coarse grained model for semiquantitative lipid simulations. _The Journal of Physical Chemistry B_ 2004, _108_ , 750–760
  * Marrink et al. 2007 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. _The journal of physical chemistry B_ 2007, _111_ , 7812–7824
  * Souza et al. 2021 Souza, P. C.; Alessandri, R.; Barnoud, J.; Thallmair, S.; Faustino, I.; Grünewald, F.; Patmanidis, I.; Abdizadeh, H.; Bruininks, B. M.; Wassenaar, T. A., et al. Martini 3: a general purpose force field for coarse-grained molecular dynamics. _Nature methods_ 2021, _18_ , 382–388
  * Khan et al. 2020 Khan, H. M.; Souza, P. C.; Thallmair, S.; Barnoud, J.; De Vries, A. H.; Marrink, S. J.; Reuter, N. Capturing choline–aromatics Cation- $\pi$ interactions in the MARTINI force field. _Journal of chemical theory and computation_ 2020, _16_ , 2550–2560
  * Weiner and Feigenson 2018 Weiner, M. D.; Feigenson, G. W. Presence and role of midplane cholesterol in lipid bilayers containing registered or antiregistered phase domains. _The Journal of Physical Chemistry B_ 2018, _122_ , 8193–8200
  * Rosetti et al. 2017 Rosetti, C. M.; Montich, G. G.; Pastorino, C. Molecular insight into the line tension of bilayer membranes containing hybrid polyunsaturated lipids. _The Journal of Physical Chemistry B_ 2017, _121_ , 1587–1600
  * Perlmutter and Sachs 2011 Perlmutter, J. D.; Sachs, J. N. Interleaflet interaction and asymmetry in phase separated lipid bilayers: molecular dynamics simulations. _Journal of the American Chemical Society_ 2011, _133_ , 6563–6577
  * De Jong et al. 2016 De Jong, D. H.; Baoukina, S.; Ingólfsson, H. I.; Marrink, S. J. Martini straight: Boosting performance using a shorter cutoff and GPUs. _Computer Physics Communications_ 2016, _199_ , 1–7
  * Xing et al. 2009 Xing, C.; Ollila, O. S.; Vattulainen, I.; Faller, R. Asymmetric nature of lateral pressure profiles in supported lipid membranes and its implications for membrane protein functions. _Soft Matter_ 2009, _5_ , 3258–3261
  * Jalili and Akhavan 2009 Jalili, S.; Akhavan, M. A coarse-grained molecular dynamics simulation of a sodium dodecyl sulfate micelle in aqueous solution. _Colloids and Surfaces A: Physicochemical and Engineering Aspects_ 2009, _352_ , 99–102
  * Marrink et al. 2010 Marrink, S. J.; Periole, X.; Tieleman, D. P.; de Vries, A. H. Comment on “On using a too large integration time step in molecular dynamics simulations of coarse-grained molecular models” by M. Winger, D. Trzesniak, R. Baron and WF van Gunsteren, Phys. Chem. Chem. Phys., 2009, 11, 1934. _Physical Chemistry Chemical Physics_ 2010, _12_ , 2254–2256
  * Lechner et al. 2006 Lechner, W.; Oberhofer, H.; Dellago, C.; Geissler, P. L. Equilibrium free energies from fast-switching trajectories with large time steps. _The Journal of chemical physics_ 2006, _124_ , 044113
  * Einstein 1905 Einstein, A. Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. _Annalen der physik_ 1905, _4_
  * Saffman and Delbrück 1975 Saffman, P.; Delbrück, M. Brownian Motion in Biological Membranes. _Proc. Natl. Acad. Sci. USA_ 1975, _72_ , 3111–3113
  * Yeh and Hummer 2004 Yeh, I.-C.; Hummer, G. System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions. _The Journal of Physical Chemistry B_ 2004, _108_ , 15873–15879
  * Vögele et al. 2018 Vögele, M.; Köfinger, J.; Hummer, G. Hydrodynamics of diffusion in lipid membrane simulations. _Physical review letters_ 2018, _120_ , 268104
  * Wassenaar et al. 2015 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. _Journal of chemical theory and computation_ 2015, _11_ , 2144–2155
  * Abraham et al. 2015 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_ , 19–25
  * Bussi et al. 2007 Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. _The Journal of chemical physics_ 2007, _126_ , 014101
  * Parrinello and Rahman 1981 Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. _Journal of Applied physics_ 1981, _52_ , 7182–7190
  * Humphrey et al. 1996 Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. _Journal of Molecular Graphics_ 1996, _14_ , 33–38
  * Berendsen et al. 1984 Berendsen, H. J.; Postma, J. v.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. _The Journal of chemical physics_ 1984, _81_ , 3684–3690
  * Bullerjahn et al. 2020 Bullerjahn, J. T.; von Bülow, S.; Hummer, G. Optimal estimates of self-diffusion coefficients from molecular dynamics simulations. _The Journal of Chemical Physics_ 2020, _153_ , 024116
  * 24 Diffusion Coefficient Fitting. https://github.com/bio-phys/DiffusionGLS
  * von Bülow et al. 2020 von Bülow, S.; Bullerjahn, J. T.; Hummer, G. Systematic errors in diffusion coefficients from long-time molecular dynamics simulations at constant pressure. _The Journal of Chemical Physics_ 2020, _153_ , 021101
  * 26 qwrap — Fast PBC wrapping and unwrapping for VMD. https://github.com/jhenin/qwrap
  * Buchoux 2017 Buchoux, S. FATSLiM: a fast and robust software to analyze MD simulations of membranes. _Bioinformatics_ 2017, _33_ , 133–134
  * Waheed and Edholm 2009 Waheed, Q.; Edholm, O. Undulation contributions to the area compressibility in lipid bilayer simulations. _Biophysical journal_ 2009, _97_ , 2754–2760
  * Ben Ishai et al. 2013 Ben Ishai, P.; Mamontov, E.; Nickels, J. D.; Sokolov, A. P. Influence of Ions on Water Diffusion—A Neutron Scattering Study. _The Journal of Physical Chemistry B_ 2013, _117_ , 7724–7728
  * Nagle and Tristram-Nagle 2000 Nagle, J. F.; Tristram-Nagle, S. Structure of lipid bilayers. _Biochimica et Biophysica Acta (BBA)-Reviews on Biomembranes_ 2000, _1469_ , 159–195
  * Venable et al. 2015 Venable, R. M.; Brown, F. L.; Pastor, R. W. Mechanical properties of lipid bilayers from molecular dynamics simulation. _Chemistry and physics of lipids_ 2015, _192_ , 60–74
  * Saeedimasine et al. 2019 Saeedimasine, M.; Montanino, A.; Kleiven, S.; Villa, A. Role of lipid composition on the structural and mechanical features of axonal membranes: a molecular simulation study. _Scientific reports_ 2019, _9_ , 1–12
  * Ingólfsson et al. 2014 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. _Journal of the american chemical society_ 2014, _136_ , 14554–14559
  * Catte et al. 2016 Catte, A.; Girych, M.; Javanainen, M.; Loison, C.; Melcr, J.; Miettinen, M. S.; Monticelli, L.; Määttä, J.; Oganesyan, V. S.; Ollila, O. S., et al. Molecular electrometer and binding of cations to phospholipid bilayers. _Physical Chemistry Chemical Physics_ 2016, _18_ , 32560–32569
  * Borges-Araújo et al. 2021 Borges-Araújo, L.; Souza, P.; Fernandes, F.; Melo, M. N. Improved parameterization of phosphatidylinositide lipid headgroups for the Martini 3 coarse grain force field. 2021,
  * Alessandri et al. 2019 Alessandri, R.; Souza, P. C.; Thallmair, S.; Melo, M. N.; De Vries, A. H.; Marrink, S. J. Pitfalls of the Martini model. _Journal of chemical theory and computation_ 2019, _15_ , 5448–5460
  * Oberhofer and Dellago 2007 Oberhofer, H.; Dellago, C. Large time-step, fast-switching free energy calculations with non-symplectic integrators. _Israel Journal of Chemistry_ 2007, _47_ , 215–223
  * Vazquez-Salazar et al. 2020 Vazquez-Salazar, L. I.; Selle, M.; De Vries, A. H.; Marrink, S. J.; Souza, P. C. Martini coarse-grained models of imidazolium-based ionic liquids: from nanostructural organization to liquid–liquid extraction. _Green Chemistry_ 2020, _22_ , 7376–7386
  * Vainikka et al. 2021 Vainikka, P.; Thallmair, S.; Souza, P. C. T.; Marrink, S. J. Martini 3 Coarse-Grained Model for Type III Deep Eutectic Solvents: Thermodynamic, Structural, and Extraction Properties. _ACS Sustainable Chemistry & Engineering_ 2021, _9_ , 17338–17350
  * Baoukina et al. 2017 Baoukina, S.; Rozmanov, D.; Tieleman, D. P. Composition fluctuations in lipid bilayers. _Biophysical journal_ 2017, _113_ , 2750–2761
  * Thallmair et al. 2021 Thallmair, S.; Javanainen, M.; Fábián, B.; Martinez-Seara, H.; Marrink, S. J. Nonconverged Constraints Cause Artificial Temperature Gradients in Lipid Bilayer Simulations. _The Journal of Physical Chemistry B_ 2021, _125_ , 9537–9546