GROMACS version:2020.1
GROMACS modification: No
Dear Sir
I am absolutely new to GROMACS simulation and I am completely lost on Glycine solvation energy in water. I used OPLSaa.ff/forcefields (I also tried charm27/forcefield.itp in separate experiment) in 3x3x3 nm^3 cube filled with water molecule. The density is close to 1g/1ml. First, I minimize the energy, then equilibrate in NVT followed by NPT for 100ps. Then, I used 10 steps to scale lamda for vdw followed by 4 step Coulomb in production run. In all experiments, I have solvation energy of 200 to 300 kj/mol where I think I should get 40kj/mol. What am I doing wrong? I have tried several TOP3P, TIP4P and SPCE. I have also tried gmx VS gmx_d, on VS off DPROSRES, charge on each atoms in Glycine, but results did not change that much. Below are cut off of top and mdp files. The itp file refered in top file is an exact copy of charm27/forcefield.itp.
-----top------
; Include forcefield parameters
#include “Charmforcefield.itp”
[ moleculetype ]
; Name nrexcl
SOUT 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
; residue 1 GLY rtp GLY q 0.0
1 NH3 1 GLY N 1 -0.3 14.007
2 HC 1 GLY H1 2 0.33 1.008
3 HC 1 GLY H2 3 0.33 1.008
4 HC 1 GLY H3 4 0.33 1.008
5 CT2 1 GLY CA 5 0.13 12.011
6 HB 1 GLY HA1 6 0.09 1.008
7 HB 1 GLY HA2 7 0.09 1.008
8 CC 1 GLY C 8 0.34 12.011
9 OC 1 GLY OT1 9 -0.67 15.9994
10 OC 1 GLY OT2 10 -0.67 15.9994 ; qtot 0
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 3 1
1 4 1
1 5 1
5 6 1
5 7 1
5 8 1
8 9 1
8 10 1
[ pairs ]
; ai aj funct c0 c1 c2 c3
1 9 1
1 10 1
2 6 1
2 7 1
2 8 1
3 6 1
3 7 1
3 8 1
4 6 1
4 7 1
4 8 1
6 9 1
6 10 1
7 9 1
7 10 1
[ angles ]
; ai aj ak funct c0 c1 c2 c3
2 1 3 5
2 1 4 5
2 1 5 5
3 1 4 5
3 1 5 5
4 1 5 5
1 5 6 5
1 5 7 5
1 5 8 5
6 5 7 5
6 5 8 5
7 5 8 5
5 8 9 5
5 8 10 5
9 8 10 5
[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3 c4 c5
2 1 5 6 9
2 1 5 7 9
2 1 5 8 9
3 1 5 6 9
3 1 5 7 9
3 1 5 8 9
4 1 5 6 9
4 1 5 7 9
4 1 5 8 9
1 5 8 9 9
1 5 8 10 9
6 5 8 9 9
6 5 8 10 9
7 5 8 9 9
7 5 8 10 9
[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3
8 5 10 9 2
; Include Position restraint file
#ifdef POSRES
#include “posre.itp”
#endif
; Include water topology
#include “charmm27.ff/tip4p.itp”
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
; Include topology for ions
#include “charmm27.ff/ions.itp”
[ system ]
; Name
GlyZ in water
[ molecules ]
; Compound #mols
SOUT 1
SOL 1338
----------mdp file------
title = Gly-Z FE
;define = -DPOSRES ; position restrain the protein
integrator = sd
nsteps = 50000
dt = 0.002
nstenergy = 10 ;originally 1000 (gromacs free energy tutorial pdf) changed to 100 to match with nstdhdl
nstlog = 5000
; cut-offs at 1.0nm
rlist = 1.0
dispcorr = EnerPres
vdw-type = pme
rvdw = 1.0
; Coulomb interactions
coulombtype = pme
rcoulomb = 1.0
fourierspacing = 0.12
; Constraints
constraints = all-bonds
; set temperature to 300K
tcoupl = v-rescale
tc-grps = system
tau-t = 0.2
ref-t = 300
; set pressure to 1 bar with a thermostat that gives a correct
; thermodynamic ensemble
pcoupl = parrinello-rahman
ref-p = 1
compressibility = 4.5e-5
tau-p = 5
; and set the free energy parameters
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; FREE ENERGY CONTROL (this one I got it from lidcaine example online)
free-energy = yes
; Standard soft-core potential parameters
sc-power = 1
sc-sigma = 0.3
sc-alpha = 1.0
; The molecule being coupled/decoupled
couple-moltype = SOUT ; special molecule
; Sub-ensemble of current simulation
init-lambda-state = 0
; Defining our target and reference state.
; My preference is for lambda = 0 to correspond to the IG state,
; where intermolecular interactions are off and for
; lambda = 1 to correspond to the target state, where intermolecular
; interactions are at full strength.
;
; Defining lambda = 0 (ideal gas state)
couple-lambda0 = none
; Defining lambda = 1 (fully interacting state)
couple-lambda1 = vdw-q
; Keep intramolecular interactions ON in all states
couple-intramol = no
; Lambda values of each sub-ensemble
; Same linear scaling sub-ensembles 0-5 as with alkanes.
; For charges we will add subensembles 6-9 where charges are scaled as sqrt(1/m_q)
; m = 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
vdw-lambdas = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.00 1.00 1.00 1.00
coul-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.71 0.87 1.00
; Print the differnces in Hamiltonians between all lambda values for MBAR
calc-lambda-neighbors = -1
; Frequency at which difference in the Hamiltonians is calculated
; (Set to same as nstenergy)
nstdhdl = 10 ; 0.20 ps
------end--------
Thank you very much and I truly appreciate your help on this.