[gmx-users] Implementing the NERD Force Field in GROMACS 2019.3
Robert Cordina
robert.cordina at strath.ac.uk
Wed Feb 19 11:05:34 CET 2020
Dear Alessandra,
Many thanks for your reply.
Re comment on nrexcl, I'm using nrexcl = 3 as the Sum paper states that "In the NERD force field, atoms/sites within a molecule that do not interact by any other intramolecular potential are also allowed to interact though the Lennard-Jones potential". In the paper all bond stretching, bond bending and torsion angle parameters are given, so is nrexcl = 3 correct? LJ potentials are 1-4 interactions, so this looks the right way to do it to me, but I stand to be corrected.
Re compressibility, that was a typo from my part in the email, it's correctly set at 0.00001.
Re vdwtype, vdw-modifier, rvdw-switch, rvdw... I have read this section multiple times but I'm still not sure how to translate the description in the Sum paper to these settings. Any suggestions as to what "Interactions were truncated and shifted at rc = 11 Å with energies shifted at this distance" would mean for these settings?
cutoff-scheme = Verlet
rcoulomb = 1.1
coulombtype = Reaction-Field (I have to use this because I'm using the epsilon-rf setting too right?)
epsilon-rf = a.bcd
vdwtype = Cut-Off (?)
vdw-modifier = Potential-shift-Verlet or Potential-switch (?)
rvdw-switch = 0 (?)
rvdw = 1.1 (will this just be the same as rcoulomb?)
Thanks again for your help.
Best regards,
Robert
-----Original Message-----
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> On Behalf Of Alessandra Villa
Sent: 19 February 2020 08:47
To: gmx-users at gromacs.org
Cc: gromacs.org_gmx-users at maillist.sys.kth.se
Subject: Re: [gmx-users] Implementing the NERD Force Field in GROMACS 2019.3
HI,
Below some suggestions assuming that the energy terms are correctly implemented (force field)
On Tue, Feb 18, 2020 at 3:13 PM Robert Cordina <robert.cordina at strath.ac.uk>
wrote:
> 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
>
>
If you do not have a pre-built structure, your equilibration time may be longer that in the original paper, the best is to ask the authors the coordinate of pre-equilibrated system.
> Interactions were truncated and shifted at rc = 11 Å with energies
> shifted at this distance
> rcoulomb = 1.1 (am I missing something here?)
>
>
here you miss the vdwtype; vdw-modifier; rvdw-switch¶; rvdw
(see for definition:
http://manual.gromacs.org/documentation/current/user-guide/mdp-options.html?highlight=vdw#Van
)
> 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
>
system
> ref-t = 350
>
> pcoupl = Berendsen
> pcoupltype = anisotropic
> compressibility = 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001
>
compressibility is 10-5 bar-1
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
>
>
check that nrexcl in topol.top is in line with the force field.
> In ffbonded.itp
> bondtypes func = 1
> angletypes func = 1
> dihedraltypes func = 5
>
>
>
Best regards
Alessandra
>
> Thanks in advance for your help.
>
> Robert
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
On Tue, Feb 18, 2020 at 3:13 PM Robert Cordina <robert.cordina at strath.ac.uk>
wrote:
> 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
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
--
Gromacs Users mailing list
* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
More information about the gromacs.org_gmx-users
mailing list