Unable to recreate data for benzene from paper in GROMACS 2020.5

GROMACS version: 2020.5
GROMACS modification:

I have been trying to simulate benzene in OPLSS-AA forcefield as per this paper: https://pubs.acs.org/doi/10.1021/ct2002122. I have been at a loss, since I believe I have everything running as per the paper, except for the vdw and coulombic cutoffs, but I am not generating the results the paper is showing.

My simulation is running bug-free on GROMACS, however, I am not getting the results they are showing. I energy minimized my system, performed an NVT equilibration for 50 ps, then I performed an NPT equilibration for a 100 ps. Then I did a production run for 5 ns.

I have rcoulomb = rvdw = 1.3 in my simulation. That is because GROMACS does not allow twin cutoff ranges anymore.

My density is coming out in the neighborhood of 823 kg/m3 when it should be around 867 kg/m3.

I obtained my GRO and ITP file from LigParGen server, however, the charges on the LigParGen Server were different from the charges in the paper, so I changed them.

The only thing the paper hasn’t provided are the harmonic bond and angle potential coefficients, and torsional coefficients.

I am really stuck as to how to go about this problem. Where could my bug be lying?



For the sake of completeness, I have attached my primary input files. This is my .gro file (coordinate):

LIGPARGEN GENERATED GRO FILE
   12
    1UNK    C00    1   0.100   0.100   0.000
    1UNK    C01    2  -0.040   0.100   0.000
    1UNK    C02    3  -0.109   0.100   0.121
    1UNK    C03    4  -0.039   0.100   0.242
    1UNK    C04    5   0.100   0.100   0.241
    1UNK    C05    6   0.170   0.100   0.121
    1UNK    H06    7   0.154   0.100  -0.094
    1UNK    H07    8  -0.094   0.100  -0.094
    1UNK    H08    9  -0.218   0.100   0.121
    1UNK    H09   10  -0.093   0.100   0.336
    1UNK    H0A   11   0.155   0.100   0.335
    1UNK    H0B   12   0.278   0.100   0.121
   1.00000   1.00000   1.00000

This my .itp file:

;
; GENERATED BY LigParGen Server
; Jorgensen Lab @ Yale University
;
[ atomtypes ]
  opls_804  C804    12.0110     0.000    A    3.55000E-01   2.92880E-01
  opls_803  C803    12.0110     0.000    A    3.55000E-01   2.92880E-01
  opls_811  H811     1.0080     0.000    A    2.42000E-01   1.25520E-01
  opls_809  H809     1.0080     0.000    A    2.42000E-01   1.25520E-01
  opls_810  H810     1.0080     0.000    A    2.42000E-01   1.25520E-01
  opls_805  C805    12.0110     0.000    A    3.55000E-01   2.92880E-01
  opls_802  C802    12.0110     0.000    A    3.55000E-01   2.92880E-01
  opls_806  H806     1.0080     0.000    A    2.42000E-01   1.25520E-01
  opls_808  H808     1.0080     0.000    A    2.42000E-01   1.25520E-01
  opls_800  C800    12.0110     0.000    A    3.55000E-01   2.92880E-01
  opls_807  H807     1.0080     0.000    A    2.42000E-01   1.25520E-01
  opls_801  C801    12.0110     0.000    A    3.55000E-01   2.92880E-01
