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

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.


Introduction
Compared to its predecessors, 1,2 the recent Martini 3 force field 3 constitutes a significant advance in biomolecular simulations.Martini 3 offers consistently parametrized coarsegrained (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-π type) can be incorporated using labels which were previously only available in the polarizable version. 4The 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 ∆t = 30 fs [5][6][7] according to the recommendation of the new-rf parameter set, 8 the current recommended time step is ∆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 ∆t = 40 fs was used. 2,9,10However, 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 ∆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.

Theory 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 ∆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. 12r 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 α and √ α, 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 where we ignored the motions of the other particles.Then, starting from a time step ∆t = 20 fs, doubling and quadrupling the ionic mass m would make it possible to use ∆t = 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. 14Both 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,16dification 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.

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 ∆t be the integration time step.With v the velocity in one dimension, the corresponding displacement is then ∆x = v∆t, ignoring the higher-order acceleration terms.We now assume that displacements ∆x > ∆x crit are critical, and lead to crashes.The velocities of a particle of mass m at inverse temperature β satisfy a Maxwell-Boltzmann distribution, p(v)dv = exp(−βmv 2 /2) βm/(2π)dv.The probability of a critical displacement during a simulation time step becomes with 1/γ 2 = βm∆x 2 crit .Here we used the asymptotic expansion of the complementary error function.In a simulation system, we expect that a certain fraction √ 2πf of the particle displacements (here, those of the "tiny" ions and their interaction partners) can result in such catastrophic clashes.The probability displacement in a given time step is approximately Assuming uncorrelated events and 0 < Q 1, the probability that a simulation run of n time steps will not crash is with t = n∆t the total time.The crash times are thus predicted to be distributed exponentially, p(t crash )dt crash = k crash exp(−k crash t crash )dt crash with a rate that depends on the simulation time step as The crash rate has two parameters, f and γ, whose product defines k 0 crash = f γ.In our systems, we expect that the dimensionless factor f depends on the concentration c ion of ions as f = (ac ion 2 + bc ion )V with coefficients a and b, and V the system volume.Based on the definition, and consistent with Eq. ( 1), we expect γ to depend on the mass of the ions approximately as γ ∝ 1/ √ m.

