Hi everyone,
Something has been bugging me for a while – I get large density errors (1.5 – 2.5 % let’s say) compared to experimental data when simulating aqueous solutions of sugars at higher temperatures (320 – 350 K). Some papers say this is normal, but some papers (one, specifically) have better predictions than I do. I ran a simulation of a solution of 50% of sucrose in water – used TIP4P water water along with the OPLS-AA force field, with the same parameters and same system that the paper I am talking about used, and got errors of around 2.5%, while the paper says their errors are of 1% at 343 K with the OPLS-AA force field. They’ve used the TIP4P/2005 water model and, when I tried using it, I got a even larger error: 4.76%. I must be doing something wrong, but I don’t know what.
These are my .mdp files:
(i) EM:
; Run control
integrator = steep
nsteps = 5000000
; EM criteria and other stuff
emtol = 100
emstep = 0.01
niter = 20
nbfgscorr = 10
; Output control
nstlog = 1
nstenergy = 1
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 1
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-05
epsilon_surface = 0
; Temperature and pressure coupling are off during EM
tcoupl = no
pcoupl = no
; Free energy control stuff
free_energy = no
; No velocities during EM
gen_vel = no
; options for bonds
constraints = h-bonds
; Do not constrain the starting configuration
continuation = no
(ii) NVT
; Run control
integrator = md
tinit = 0
dt = 0.001
nsteps = 1000000 ; 1 ns
nstcomm = 100
; Output control
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 20000
nstenergy = 20000
nstxout-compressed = 20000
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-05
epsilon_surface = 0
; Temperature coupling
Tcoupl = nose-hoover
tc_grps = system
tau_t = 1.0
ref_t = {temperature}
; Pressure coupling is off for NVT
Pcoupl = no
; Free energy control
free_energy = no
; Do NOT generate velocities — continue from EM
gen_vel = no
gen_temp = {temperature}
gen_seed = -1
; Bond constraints
constraints = h-bonds
continuation = yes
(iii) NPT
; Run control
integrator = md
tinit = 0
dt = 0.001
nsteps = 2000000 ; 2 ns
nstcomm = 100
; Output control
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 20000
nstenergy = 20000
nstxout-compressed = 20000
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-05
epsilon_surface = 0
; Temperature coupling
Tcoupl = nose-hoover
tc_grps = system
tau_t = 1.0
ref_t = {temperature}
; Pressure coupling is on for NPT
Pcoupl = Berendsen
pcoupltype = isotropic
tau_p = 2.0
compressibility = 4.5e-05
ref_p = 1.0
; Free energy control
free_energy = no
; Do NOT generate velocities — continue from NVT
gen_vel = no
gen_temp = {temperature}
gen_seed = -1
; Bond constraints
constraints = h-bonds
continuation = yes
(iv) Production
; Run control
integrator = md
tinit = 0
dt = 0.001 ; 1 fs
nsteps = 2500000 ; 2.5 ns
nstcomm = 100
; Output control
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 20000
nstenergy = 20000
nstxout-compressed = 10000
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2 ; Short-range neighbor list cutoff
; Electrostatics
coulombtype = PME
rcoulomb = 1.2 ; Short-range electrostatics cutoff (12 Å)
fourierspacing = 0.12
pme_order = 4
ewald_rtol = 1e-05
epsilon_surface = 0
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2 ; LJ cutoff at 12 Å
; Long-range dispersion correction
DispCorr = EnerPres ; Apply long-range corrections to E and P
; Temperature coupling
Tcoupl = nose-hoover
tc_grps = system
tau_t = 1.0
ref_t = {temperature}
; Pressure coupling
Pcoupl = C-rescale
pcoupltype = isotropic
tau_p = 5.0
compressibility = 4.5e-05
ref_p = 1.0
; Free energy control
free_energy = no
; Do NOT generate velocities — continue from NPT
gen_vel = no
gen_temp = {temperature}
gen_seed = -1
; Bond constraints
constraints = h-bonds
continuation = yes
What do you think I could possibly be messing up in this case?
I’ve constructed the system exactly as they did – with LigParGen. Same amount of water molecules and sucrose molecules. However, they’ve used a 12.5 A value for LJ cut-off and I’ve used 12 A because I think GROMACS doesn’t allow rvdw to be bigger than rcoulomb in this case. Also, they’ve used the particle–particle particle-mesh method, and I think it is not implemented on GROMACS. They’ve also used the Nosé-Hoover barostat, I don’t see it in the manual as an option in GROMACS. Anyways, these are the differences I could spot. Any help is appreciated :)