11
$\begingroup$

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 

```
$\endgroup$
10
  • 2
    $\begingroup$ 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? $\endgroup$ – Shoubhik R Maiti Mar 3 at 14:53
  • 2
    $\begingroup$ 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? $\endgroup$ – Shoubhik R Maiti Mar 3 at 15:15
  • $\begingroup$ 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. $\endgroup$ – megamence Mar 3 at 15:40
  • $\begingroup$ 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 $\endgroup$ – megamence Mar 3 at 15:47
  • 1
    $\begingroup$ @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 $\endgroup$ – megamence Mar 3 at 17:46
4
$\begingroup$

Thanks to the generous time and comments of @ShoubhikRMaiti, I was able to recreate the simulation.

The problem in my code was 2-fold: I was using the topology filed created by LigParGen, which had different charges on the C and H atoms in my simulation as in the paper. I did not think to check LigParGen's output, and that is my own fault.

Also, I was using LINCS constraint on only the H-bonds in my simulation, when instead I was supposed to be using LINCS constraint on all-bonds in my simulation. This led to almost a perfect reproduction of the properties mentioned in the paper.

Thank you everyone who chipped in!!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.