Methods Simulation of NaCl solutions
To demonstrate the effect of changing ionic masses (or the lack thereof) on the viscosity µ 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.py 17and 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 field 3 and the Gromacs package 18 (version 2020.1).To eliminate differences stemming from finite size effects between these systems, 15 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×m, 2×m, 4×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. 19The pressure was maintained at 1 bar by an isotropic Parrinello-Rahman barostat 20 (p ref = 1 bar and τ 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 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 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.

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.py. 17All bilayers were constructed in an initial box of 18×18×18 nm 3 and solvated by 14405 CG waters (pure W for Martini 3, 10% WF antifreeze 2 for Martini 2).The equilibrium membrane area in the simulation boxes was about 17x17 nm 2 , larger than the ≈8x8 nm 2 used in the Martini 3 validation test.Every simulation system contained 0.15 M NaCl, as added with insane.pyusing 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×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 thermostat. 19The pressure of the systems was controlled in a semi-isotropic manner with a compressibility of 3 × 10 Figure 1 shows illustrative snapshots of the simulated NaCl solutions and lipid bilayers.

Analysis
In our numerical tests, we ran 40 simulations for each system.The time of t i is the targeted run length, t i = ∆t × 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, where 0 ≤ n crash ≤ 40 is the number of crashes observed in the 40 runs.We assessed the uncertainty of the estimator using the standard error σ = k crash / √ n 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 crash ∝ c ion , where c ion is the ion concentration.By contrast, if ion-ion interactions were dominant, we would expect k crash ∝ (c ion ) 2 , as discussed in the Theory section.
From the trajectories of the bulk NaCl solutions, we determined the diffusion coefficients D 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 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 FATSLiM 27 ), area compressibility modulus (computed from the projected area fluctuations K A , without correction for finite size effect 28 ), 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 lat of the lipid centers of mass within the membrane.

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 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 crash values.As shown in Fig. 2, the slope of the curves in the semilogarithmic representation is constant.Accordingly, a single γ parameter accurately captures the time-step dependence of the crash tendency independent of ion concentration.According to the theory, the prefactor k 0 crash = f γ in Eq. ( 5) accounts for the ion concentration dependence.Indeed, the values of the dimensionless concentration factor f = k 0 crash /γ 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 0 crash (c ion , V ) = γ (ac 2 ion + bc ion ) V .The time step, concentration, and volume dependent crash rate then becomes where a = 7.5820×10 −7 mM −2 nm −3 and b = 3.2043×10 −4 mM −1 nm −3 and γ = 5.2568 ps −1 .
Note that a rearrangement of the prefactor gives k 0 crash (c ion , V ) = γn ions (ac ion + b), where n ions /2 is the total number of positive and negative TQ5 ions in the system.Equation ( 5) with a linear-quadratic concentration dependence accurately accounts for the observed crash rates k 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 crash depends linearly on the overall volume of the system (see Fig. S1).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 crash (see Fig. S2).Besides estimating the crash rate, our model also allows to determine the particle types involved in the crashes.At physiological ion concentrations, c NaCl = c ion /2 ≈ 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 ion /b ≈ 0.71.Supported by the data, the theory also predicts an extremely steep increase in k crash as a function of the time step, which was also the main reason for the use of such a seemingly narrow range of ∆t between 29-32 fs .With relatively short (10 8 steps) runs of a small (∼6000 particles) test system, the theory predicts that one would have to perform ∼10 000 seconds of simulation time at ∆t = 20 fs to reach a single failure even in the most concentrated (500 mM) system.On the other extreme of the scale, ∆t = 40 fs results in a failure after every ∼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 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 crash as a function of the time step is presented in Fig. 5. Using ∆t = 25 fs, one can expect one crash about every 366 µs of simulation time, while ∆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 crash , provided they do not introduce additional numerical instabilities.As a consequence, our scaling law limits the largest accessible time step to ∆t ≤ 25 fs in any practical application.

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 /m) 2 = 40 2 /30 2 ≈ 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 γ.For a modified ionic mass m , we have γ = γ m/m .In practice, the mass m may have to be adjusted to account for the more complex coupled motions of the different particles.As shown in Fig. S3, 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 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 ∆t = 30 fs reduces the estimated crash rate k crash from ≈ 10/µs to ≈ 10 −8 /µs, while quadrupling the masses at ∆t = 40 fs reduces the estimated k crash from ≈ 10 5 /µs to ≈ 10 −10 /µs.

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. S7).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 3D of CG water obtained from these systems are collected in Table 1.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. S4).
The observed decrease in diffusivity with increasing NaCl concentration is in agreement with the tendency observed in experiments. 29However, we also noticed a mass-independent small decrease in D 3D with increasing time step ∆t.When compared to the conventional approach of fitting a straight line to an ad hoc interval of the MSD curve, the D 3D values calculated with a GLS estimator show smaller errors, more consistent estimates between replicas and clearer tendencies across the simulations (see Table S1 along with Figs.S4 and S6 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. 23ble 1: Estimated D 3D diffusion coefficient [nm 2 /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 nm 2 /ns.These values were obtained from the correctly unwrapped trajectories 25 using DiffusionGLS. 23,24Entries are missing where simulations suffered from a high rate of crashing k crash .Results are listed without corrections for finite system size.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×m) and subsequently increasing ∆t has no detectable effect on either of these quantities. in excellent agreement with the experimentally measured 231±20.0mN/m. 30The K 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 A values are both sensitive to the length of the time steps 31 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 membrane. 28A 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 leaflets. 16,33Such "pinning" of a subset of lipids suppresses long wavelength undulations.The suppression of undulations should be taken into account when comparing K A from simulations of different sizes.As we saw with A l and d, doubling the ionic masses does not influence K A .However, doubling ∆t from 20 to 40 fs resulted in a noticeable decrease (10-20%) of K A for the M2 and M3-RQ5 systems, where runs with ∆t = 40 fs were stable.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 ∆t causes a minute increase in S n , which matches our expectation of straighter lipid chains based on the computed membrane thickness values.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 lat of the lipid centers of mass.Similar to the radial distribution profiles of the NaCl solutions (Fig. S7), 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. S8).
Compared to Martini 2, the lateral diffusion coefficient D 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 S2 and Fig. S5 of the Supporting Information).Note, however, that these diffusion coefficients were not corrected for finite-size effects. 16From the computed D lat values it is clear that simply changing the mass of the TQ5 bead without increasing ∆t does not alter the diffusion coefficient of the lipids, which is in accord with our results for the NaCl solutions.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 ∆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. 34The 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 ∆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)P 2 (phosphatidylinositol 4,5-bisphosphate) 35 using the Martini 3 force field.The results are presented in Tables S3 to S7 and Figs.S9 and S10.The changes in the computed quantities were commensurate with those presented here for neat DPPC membranes.However, the PI-(4,5)P 2 -containing simulations were limited to a maximum 30 fs time step.

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 crossinteraction 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 crash (c ion , ∆t).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 ∆t just short of the stability limit. 12,37Another example is the setup of in silico high throughput campaigns such as the optimization of compositions of ionic liquids 38 or deep eutectic solvents 39 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 ∆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 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 ∆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.precludes the use of 40 fs (see Supporting Information).Our analysis indicated that to reach ∆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. 41We plan to address this problem in a future publication.

