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