Hello everyone, I am doing free energy calculations upon aminoacid mutation for the first time and got a problems. I want to calculate equilibrium trajectories for fast-growth method.
Here is what I did.
- Downloaded PDB structure from RCSB
- Splitted it into several chains(because I dont know how to rebuild missing residues in in multimeric protein)
- Rebuilt missing residues in every chain using Modeller
- Combined all chains together
- Created hybrid topology using PMX
- Added box, water and ions in Gromacs
During minimization step there was 1 atom that Gromacs tried to minimize 90% of time.
Then, after gmx grompp -f nvt.mdp -c em.gro -r em.gro -p hybrid.top -o nvt.tpr -maxwarn 2
I’ve got 2 warnings
WARNING 1 [file hybrid_posre_Protein_chain_A.itp, line 2692]:
Some parameters for bonded interaction involving perturbed atoms are
specified explicitly in state A, but not B - copying A to B
WARNING 2 [file hybrid.top, line 19]:
The bond in molecule-type Protein_chain_A between atoms 5324 DOH and 5325
DHH has an estimated oscillational period of 6.5e-03 ps, which is less
than 5 times the time step of 2.0e-03 ps.
Maybe you forgot to change the constraints mdp option.
And then after gmx mdrun -deffnm nvt
Fatal error:
There are perturbed non-bonded pair interactions beyond the pair-list cutoff
of 1.111 nm, which is not supported. This can happen because the system is
unstable or because intra-molecular interactions at long distances are
excluded. If the latter is the case, you can try to increase nstlist or rlist
to avoid this.The error is likely triggered by the use of couple-intramol=no
and the maximal distance in the decoupled molecule exceeding rlist.
Could you please help me to solve this issue. It seems for me(maybe I am wrong) that the source of the problem are dummy atoms, because without hybrid topology(without dummy atoms) I dont have any problems. Maybe dummy atoms make clashes with another residues. If so, is there right way how to construct hybrid topology without clashes? Thanks.
Here are my mdp files
minimization.mdp
;Run parameters
define = -DFLEXIBLE
nsteps = 1500
integrator = steep
emtol = 50.0
emstep = 0.01
;Neighbour searching
cutoff-scheme = Verlet
ns-type = grid
rcoulomb = 1.0
vdwtype = cut-off
rvdw = 1.0
nstlist = 10
rlist = 1
;Eletrostatics
coulombtype = PME
; FREE ENERGY
;----------------------------------------------------
free-energy = yes
init-lambda = 0 ; start from state A
delta-lambda = 0 ; stay in state A
;
nstcomm = 100
constraints = none
coulomb-modifier = Potential-shift-Verlet
pme-order = 4
fourierspacing = 0.12
ewald-rtol = 1e-5
vdw-modifier = Potential-shift-Verlet
verlet-buffer-tolerance = 0.005
DispCorr = EnerPres
Tcoupl = no
Pcoupl = no
gen_vel = no
nvt.mdp
;Type of MDP: nvt
;Position restrain
Define = -DPOSRES
;Run parameters
nsteps = 50000
integrator = sd
dt = 0.002
comm-mode = Linear ; remove center of mass translation
;Output control
nstxout = 500
nstvout = 500
nstenergy = 500
nstlog = 500
nstcalcenergy = 100
nstcomm = 100
nstxout-compressed = 1000
compressed-x-precision = 1000
compressed-x-grps = System
;Bond parameters
constraint-algorithm = lincs
constraints = h-bonds
lincs-iter = 1
lincs-order = 4
continuation = no
;Neighbour searching
cutoff-scheme = Verlet
ns-type = grid
rcoulomb = 1.0
vdwtype = cut-off
rvdw = 1.0
nstlist = 10
rlist = 1
;Eletrostatics
coulombtype = PME
pme-order = 4
fourierspacing = 0.12
fourier-nx = 0
fourier-ny = 0
fourier-nz = 0
ewald-rtol = 1e-5
;Temperature coupling
tcoupl = V-rescale
tc-grps = Protein Non-Protein
tau-t = 0.1 0.1
ref-t = 300 300
;Pressure coupling
pcoupl = no
;Dispersion correction
DispCorr = EnerPres
;Velocity generation
gen-vel = yes
gen-temp = 300
gen-seed = -1
;Periodic boundary conditions
pbc = xyz
ld-seed = 17
; FREE ENERGY
;----------------------------------------------------
free-energy = yes
init-lambda = 0 ; start from state A
delta-lambda = 0 ; stay in state A
;
vdw-modifier = Potential-shift-Verlet