[gmx-users] Crashes when doing free energy due to pressure and LINCS warnings
matteo.aldeghi at lincoln.ox.ac.uk
Thu Jan 16 14:18:44 CET 2014
I am trying to carry out some alchemical free energy calculations for the
binding of a small molecule. I run energy minimisation (10000 steps), NVT
(10 ps) and NPT (100 ps) equilibration steps before the NVT production run
(1 ns). The problem is that I get crashes at the beginning of the NPT, after
receiving a series of pressure scaling and LINCS warnings. I am trying to
modify the restraints only.
I am having troubles figuring out what I'm doing wrong, since the crashes do
not happen at low lambda values (e.g. lambda <= 0.1) or if I run the same
simulation with full restraints but without using the free energy code
(free-energy = no). I ran the simulations more times obtaining the same
results. I apply the distance restraint using the pull code, whereas I
specify 2 angle and 3 dihedral restraints in the topology (as shown below).
Also, I tried running simulations for the decoupling of the ligand charges
(defining a B state with zero charges) using:
fep-lambdas = 0.0 0.25 0.5 0.75 1.0
restraint-lambdas = 1.0 1.0 1.0 1.0 1.0
Now, I ran a simulation at lambda = 0.5 where I commented out all the
restraints from the topology and mdp files, and it crashed. On the other
hand, I ran the same simulation (restraints commented out) twice, but
filling the restraint-lambdas vector with 0.0, and they both ran fine until
I noticed that the crashed simulations have large negative pressures (-10^5
bar) and large positive potentials.
- gromacs 4.6.3 and run mdrun_mpi with particle decomposition.
- Amber99SB-ILDN and GAFF with AM1-BCC charges for the ligand.
Any suggestion would be much appreciated,
* I tried playing around with the tau_p value with no success.
* I ran a couple of simulations with 4.6.1 and they crashed too at the same
* the B-state constant for the pull code is not defined as it should, as I
mean to have k1=0 and kB1=4184 (rather than only k1=4184). This is a
mistake, however I believe this shouldn't affect the crashes.
* I couldn't find information about whether there should be a B-state
specified for the angle and dihedral restraints, am I right to assume the
lambda value scales the force constant directly?
Example of warnings
Step 17, time 0.034 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 349.859101, max 17941.974609 (between atoms 2297 and 2298)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
2084 2086 89.6 0.1330 0.1240 0.1335
2086 2087 89.8 0.1006 0.1440 0.1010
# many more bonds follow
Step 18 Warning: pressure scaling more than 1%, mu: 156.752 156.752 156.752
# a mix of these 2 warnings continues until the simulation crashes with a
LINCS fatal error due to too many warnings.
in the topology file
[ angle_restraints ]
; ai aj ak al type th0 fc multiplicity
1393 1391 2615 1391 1 88.8 4.1840e+01 1
1391 2615 2614 2615 1 32.9 4.1840e+01 1
[ dihedral_restraints ]
; ai aj ak al type phi dphi kfac
1410 1393 1391 2615 1 -159.7 0.0 4.1840e+01
1393 1391 2615 2614 1 122.6 0.0 4.1840e+01
1391 2615 2614 2610 1 12.8 0.0 4.1840e+01
View this message in context: http://gromacs.5086.x6.nabble.com/Crashes-when-doing-free-energy-due-to-pressure-and-LINCS-warnings-tp5013820.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
More information about the gromacs.org_gmx-users