Gromacs error messages
Diagnosing issues when using highly optimized MD codes is not a trivial task.With the exception of LINCS-related issues, Gromacs generally does not print out the steps leading to a crash, therefore problems are generally elusive and difficult to track down.Here, we mention the encountered error messages to provide guidance in detecting the above discussed problem.

Errors in the solvent simulations
All solvent simulations were preformed using 1 node with 1 GPU.The 3 types of error observed were the following: • the simulation stalled: the process was active, but nothing was written to the logfile.
• the simulation stopped with an error message after several Warning: Pressure scaling more than 1%.
• one update resulted in "not-a-number" (NAN), but the simulation continued until the prescribed number of timesteps.Most of the logged values and final averages were NAN.
We considered a simulation as crashed at the first occurrence of any of these error.

Errors in the membrane simulations
An invalid displacement of the ions can make the interacting partner explode.Such explosions result in the atoms being scattered way outside the simulation box leading to non-physical geometries.In the N P T ensemble, this is usually followed by a huge scaling of the simulation box and a warning from the barostat.Finally, if constraints are present, this leads to LINCS issues LINCS WARNING relative constraint deviation after LINCS ... and eventual failure One or more atoms moved too far between two domain decomposition steps.

Effect of increasing the system volume on k crash
To verify that the rate of crashing k crash is an extensive property of the system, we increased the linear system size of the NaCl solution by a factor of 2, scaling the volume by a factor of 8.The k crash values computed using ∆t = 29 fs are presented in Fig. S1, where the k crash values of the larger system were divided by 8 to account for the volume difference.

Influence of rlist on k crash
To study k crash as a function of concentration and time step, we kept constant all the other molecular dynamics settings.However, as Fig. S2 illustrates, a fixed value of verlet-buffer-tolerance resulted in a step-like change in k crash as a function of the concentration.The "inner list" value of the "equivalent classical 1x1 list", as listed in the log file, exhibited prefect correlation with the step-like behavior of k crash .Moreover, fixing the rlist value itself the mdp file eliminated the step-like changes, and we obtained smoothly varying values of k crash , as shown in the main text.

Effect of increased ion mass on crash rate
According to Eq. ( 1) in the main text, scaling the masses by ∆t 2 /∆t 2 should reduce k crash from those observed at ∆t to ∆t.We tested this relationship by increasing the masses of the ions by a factor of 40 2 /30 2 ≈ 1.7778, expecting the resulting rate to be similar to the ∆t = 30 fs simulations (see Fig. S3).The resulting k crash (red markers) was somewhere between that of the runs with unmodified ion masses and time steps of 30 and 31 fs (blue and green).A direct fit of Eq. ( 7

]
Figure S3: Rate of crashing k crash in the solvent simulations using ions with a mass scaled by a factor 40 2 /30 2 ≈ 1.7778 and a time of ∆t = 40 fs (red markers).The fit of Eq. ( 5) in the main text gave an effective mass increase of m = 1.7188m, about 3.5% lower than expected from the ratio of the nominal ion masses.For reference, the blue and green symbols show the values of k crash obtained with unmodified ion masses at time steps of 30 and 31 fs, respectively, while the continuous lines are fits of Eq. 7.
6 Diffusion coefficients and quality factors from the GLS estimator

