Input parameters

Which mdp options should I use?

See the example mdp file for recommended parameter values, and the answers to specific parameters below.

How large can the time step be?

Martini has been parameterized using time steps in the range of 20-40 fs. Whether you can use 40 fs or have to settle for a somewhat smaller time step depends on your system, and on your attitude toward coarse-grained modeling, as explained below.

First, the Martini force field is not an atomistically detailed force field. Many assumptions underlie the model, the major one being the neglect of some of the atomistic degrees of freedom. As a result, the interactions between particles are effective ones and the energy landscape is highly simplified. This simplified energy landscape allows for a greatly increased sampling speed at the cost of a loss of detail. This makes CG models in general so powerful. The emphasis, therefore, should not be to sample the energy landscape as accurately as possible, but rather, as effectively as possible. This is in contrast to traditional all-atom models, for which the energy landscape is more realistic and an accurate integration scheme is more important. In practice, the inherent ‘fuzziness’ of the Martini model makes the presence ofartificial small energy sinks or sources a less critical problem than in accurate atomistic simulations. Second, most importantly, structural properties are rather very robust with respect to time step; For a time step of 40 fs, there are no noticeable effects on structural properties of the systems investigated. Moreover, thermodynamic properties such as the free energy of solvation also appear insensitive to the size of the time step. Thus, if the goal is to generate representative ensembles quickly, large time steps seem acceptable.

Based on the two arguments discussed above, we conclude that time steps exceeding 10 fs can be used in the Martini force field. Whereas one can debate the first argument (i.e. the ‘idealist’ versus ‘pragmatic’ view of the power of CG simulations), the second argument (i.e. the insensitivity of both structural and thermodynamic properties to the magnitude of the time step) implies that a reduction of the time step to 10 fs or below is a waste of computer time. Nevertheless, time steps of 40 fs and beyond may be pushing the limits too far for certain systems. We therefore recommend a time step of 20-30 fs, in combination with an enlarged neighbourlist cut-off (to 1.4 nm) to be on the safe side.

Of course, one should always check whether or not results are biased by the choices made. Given that the largest simplifications are made at the level of the interaction potentials, this can best be done by comparing to results obtained using more detailed models.

Please see the following two papers for a recent discussion on the use of large time steps in coarse-grained simulations:

[1] M.Winger, D. Trzesniak, R. Baron, W.F. van Gunsteren. On using a too large integration time step in molecular dynamics simulations of coarse-grained molecular models, Phys. Chem. Chem. Phys., 2009, 11, 1934-1941.

[2] S.J. Marrink, X. Periole, D.P. Tieleman, A.H. de Vries. Comment on using a too large integration time step in molecular dynamics simulations of coarse-grained molecular models. Phys. Chem. Chem. Phys., 12:2254-2256, 2010.

Which coupling time should I use for temperature and pressure control?

Good temperature control can be achieved with the Berendsen thermostat, using a coupling constant of the order of τ = 1 ps. Even better temperature control can be achieved by reducing the temperature coupling constant to 0.1 ps, although with such tight coupling (τ approaching the time step) one can no longer speak of a weak-coupling scheme.

Similarly, pressure can be controlled with the Berendsen barostat, with a coupling constant in the range 1-5 ps and typical compressibility in the order of 10-4 - 10-3 bar-1. Note that, in order to estimate compressibilities from CG simulations, you should use Parrinello-Rahman type coupling.

How often do I need to update the pairlist and how large should my pairlist cutoff be?

Due to the use of shifted potentials, the noise generated form particles leaving/entering the neighbour list is not so large, even when large time steps are being used. In practice, once every ten steps works fine with a neighborlist cutoff that is equal to the non-bonded cutoff (1.2 nm). However, to improve energy conservation or to avoid local heating/cooling, you may increase the update frequency and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option is computationally less expensive and leads to improved energy conservation (see [1]).

[1] S.J. Marrink, X. Periole, D.P. Tieleman, A.H. de Vries. Comment on using a too large integration time step in molecular dynamics simulations of coarse-grained molecular models, Phys. Chem. Chem. Phys, 12:2254-2256, 2010).

Which cutoffs should I use?

Standard cut-off schemes are used for the non-bonded interactions in the Martini model: LJ interactions are shifted to zero in the range 0.9-1.2 nm, and electrostatic interactions in the range 0.0-1.2 nm. For details about the shift function, see the answer to the question below. The treatment of the non-bonded cut-offs is considered to be part of the force field parameterization, so we recommend not to touch these values as they will alter the overall balance of the force field.

Which shift function is used?

The shift function Φ is the following, as implemented in Gromacs:

shifted-potentials

Here α denotes the power of the respective LJ (α=6,12) or Coulombic (α=1) terms. Both potential and force are continuous and smoothly decay to zero between rshift and the cut-off distance rcut. The Martini force field has been parameterized with rshift=0.9 (LJ) or 0.0 (Coulomb) and rcut=1.2 nm.

Can I use PME or reaction field with Martini?

In principle you can, although the Martini force field has been parameterized with short range shifted electrostatic interactions. The use of reaction field (which is effectively also a shifted potential) is not likely to affect the behavior very much. PME, on the other hand, may lead to significantly different behavior, and could be more realistic in certain applications (e.g. realistic membrane pores were seen in simulations of both dendrimers [1] and anti-microbial peptides [2] attacking lipid bilayers). Please realize that electrostatic interactions in the Martini model are not considered to be very accurate to begin with, especially as the screening in the system is set to be uniform across the system with a screening constant of 15. When using PME, please make sure your system properties are still reasonable.

However, in combination with the polarizable Martini water model [3], PME can be used and makes sense as the elctrostatic interactions are more realistic.

[1] H. Lee and R. G. Larson, J. Phys. Chem. B, 2008, 112, 7778–7784

[2] A.J. Rzepiela, D. Sengupta, N. Goga, S.J. Marrink, Faraday Discuss., 2010, 144, 431-443.

[3]  S.O. Yesylevskyy, L.V. Schäfer, D. Sengupta, S.J. Marrink, PLoS Comp. Biol, 6:e1000810, 2010. open access

Can I do free energy calculations with Martini?

Yes you can!