# Unable to recreate data for benzene from paper in GROMACS 2020.5

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 is 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


My nvt equilibration mdp:

title                   = liquid benzene nvt equilibration
; define                  = -DPOSRES  ; position restrain the protein
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 50000     ; 2 * 50000 = 100 ps
dt                      = 0.002     ; 2 fs
; Output control
nstxout                 = 500       ; save coordinates every 1.0 ps
nstvout                 = 500       ; save velocities every 1.0 ps
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
; Bond parameters
continuation            = no        ; first dynamics run
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
; Nonbonded settings
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = System                ; one coupling groups - more accurate
tau_t                   = 0.1               ; time constant, in ps
ref_t                   = 298               ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl                  = no        ; no pressure coupling in NVT
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Velocity generation
gen_vel                 = yes       ; assign velocities from Maxwell distribution
gen_temp                = 300       ; temperature for Maxwell distribution
gen_seed                = -1        ; generate a random seed



My npt equilibration file:

title                   = liquid-benzene-simulation
; define                  = -DPOSRES  ; position restrain the protein
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 100000     ; 1 * 100000 = 100 ps
dt                      = 0.001     ; 1 fs
; Output control
nstxout                 = 500       ; save coordinates every 1.0 ps
nstvout                 = 500       ; save velocities every 1.0 ps
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
; Bond parameters
continuation            = yes       ; Restarting after NVT
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
; Nonbonded settings
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; 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.1                   ; 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
refcoord_scaling        = com
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Velocity generation
gen_vel                 = no        ; Velocity generation is off

$$$$
`
• One thing I noticed is that you have 600 benzene molecules in the topology file, but the .gro file only has one benzene molecule defined. Did you create a new .gro file with more than one benzene? Mar 3 at 14:53
• I can't edit my old comment now, but after digging up the paper, I can see that there are some differences in the cutoffs and other options in your md.mdp file. The paper used a cutoff of 1.5, you used 1.3 nm, the paper used Nose-Hoover thermostat, you used modified Berendsen, and there are some other minute differences. Would it be possible for you to post a link to the actual .gro file that you used for the NVT or NPT steps, so that everything is clear? Mar 3 at 15:15
• right, yes, I can do that. The thing I posted was my initial gro file, which I editconf-ed, then I applied insert-molecule to get the 600 molecules. Mar 3 at 15:40
• initial 600 molecule gro file: ctxt.io/2/AACgtax-FA, 600 molecule gro file after NVT: ctxt.io/2/AACgtax-Eg, 600 molecule gro file after NPT: ctxt.io/2/AACgpV-XEg Mar 3 at 15:47
• @ShoubhikRMaiti, I agree about the 5nm. I simply did that to ensure that all the molecules fit in the box, since i was going to equilibrate anyway with NPT later. But I will change the bond-constraints Mar 3 at 17:46