Figure 1 :
Figure1: 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

2 ∆t − 2 [ 1 ]Figure 2 :
Figure2: Rate of crashing in the NaCl solutions as a function of the time step size.Shown is k crash on a semi-logarithmic scale as a function of 1/∆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 γ value.

Figure 3 :
Figure3: Dependence of rate of crashing on the ion concentration.Shown is the dimensionless concentration factor f = k 0 crash /γ = (ac 2 ion + bc ion ) V , c ion = 2c NaCl extracted from the intercept of the linear fits in Fig.2as a function of the ion concentration in NaCl water solutions.The solid line shows a fit of a linear-quadratic dependence on c ion with a fixed value of γ.

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

1 ]Figure 5 :
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 crash for runs with time steps of ∆t = 25, 26, and 27 fs (bottom to top).The lines are labeled by the expected mean time to a first crash.

Figure 6 :
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 ∆t = 20 fs, and are representative of all time steps.
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 ∆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 ∆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 ∆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 A are listed in Tables2, 3, and 4, respectively, in columns M3-SQ5 and M3-RQ5.Similarly to changing the ionic mass, modifying the bead type while keeping ∆t constant has virtually no effect on these values.Whereas increasing ∆t does not affect A l and d values, we found the K A values to decrease noticeably, by 10-20%.

1 ]
FigureS1: Volume-scaled rate of crashing k crash V 0 /V in the original solvent simulations of V 0 (blue markers) and in a system with volume V = 8V 0 (red markers) at ∆t = 29 fs.Note that c ion = 2c NaCl .

Figure S2 :
Figure S2: Rate of crashing in the solvent simulations as a function of concentration, at ∆t = 30 fs.The red line corresponds to the rlist simulation parameter, as described in the text.
) with variable γ = γ m/m gave an effective mass of m = 1.7188m, which is very close to the nominal ion mass of 1.7778m.These results indicate that the scaling of the masses indeed closely follows Eq. (1).c NaCl [mM] Rate of crashing k crash [µs −1

Figure
Figure S4: A representative plot of the Generalized Least Squares (GLS) estimation of the water 3D diffusion coefficient of the Martini 3 NaCl solutions.Optimal diffusion coefficients as a function of lag time (TOP).Quality factor a function of lag time (BOTTOM).

Figure S5 :
Figure S5: Generalized Least Squares (GLS) estimation of the lateral diffusion coefficient of the Martini 3(-TQ5) DPPC lipids.Optimal diffusion coefficients as a function of lag time (TOP).Quality factor as a function of lag time (BOTTOM).

Figure S6 :
Figure S6: A representative plot of the 3D diffusion coefficients of CG water beads obtained with gmx msd (blue triangles) and the Generalized Least Squares (GLS) estimator (red squares) computed over 4 replicas.The small shifts of the symbols in ∆t from 20, 30, and 40 fs are for better visibility.The dashed line connect the mean values of the respective data sets, and are guides for the eye.
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 a These simulations had k crash around O (1/simulation)

Table 2 :
Average area-per-lipid A l [ Å2 ] from simulations of DPPC bilayers.The error is less than ±0.5Å2 .Entries are missing where simulations failed to run properly.M2 and M3 refer to Martini 2 and Martini 3, respectively.

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

Table 4 :
28ea compressibility modulus K 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 size.28

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.

Table 6 :
Lateral diffusion coefficients D lat [10 −2 nm 2 /ns] from simulations of DPPC bilayers.The error is smaller than 0.1×10 −2 nm 2 /ns.Entries are missing where simulations failed to run properly.
5,40Although 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 ∆t = 40 fs, besides the general consensus, 11 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