Project

General

Profile

Bug #68

bilayer crashing in 3.3.1 but not 3.2.1 mdrun

Added by Alan Dodd over 13 years ago. Updated about 13 years ago.

Status:
Closed
Priority:
Normal
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

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.

topl3.2.1.tpr (6.53 MB) topl3.2.1.tpr 3.2.1 tpr file that only crashes with 3.3 onwards Alan Dodd, 04/10/2006 01:21 PM
ener files.tar (530 KB) ener files.tar tar of ener files from 3.2.1 and 3.3 Alan Dodd, 04/10/2006 01:32 PM
lipidtest.tar (360 KB) lipidtest.tar tar of tpr and run outputs from 3.3 and 3.2.1 Alan Dodd, 08/24/2006 12:44 PM
dopc.itp (13.6 KB) dopc.itp .itp for DOPC Alan Dodd, 08/31/2006 02:03 PM

History

#1 Updated by Alan Dodd over 13 years ago

Created an attachment (id=33)
3.2.1 tpr file that only crashes with 3.3 onwards

#2 Updated by Alan Dodd over 13 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 about 13 years ago

Hi,

I'm getting enormous LJ-14 energies with 3.3 when testing this:

Energies (kJ/mol)
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 about 13 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 about 13 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
two programs.

#6 Updated by David van der Spoel about 13 years ago

Just as a further comment, I have run thr tpr with mdrun 3.3 (CVS) and mdrun
3.2.1 and there are very minor differences in the forces as seen from the virial.

#7 Updated by David van der Spoel about 13 years ago

I have reproduced the problem with "standard" 3.3 code, will now retry with the
latest CVS.
Did you use surface tension coupling a lot before? THis might be the culprit
although I could not find any indication that something has changed.

#8 Updated by Alan Dodd about 13 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.

#9 Updated by Alan Dodd about 13 years ago

Created an attachment (id=67)
tar of tpr and run outputs from 3.3 and 3.2.1

#10 Updated by David van der Spoel about 13 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 about 13 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:
http://www.gromacs.org/gromacs/updates/changes-from-3.2-to-3.3-under-construction.html
Changed on 17 Jun 2005

In your input the Cut-off method is used for Van der Waals so this does not
affect you.

Could you post the top file (or itp file)?

#12 Updated by Alan Dodd about 13 years ago

Created an attachment (id=72)
.itp for DOPC

You're quite right, there is a typo in the 'pairs' section. Am testing with
this fixed.

#13 Updated by Alan Dodd about 13 years ago

With the typo fixed, the bilayers have not seperated for at least 750ps. The
extra stresses must have been causing the problem.

#14 Updated by David van der Spoel about 13 years ago

Since the system is stable I'm closing the bug.

Also available in: Atom PDF