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

Question

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?:

Answer 1

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

Answer 2

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

df = pd.read_csv("/home/nick/GROMACS/martini3/default-water-ff-params.txt", delim_whitespace=True, names=['bead-1', 'bead-2', 'x', 'sigma', 'epsilon'])

for index, row in df.iterrows():
    
    if 'W' not in row['bead-2']:
        
        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.