[ moleculetype ]
; Name               nrexcl
UNK                   3
[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass
     1   opls_800      1    UNK   C00      1    -0.115    12.0110
     2   opls_801      1    UNK   C01      1    -0.115    12.0110
     3   opls_802      1    UNK   C02      1    -0.115    12.0110
     4   opls_803      1    UNK   C03      1    -0.115    12.0110
     5   opls_804      1    UNK   C04      1    -0.115    12.0110
     6   opls_805      1    UNK   C05      1    -0.115    12.0110
     7   opls_806      1    UNK   H06      1     0.115     1.0080
     8   opls_807      1    UNK   H07      1     0.115     1.0080
     9   opls_808      1    UNK   H08      1     0.115     1.0080
    10   opls_809      1    UNK   H09      1     0.115     1.0080
    11   opls_810      1    UNK   H0A      1     0.115     1.0080
    12   opls_811      1    UNK   H0B      1     0.115     1.0080
[ bonds ]
    2     1     1      0.1400 392459.200
    3     2     1      0.1400 392459.200
    4     3     1      0.1400 392459.200
    5     4     1      0.1400 392459.200
    6     1     1      0.1400 392459.200
    7     1     1      0.1080 307105.600
    8     2     1      0.1080 307105.600
    9     3     1      0.1080 307105.600
   10     4     1      0.1080 307105.600
   11     5     1      0.1080 307105.600
   12     6     1      0.1080 307105.600
    6     5     1      0.1400 392459.200

[ angles ]
;  ai    aj    ak funct            c0            c1            c2            c3
    1     2     3     1    120.000    527.184
    2     3     4     1    120.000    527.184
    3     4     5     1    120.000    527.184
    2     1     6     1    120.000    527.184
    2     1     7     1    120.000    292.880
    1     2     8     1    120.000    292.880
    2     3     9     1    120.000    292.880
    3     4    10     1    120.000    292.880
    4     5    11     1    120.000    292.880
    1     6    12     1    120.000    292.880
    6     5    11     1    120.000    292.880
    6     1     7     1    120.000    292.880
    5     6    12     1    120.000    292.880
    5     4    10     1    120.000    292.880
    1     6     5     1    120.000    527.184
    4     3     9     1    120.000    292.880
    4     5     6     1    120.000    527.184
    3     2     8     1    120.000    292.880

[ dihedrals ]
; IMPROPER DIHEDRAL ANGLES
;  ai    aj    ak    al funct            c0            c1            c2            c3            c4            c5
    12     6     1     5    4        180.000     10.460     2
    11     5     4     6    4        180.000     10.460     2
    10     4     3     5    4        180.000     10.460     2
     9     3     2     4    4        180.000     10.460     2
     8     2     1     3    4        180.000     10.460     2
     7     1     2     6    4        180.000     10.460     2

[ dihedrals ]
; PROPER DIHEDRAL ANGLES
;  ai    aj    ak    al funct            c0            c1            c2            c3            c4            c5
    4    3    2    1        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    5    4    3    2        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    6    1    2    3        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    4    5    6    1        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    6    5    4    3        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    5    6    1    2        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
   12    6    1    2        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    7    1    6    5        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
   11    5    6    1        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
   10    4    5    6        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
   11    5    4    3        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
   12    6    5    4        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    9    3    2    1        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    7    1    2    3        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    8    2    1    6        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
   10    4    3    2        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    8    2    3    4        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    9    3    4    5        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
   10    4    3    9        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
   12    6    1    7        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
   12    6    5   11        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    9    3    2    8        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
    8    2    1    7        3      30.334   0.000  -30.334  -0.000  -0.000   0.000
   11    5    4   10        3      30.334   0.000  -30.334  -0.000  -0.000   0.000

[ pairs ]
     1     4    1
     2     5    1
     3     6    1
     3     7    1
     1     9    1
     5     7    1
     4     8    1
     2    10    1
     1    11    1
     6     8    1
     5     9    1
     3    11    1
     2    12    1
     7     8    1
     6    10    1
     4    12    1
     8     9    1
     9    10    1
     7    12    1
    10    11    1
    11    12    1

This is my topol.top file:

; include forcefield parameters
#include "oplsaa.ff/forcefield.itp"

; Include mol topology
#include "UNK_5D626B_rev.itp"

[ system ]
; Name
liquid-benzene-simulation

[ molecules ]
;compound   number of molecules
UNK        600

And this is my final production run md.mdp file:

title                   = liquid benzene oplss simulation
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 5000000    ; 1 * 1000000 = 1000 ps (1 ns)
dt                      = 0.001     ; 1 fs
; Output control
nstxout                 = 0         ; suppress bulky .trr file by specifying
nstvout                 = 0         ; 0 for output frequency of nstxout,
nstfout                 = 0         ; nstvout, and nstfout
nstenergy               = 5000      ; save energies every 5.0 ps
nstlog                  = 5000      ; update log file every 5.0 ps
nstxout-compressed      = 5000      ; save compressed coordinates every 5.0 ps
compressed-x-grps       = System    ; save the whole system
; Bond parameters
continuation            = yes       ; Restarting after NPT
constraint_algorithm    = lincs     ; holonomic constraints
constraints             = h-bonds   ; bonds involving H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Neighborsearching
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb                = 1.3       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.3       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 8         ; order-8 interpolation
fourierspacing          = 0.12      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = System                ; two coupling groups - more accurate
tau_t                   = 0.5                   ; time constant, in ps
ref_t                   = 298                   ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
pcoupltype              = isotropic             ; uniform scaling of box vectors
tau_p                   = 2.0                   ; time constant, in ps
ref_p                   = 1.0                   ; reference pressure, in bar
compressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Dispersion correction
DispCorr                = Ener  ; account for cut-off vdW scheme
; Velocity generation
gen_vel                 = no        ; Velocity generation is off

Please see my answer to your previous post
\Alessandra