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$
    – S R Maiti
    Commented Mar 3, 2021 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$
    – S R Maiti
    Commented Mar 3, 2021 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
    Commented Mar 3, 2021 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
    Commented Mar 3, 2021 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
    Commented Mar 3, 2021 at 17:46

1 Answer 1

5
$\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$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .