[gmx-developers] Reaction Filed crash!

[David van der Spoel]

On 2011-10-03 20:38, Miguel Machuqueiro wrote:
>
> Hi again,
>
> Thank you Justin and Chris for your comments.
>
> We did follow that thread and we did prepare very complex initialization
> procedures using "nstpcouple=1" which increased the stability of our
> simulation. Nevertheless, in all my production runs the simulations
> crashed. We have a few tests running at the moment (already 20 ns) that
> are surviving for now. But you have to agree that something must be
> wrong when, using such a long/smooth initialization (and with
> nstpcouple=1), some systems are surviving and others are already dead.

So the production is with default nstpcouple (which is not 1)? Have you 
tried to set it to 1 again, or increase tau_P? 20 ps is not very long 
for long sims, in particular since the RF/Cut-off will lead to heating 
due to cut-off artifacts.

>
> I don't know if Itamar Kass managed to solve all his problems just by
> "doing a better initialization", but I can assure you that we can't.
>
> I am really curious that noone else has bumped into this one. Is there
> noone else using Reaction Field this days in the GMX community?
>
> I would appreciate any examples from people that have successfully
> simulate MD of a protein with RF.
>
> Regards
> Miguel
>
>
> On 03-10-2011 18:33, Chris Neale wrote:
>> You might be experiencing the same thing that was discussed on the
>> users list in September:
>>
>> http://lists.gromacs.org/pipermail/gmx-users/2011-September/064236.html
>> http://lists.gromacs.org/pipermail/gmx-users/2011-September/064237.html
>>
>> The conclusion there was to use a small nstpcouple (perhaps =1) during
>> equilibration.
>>
>> Chris.
>
> On 03-10-2011 18:34, Justin A. Lemkul wrote:
>> Miguel Machuqueiro wrote:
>>> Hi to all developers,
>>>
>>> We have been using GMX here in the group for almost a decade and we are
>>> familiar with almost all versions of the package.
>>> Recently, when testing the newer versions (4.5.x) using Reaction Field
>>> in a simple lysozyme system (4LZT), my system crashed constantly.
>>>
>>> System: 4LZT (129 a.a.) directly from PDB databank.
>>> Water: ~6250 molecules
>>> GMX: 4.5.x (all versions tested)
>>> ForceField: Gromos 53A6
>>> Long range elect.: RF, RF-nec, GRF or simply "Cutoff" (all with same
>>> fate.. crash!)
>>> Ensemble: NPT or NVT
>>>
>>> So, what do we already know about the problem:
>>> - The crash always starts in the NH bond present in the amides (main
>>> chain or side chains of GLN and ASN). In the stepxxxx.pdb, we observe
>>> the hydrogen jumping out of place.
>>> - The stability of the simulation is dependent on the quality of the
>>> minimization/initialization procedure. Simple/regular procedures leads
>>> to crashs within 10 ns of simulation. More elaborate protocols of
>>> minimization/initialization can delay the crashs a little longer.
>>> - The system crashes independently of the number of "threads" used
>>> (1-12).
>>> - The instability only leads to crashs when dealing with proteins>100
>>> a.a.. Simulating a 15 a.a. peptide for over 1 microsecond does not lead
>>> to a crash.
>>> - No problem/crash is observed when using PME.
>>> - No problem/crash is observed when using GMX versions prior to 4.0.7.
>>>
>>> My question is this, what might have happen with the RF code between
>>> versions 4.0.x and 4.5.x that leads to this problem? Or is there any
>>> specific parameter to be used with RF/Cutoff to avoid these problems?
>>>
>>> Even though RF is not very popular with GMX developers, it has been
>>> essential for the development of our in house Constant-pH MD method
>>> (http://www.gromacs.org/Documentation/How-tos/Constant_pH_Simulation
>>> <http://www.gromacs.org/Documentation/How-tos/Constant_pH_Simulation?highlight=constant+pH>).
>>>
>>>
>>> On a typical MD run we use the following MDP file, or a similar one.
>>>
>>> ________________
>>> integrator = md
>>> tinit = 0
>>> dt = 0.002
>>> nsteps = 5000000
>>> comm-mode = Linear
>>> nstcomm = 1
>>> nstxout = 0
>>> nstvout = 0
>>> nstfout = 0
>>> nstlog = 5000
>>> nstenergy = 5000
>>> nstxtcout = 5000
>>> xtc-precision = 1000
>>> energygrps = Protein SOL
>>> nstlist = 5
>>> ns_type = grid
>>> pbc = xyz
>>> periodic_molecules = no
>>> rlist = 0.9
>>> coulombtype = Reaction-Field
>>> rcoulomb-switch = 0
>>> rcoulomb = 1.4
>>> epsilon_r = 1
>>> epsilon_rf = 54
>>> vdwtype = Cut-off
>>> rvdw-switch = 0
>>> rvdw = 1.4
>>> DispCorr = No
>>> tcoupl = berendsen
>>> tc-grps = Protein SOL
>>> tau-t = 0.1 0.1
>>> ref-t = 298.15 298.15
>>> Pcoupl = berendsen
>>> Pcoupltype = Isotropic
>>> tau-p = 0.5
>>> compressibility = 4.5e-5
>>> ref-p = 1
>>> refcoord_scaling = No
>>> gen_vel = no
>>> constraints = all-bonds
>>> constraint_algorithm = lincs
>>> lincs-order = 4
>>> lincs-iter = 1
>>> lincs-warnangle = 90
>>> ___________________________
>>>
>>> Thank you for any help you can assist us.
>>>
>> Sporadic crashes have been reported for systems that are stable under
>> 4.0.x and
>> not 4.5.x. The solution that seems to have worked is to adjust the
>> nsttcouple
>> and nstpcouple values. See the discussion here (and subsequent posts
>> in the
>> thread):
>>
>> http://lists.gromacs.org/pipermail/gmx-users/2011-September/064232.html
>>
>> -Justin
>>
>>> Regards,
>>> Miguel
>>>
>>>
>>> --
>>> ============================================
>>> Miguel Machuqueiro
>>> Department of Chemistry and Biochemistry
>>> Faculty of Sciences, University of Lisbon
>>> Campo Grande, Edifício C8 (sala 8.5.47)
>>> 1749-016 Lisboa, Portugal
>>> Tel. : +351 217500112 (int.ext.28547)
>>> Mobile: +351 967562285
>>> E-mail: machuque at fc.ul.pt
>>> www1: http://webpages.fc.ul.pt/~mamachuqueiro
>>> www2: http://intheochem.fc.ul.pt
>>> ______________________________________________
>>>
>>>
>
>


-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se