[gmx-users] Stuck --> unstuck ...

[Stéphane Téletchéa]

I've found out where the problem lie (at least one of ...).

Small reminder: i've tried to setup a system consisting of HIV protease 
+ one ligand on a water box. I used the excellent protocol from John E. 
Kerrigan (Drug-enzyme tutorial) in order to set up my simulation and 
encountered some problems with it (basically: Box was shifted at least 
10 times. Please see log-file.).

I add some constructive comments on ways to explore where the problem 
was (changing the protocol, including ions into the water to get a 
better coupling bath, etc.). Thanks to all of them, even if sometimes 
humour is hard to catch by email!

I'm not interested in the system per se (i've done other dynamics 
working fine with this protocol) but was rather interested in providing 
an updated way of doing a drug-enzyme tutorial with latest force fields 
(g45a3, g53a6) and using the latest prodrg (since it allow the output of 
  g43a1-type atoms, mostly transferable to the newer ones).
I'll put this on the wiki soon.

I've checked many parameters except one: all-bons vs hbonds.

I do not recall a 'rationale' reason for using the all-bonds parameter 
(instead of hbonds where it worked for other dynamics simulations) but 
the fact is that if i run all-bonds simulations (whatever the starting 
system) with NVT+NPT, i see a crash with the error message like this:

Step 1  Warning: pressure scaling more than 1%, mu: -1.59647e+20 
-1.59647e+20 -1.59647e+20
Correcting invalid box:
old box (3x3):
    old box[    0]={-9.82630e+20,  0.00000e+00, -0.00000e+00}
    old box[    1]={ 0.00000e+00, -8.90833e+20, -0.00000e+00}
    old box[    2]={ 0.00000e+00,  0.00000e+00, -1.12456e+21}
new box (3x3):
    new box[    0]={-9.82630e+20,  0.00000e+00, -0.00000e+00}
    new box[    1]={ 0.00000e+00, -8.90833e+20, -0.00000e+00}
    new box[    2]={ 0.00000e+00,  8.90833e+20, -1.12456e+21}

using the *same* protocol (same starting structures, etc) and only 
changing all-bonds to h-bonds leads to a 'normal' trajectory (no error, 
no warning, no crash).

I remember gromacs 3.3.0 had a problem where pme=4 (that lead to the 
3.3.1 rapid release). There are a lot of reports on this mailing-list 
with crash/segfaults when using all-bonds, and it seems more when 
all-bonds are used instead of h-bonds.

Of course it may come from not completely converged dynamics (my setup 
is very rapid, the 4 PR steps last only 100 ps), from not enough 
minimsed structures, etc.

Could it be also possible there is a problem while using all-bonds (even 
with a 2fs time step) instead of in gromacs, and pme_order was just 
hiding it ?

I've some pointers to look at on the system (equilibration, rmsd, 
temperature, better minimisation if possible, and since it'll be on the 
wiki, it will be easier to comment), but i think this was important to 
share already.

Thanks a lot in advance for comments,

Stéphane

-- 
Stéphane Téletchéa, PhD.                  http://www.steletch.org
Unité Mathématique Informatique et Génome http://migale.jouy.inra.fr/mig
INRA, Domaine de Vilvert                  Tél : (33) 134 652 891
78352 Jouy-en-Josas cedex, France         Fax : (33) 134 652 901