I am trying to run a simulation of protein and ligand and I getting a LINCS error as
Step 33, time 0.033 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.002814, max 0.053588 (between atoms 1553 and 1552)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1553 1552 90.0 0.1000 0.1054 0.1000
Step 34, time 0.034 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000318, max 0.005474 (between atoms 1553 and 1552)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1553 1552 70.0 0.1054 0.1005 0.1000
Step 35, time 0.035 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000320, max 0.005247 (between atoms 1553 and 1552)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1553 1552 79.9 0.1005 0.0995 0.1000
nvt.mdp
title = 123
define = -DPOSRES -DPOSRES_LIGAND ; Run parameters
integrator = md-vv ; leap-frog integrator
nsteps = 2500000 ; 0.001 * 2500000 = 2500 ps
dt = 0.001 ; 1 fs
; Output control
nstxout = 5000 ; suppress .trr output
nstvout = 5000 ; suppress .trr output
nstenergy = 5000 ; save energies every 1000.0 ps
nstlog = 5000 ; update log file every 1000.0 ps
nstxtcout = 5000
xtc_precision = 100
energygrps = Protein LIG SOL
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
ns_type = grid ; search neighboring grid cells
nstlist = 5 ;
rcoulomb = 0.9
rvdw = 1.4
coulombtype = PME ; Particle Mesh Ewald for long-
rlistlong = 1.4
rlist = 0.9
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
periodic_molecules = yes
; Temperature coupling
tcoupl = nose-hoover ; modified Berendsen
thermostat
tc-grps = Protein_LIG Water_and_ions ; two coupling groups - more
accurate
tau_t = 2.0 2.0 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one
for each group, in K
; Pressure coupling
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
Please help in resolving this issue.
How many steps did you minimize? What was the largest force? Where are your ligand parameters from? What exactly are your position restraints? What atoms are causing the error? Ligand? Binding site residue atoms? What water model do you use?
Seriously, this error could have many reasons.
Have you tried to increase lincs warn angle to see if the error equilibrates out after a few hundreds of steps? You could also try to equilibrate only a few steps, and then minimize again, and repreat a couple of times. Stuff heavily depends on the quality of your initial structure.
Your mdp file seems ok, though.
Thanks Eric for your response.
Performed em and it converged in 169 steps.
Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 1000 (which may not be possible for your system).
It stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters.
You might need to increase your constraint accuracy, or turn
off constraints altogether (set constraints = none in mdp file)
Steepest Descents converged to machine precision in 169 steps,
but did not reach the requested Fmax < 1000.
Potential Energy = -2.20600345422483e+13
Maximum force = 3.53020624912800e+25 on atom 40
Norm of force = 5.52100420408231e+23 (from energy minimization)
Obtained ligand parameters from PRODRG.
Water model SPC, Usually the ligand atoms are blowing up after a couple of steps.
How do I increase the LINCS warn angle?
I will try the re-minimization again and see what results I obtain.
Using GROMOS53a6 force field.
Awaiting your response.
This kind of problem happened to me once. It might be that the topology of your ligand does not match the ligand description in your structure file. The atom order and name in each description has to be absolutely identical or they will be badly assigned, leading to the ligand explosion you describe.
Did you check for an exact match between your topology and your ligand structure?
Thanks Laurent for your reply.
I haven't checked the ligand structure and topology, thanks for pointing it out.
Your maximum force is insane! Before anythin else, double and triple check the outcome of PRODRG. I am am 99% sure the error is somewhere in your topology. Either PRODRG screwed up, or there is a description fail as Laurent said.
Good luck!
- Badges
- Science topic