[gmx-developers] Reaction Filed crash!

Miguel Machuqueiro machuque at fc.ul.pt
Thu Nov 3 18:52:42 CET 2011 

Hi there,

Before I start, let me apologize for taking this long to reply. In fact, 
we wanted to test several of the suggestions before addressing this 
problem in the list.

I would also like to thank David van der Spoel and Berk Hess for your 
comments.

Davis, has Berk already mentioned, the crashes observed in the NVT 
simulations kind of rules out a barostat problem.

Berk, I think you really identified the source of the problem. But here 
is what else we know now:

  - MD with RF and a smaller protein does not crash, even if we increase 
the number of water molecules to achieve the same number of atoms as in 
the lysozyme system (which crashed promptly).
  - Better initiations (including nstpcouple=1) can delay the crashes 
several tens of nanoseconds but cannot eliminate them.
  - The use of an epsilon_rf=0 (infinity) does not influence the rate of 
crashes in our test simulations.

More importantly:

  - The lysozyme system with a twin range cut-off (0.8-1.4) together 
with regular RF settings and trivial initialization crashes within a 
couple of nanoseconds.
  - The same system with a single cut-off (1.4) runs for long 
simulations without crashing or even reporting a single LINCS warning.

It became clear to us that the leap-frog version of the reversible 
Trotter decomposition scheme is the origin of the problem. As of Gromacs 
version 4.5, it is not possible to do MD using the reaction field 
together with a *twin range* scheme. As you can imagine, this is a major 
limitation taking in consideration that the simulations with a single 
cut-off take 3-4 times longer then the twin range ones.

I wonder if there is any kind of bug affecting this that could be easily 
corrected in the following GMX versions. We would really appreciate it, 
specially because we rely heavily on the RF treatment of the long range 
electrostatics in Gromacs for our Constant-pH MD package.

Once again, thank you for your comments and your time.
Miguel


On 03-10-2011 20:50, Berk Hess wrote:
> On 2011-10-03 20:49, David van der Spoel wrote:
>> 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.
> NVT also crashes according to the information in the first mail.
>
> I changed the way the twin-range interaction scheme works for 4.5. Now a
> proper, reversible
> Trotter scheme is used, see manual section 3.4.7. The energy
> conservation is now an order
> of magnitude better. But there is still the (by default) unbuffered
> reaction-field. The lower
> energy drift of the twin-range interactions will lead to less kinetic
> energy being coupled off.
> But this in turn could possible lead to less cut-off heating being
> coupled off, which I would
> expect to have the largest effect exactly on the protons in the protein.
>
> I would always use an epsilon_rf=0 (infinity). This drastically
> decreases the energy drift,
> as the force is zero at the cut-off. According to me any non-zero rf is
> nonsense, as the cut-off
> artifacts are always worse than the "more correct" epsilon.
> If that doesn't help enough, you could try using reaction-field-zero
> with a buffer (which will be
> much slower).
> NB: for version 4.6 I have implemented new, efficient, proper cut-off
> reaction-field loops,
> although currently without twin-range cut-offs.
>
> Berk
>>> 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
>>>>> ______________________________________________
>>>>>
>>>>>
>>>
>>


-- 
============================================
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
______________________________________________

More information about the gromacs.org_gmx-developers mailing list