# How can I modify the Martini v2 or v3 force field to increase the Lennard Jones potential well depth (epsilon)?

I am trying to run some coarse-grained simulations of an Intrinsically Disordered Protein (IDP) using the Martini force field, and I want to increase the strength of interaction between water and protein. Specifically I want to alter the $$\epsilon$$ value to change the Lennard-Jones potential calculation at each update:

$$V_{\text{LJ}}(r_{ij})=4\epsilon_{ij}\Bigg[\bigg(\frac{\sigma_{ij}}{r_{ij}}\bigg)^{12}-\bigg(\frac{\sigma_{ij}}{r_{ij}}\bigg)^{6}\Bigg]$$

From my understanding, the potential parameters are stored in this .itp file in the form of an interaction matrix:

martini_v3.0.b.3.2.itp:

[ nonbond_params ]
P5    P5  1 4.700000e-01    5.000000e+00
P4    P4  1 4.700000e-01    4.750000e+00
P2    P2  1 4.700000e-01    4.300000e+00
P1    P1  1 4.700000e-01    4.100000e+00
N3    N3  1 4.700000e-01    4.100000e+00
...


martini_v2.2.itp:

[ nonbond_params ]
P5    P5      1   0.24145E-00     0.26027E-02 ; supra attractive
SP5   SP5     1   0.10620E-00     0.67132E-03 ; 75supra attractive, s=0.43
P4    P4      1   0.21558E-00     0.23238E-02 ; attractive
BP4   BP4     1   0.21558E-00     0.23238E-02 ; attractive
SP4   SP4     1   0.94820E-01     0.59939E-03 ; 75attractive, s=0.43
...


There are no headers so I don't know which value is which, but I imagine the first is $$\sigma$$ and the second is $$\epsilon$$ ? I also can't find what each of the letters represents except for some definitions in the 2008 Martini paper.

Do I need to write a parser in python that modifies these column values for the relevant interaction pairs? Or is there an easier way to do this?:

• +1 Welcome to our new community and thank you for contributing your question here!!! We hope to see much more of you in the future! How did you find us? Jul 17 at 17:59
• Hi, Thank you, I just googled "GROMACS Stack Exchange" Jul 17 at 18:13

Since you seem to have at least a start on automating the process of modifying the file, I will just address the other two parts of your question: the format of the nonbond_params section of a GROMACS .itp file and the meaning of the symbols used in the Martini force field.

## GROMACS format

This is given in the GROMACS manual. To summarize:

Atom1 Atom2    Func (1)     C6/sigma  C12/epsilon


The first two columns give the types of particles interacting, the third column being 1 means it is a Lennard-Jones interaction. The last columns are the one of two different parameterizations. The interpretation of these values depends on the combination rule specified in the earlier [defaults] section of the .itp file. It determines the meaning of the values stored in the these two columns and how to generate missing cross-interactions. If this combination rule equals 1, the next two columns are \begin{align} C_{ij}^{(6)}&=4\epsilon_{ij}\sigma_{ij}^6\\ C_{ij}^{(12)}&=4\epsilon_{ij}\sigma_{ij}^{12} \end{align} And missing interactions between atoms of different types are formed by: \begin{align} C_{ij}^{(n)}&=\sqrt{C_{ii}^{(n)}C_{jj}^{(n)}} \end{align}

A value of 2 or 3 will make the next columns $$\sigma$$ and $$\epsilon$$, with different rules for creating missing interactions. Combination rule 3 will use the formula above, forming C6/C12 from $$\sigma$$ and $$\epsilon$$ while combination rule 2 uses the following formulas. \begin{align} \sigma_{ij}&=\frac{1}{2}(\sigma_{ii}\sigma_{jj})\\ \epsilon_{ij}&=\sqrt{\epsilon_{ii}\epsilon_{jj}} \end{align}

I believe the Martini values are specified as $$\sigma$$ and $$\epsilon$$ and all pairs of interactions are specified, so you should be able to use either 2 or 3 for the combination rule. This may is probably the default setting in these files, so you ideally shouldn't need to change anything.

## Martini

The atom types are described in [1][2]. The first reference is the original paper, while the 2nd is about new additions in Martini 3.

