bilayer crashing in 3.3.1 but not 3.2.1 mdrun
Versions of Gromacs after 3.2.1 crash due to an expanding box, due to
seperation of a lipid bilayer. Changing electrostatics does not stop the
effect, nor does increasing tau_p, suggesting another cause? Making the
leaflets more "sticky" by increasing rvdw slows the seperation down, but even
at 3nm the bilayers still seperate. The bilayer has been previously
equilibrated under 3.2.1 and run for a total of hundreds of ns without anything
similar occuring, suggesting something in the code has changed to cause this.
I'll attach an example 3.2.1 .tpr file that crashes, along with the .edr files
to compare, though I've tried several systems with almost identical results.
If there are any specific comparisons you wish me to run for you, let me know.
#2 Updated by Alan Dodd over 14 years ago
Created an attachment (id=34)
tar of ener files from 3.2.1 and 3.3
This example is using 3.3 (patched) rather than 3.3.1, though the same occurs,
and unfortunately my rather ancient cluster crashed on the 3.2.1 run. However,
140ps is long enough to see that the bilayers seperate.
#3 Updated by David van der Spoel almost 14 years ago
I'm getting enormous LJ-14 energies with 3.3 when testing this:
Angle Proper Dih. Ryckaert-Bell. Improper Dih. LJ-14
2.87184e+04 9.63811e+03 1.25768e+04 1.73143e+03 3.27403e+06
Coulomb-14 LJ (SR) Coulomb (SR) Coul. recip. Potential
2.42737e+04 1.68884e+03 -5.72327e+05 -2.08088e+05 2.57224e+06
Kinetic En. Total Energy Temperature Pressure (bar)
1.06528e+05 2.67877e+06 3.00081e+02 1.60075e+03
Is this what you get too? On the first step of the tpr?
#4 Updated by Alan Dodd almost 14 years ago
Yes. These values are comparable to other stable runs I have. In retrospect,
the forcefield I'm using (ffgmx modified for lipids) appears to be relatively
fragile and prone to producing errors (hence the odd 1.5fs timestep), although
I never experience a problem with 3.2.1 in normal use.
#5 Updated by David van der Spoel almost 14 years ago
Is your mdrun patched for the problem with PME order != 4 ? Otherwise, retry it
with pme-order =4 and ewald_tol = 1e-6
If you still have problems then please upload a new tpr file for a single
molecule that runs 20 steps, and writes the energies and forces at every step.
You are welcome to run gmxcheck on the output files from this tpr file with the
#8 Updated by Alan Dodd almost 14 years ago
I've always used surface-tension. pme_order was definitely patched, one of the
first things I checked.
Am attaching .tpr files for 20-step run as requested, along with the outputs
from the runs. Both runs were created using the 3.2.1 .tpr (topol.tpr), though
I have included the 3.3 version too. Am fairly sure this is using the patched
version, though I'm not sure of a simple way to check - I haven't used any
3.3.* versions for several months.
#10 Updated by David van der Spoel almost 14 years ago
There definitely is an error in your topology: you have specified a pair
interaction between atoms C36 and C38, which are too close together. This leads
to very large 1-4 interactions (try gmxcheck -f traj.trr) and hence instability.
This does not explain why gromacs 3.3 responds different than 3.2.1, but please
retest with this pair interaction taken out.
#11 Updated by David van der Spoel almost 14 years ago
The bahavior of 1-4 potentials has changed in 3.3, but only for shift and switch
potentials. This does include PME, but between the two atoms in your lipid there
is no Coulomb interation. This is documented here:
Changed on 17 Jun 2005
In your input the Cut-off method is used for Van der Waals so this does not
Could you post the top file (or itp file)?