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-requ...@gromacs.org.