Martini 3 defines five main types of coarse grained sites: Q (charged), P (polar), N (nonpolar), C (apolar), and W (water). Each of these has further subtypes depending on what particular atoms are included in the site and their relative interaction energy. They generally map 4 atoms to a site, but the S (small) and T (tiny) prefixes denote 3:1 and 2:1 mappings for a slightly less coarse depiction. The LJ $$\sigma$$ and $$\epsilon$$ values among all these CG sites were reparameterized in Martini 3 for each possible interaction. Note, if you wind up working with Martini 2, there was no W site and instead water was included in the P group (P4 specifically).

1. Siewert J. Marrink, H. Jelger Risselada, Serge Yefimov, D. Peter Tieleman, and Alex H. de Vries The MARTINI Force Field:  Coarse Grained Model for Biomolecular Simulations The Journal of Physical Chemistry B 2007 111 (27), 7812-7824 DOI: 10.1021/jp071097f
2. Souza, P.C.T., Alessandri, R., Barnoud, J. et al. Martini 3: a general purpose force field for coarse-grained molecular dynamics. Nat Methods 18, 382–388 (2021). https://doi.org/10.1038/s41592-021-01098-3
• Thanks for the detailed answer -- I eventually came across these papers that you have cited and managed to piece together a better understanding. I still don't fully understand the format column... Do you mean that the format column allows the identity of the next two columns to serve as different parameters in say a force field that uses a potential with more/different variables? Was I wrong to simply multiply all of the $\epsilon$ values by a constant? In my martini_v3.0.0.itp force field file, all of the format' values are 1, so I didn't change these. Jul 18 at 13:42
• I believe I misinterpreted that section. I believe you are correct that this would be epsilon. I'm updating my answer.
– Tyberius
Jul 18 at 15:19
• @Vranvs I updated the answer. The value in that column seems to just specify the form of the potential, with 1 indicating LJ. I mistook it for the combination rule, which should be specified once earlier in the topology file. The Martini force field seems to be specified with sigma and epsilon, so I believe the approach described in your answer is correct and is modifying epsilon.
– Tyberius
Jul 18 at 15:38

I believe I have figured it out.

For those trying to martinize their all-atom files into coarse grain representations for martini 3, use the following package:

https://github.com/marrink-lab/vermouth-martinize

With respect to modifying the water for IDP simulations, I pulled out the W, SW and TW interactions from the martini v3.0.0 forcefield, like so:

 W    P6  1 4.650000e-01    4.750000e+00
SW    P6  1 4.350000e-01    3.880000e+00
TW    P6  1 3.950000e-01    3.410000e+00
W   SP6  1 4.250000e-01    4.530000e+00
SW   SP6  1 4.050000e-01    3.730000e+00
TW   SP6  1 3.780000e-01    3.290000e+00
W   TP6  1 3.850000e-01    4.630000e+00
SW   TP6  1 3.750000e-01    3.650000e+00
...


The last column is the epsilon value. I saved all these values into a separate file called default-water-ff-params.txt, and wrote a short script in python:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

#How much to increase the strength? (Between 1 and 1.5 typically)
scalar = 1.06 #6% increase in strength between water and protein

for index, row in df.iterrows():

df.loc[index, 'epsilon'] = row['epsilon'] * scalar

print(df.to_string(index=False))


Output from this script is like so:

bead-1 bead-2  x  sigma    epsilon
W     P6  1  0.465   5.996766
SW     P6  1  0.435   4.898411
TW     P6  1  0.395   4.305046
W    SP6  1  0.425   5.719021
SW    SP6  1  0.405   4.709039
TW    SP6  1  0.378   4.153549
W    TP6  1  0.385   5.845268
SW    TP6  1  0.375   4.608041
...
`

Make a new copy of the martini_v3.0.0.itp forcefield file and name it something like martini_v3.0.0_WATER-PROTEIN-1.06.itp to specify the increased strength of interaction. Copy the values from the script to replace the originals, and run your MD simulation.

It would be good if someone can confirm that this is the right approach. My protein (which in previous sims was very compact) quickly unfolded into a noodle, consistent with experimental data, and the $$R_g$$ fluctuated around the measured value. So, I am pretty sure this changed the forcefield parameters, but would like confirmation that I didn't modify an interaction I shouldn't have.