[gmx-users] Implementing the NERD Force Field in GROMACS 2019.3

[Robert Cordina]

Hi, I’m trying to use the NERD forcefield (Sum et al. J. Phys. Chem. B 2003, 107, 14443-14451) in GROMACS 2019.3 for a pure triacylglyceride system, however I’m not entirely sure that I’m setting up the equilibration parameters in the mdp file correctly as I’m getting different simulation results from others who have published papers using this force field.  All the bonding and non-bonding interactions have been typed up in the respective force-field files, and assuming those are correct (I’ve checked and re-checked them multiple times), the only thing I can think that is not correct are the mdp parameters.

The Sum et al. paper, and another paper (Pizzirusso et al. J.Am.Chem.Soc.2018, 140, 12405−12414) state the following (text from Sum paper in italics, text from Pizzirusso paper in bold, what I’ve put in the mdp file in regular font)

Both NVT and NPT calculations were performed using 40 molecules in a cubic box with periodic boundary conditions
pbc                         = xyz

Interactions were truncated and shifted at rc = 11 Å with energies shifted at this distance
rcoulomb             = 1.1 (am I missing something here?)

The system evolved with a leapfrog algorithm using a 2 fs timestep
integrator            = md
dt                            = 0.002

Constant temperature and constant pressure simulations were maintained with the Berendsen’s thermostat and barostat, respectively, by uniform scaling of the atomic velocities (temperature) and by uniform scaling of the atomic positions and box length (pressure)
Atomistic MD was performed using GROMACS 4.5 under an NPT ensemble. In this ensemble, temperature and pressure were kept constant using a Berendsen thermostat/barostat.  Temperature and pressure couplings of 1 and 10 ps were used, respectively. Compressibility was fixed at 1 × 10-5 bar-1, and anisotropic pressure coupling was used in the MD simulations.
tcoupl                   = Berendsen
tau-t                      = 1.0
tc-grps                  = XXX
ref-t                       = 350

pcoupl                  = Berendsen
pcoupltype         = anisotropic
compressibility  = 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001
tau-p                     = 10
ref-p                      = 1.01325 1.01325 1.01325 1.01325 1.01325 1.01325

For all the calculations, a reaction-field correction was applied with continuum dielectric εRF to correct for long-range interactions due to electrostatics
cutoff-scheme  = Verlet
coulombtype     = Reaction-Field
epsilon-rf            = a.bcd (I used actual numbers here of course)


In addition I am using
gen-vel                 = yes
gen-temp            = 350
gen-seed             = -1
continuation      = no


Should I be adding anything with respect to vdw?


In the forcefield.itp file I’m setting the following
nbfunc                  = 1
comb-rule           = 2 (Sum paper states that they use Lorentz-Berthelot combining rule for the LJ parameters)
gen-pairs             = yes
fudgeLJ                = 1
fudgeQQ             = 1

In ffbonded.itp
bondtypes func = 1
angletypes func = 1
dihedraltypes func = 5



Thanks in advance for your help.

Robert