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

[Stéphane Téletchéa]

David van der Spoel a écrit :
> Stéphane Téletchéa wrote:
>> 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 don't think so. All-bonds with 2 fs should work reliably with both 
> lincs and shake, provided your system is equilibrated to some extent.

I'm ready to provide again detailed steps i've followed, but on 
different systems, including systems coming from a working gromacs run 
(for instance a snapshot after 10 ns of production) the problem arises 
when i switch from h-bonds to all-bonds.

There are many reports about segfaults in gromacs while using all-bonds 
constraints, this was assumed to be related to pme (and gromacs 3.3.1 
was released to correct this). I *suppose* this may hide another problem 
and pme is not alone the culprit.

I'll try to continue this on bugzilla with input files, procedure on the 
wiki (since the procedure in itself works with h-bonds) so we can 
decipher where this comes from..

Cheers,
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