[gmx-developers] Reaction Filed crash!

Miguel Machuqueiro machuque at fc.ul.pt
Mon Dec 12 21:27:21 CET 2011 

Hi Berk, and all the devel guys.

Once again, thank you for not giving up on us.

It took me a while to do the tests to check your suggestions and my 
results are not very famous.

A decrease in the nstlist from 5 to 4 was not able to provide stable MD 
simulations in my HEWL system. All attempts (several by now) ended up in 
lincs crashed with several stepxxxxxx.pdb files being generated.
I also tried simulations using nstlist values of "2" and "1" to be sure:
  - The system with a value of "2" held on for a while, but eventually 
crashed after ~15 ns. I reproduced this crash in 3 replicates.
  - The system with a nstlist of "1" went through the 50 ns test without 
any problems. Unfortunately, this simulation was very slow; 4 times 
slower than PME (remember that nstlist=5 gives RF simulations much 
faster than PME).

The idea that RF in previous versions of Gromacs was stable probably due 
to error cancellation raises in me several concerns. Unfortunately, with 
the correction in the current implementations, I am not able to use the 
RF option for the treatment of the long range electrostatics.

If you are interested, I could send you the files used to run my tests.
Even though we modified our code in order to use GMX version 4.0.7 (in 
use at the moment), we are still very interested in using the latest and 
faster versions of Gromacs.

Thank you and all the guys doing the hard work!
Miguel



On 27-11-2011 14:59, Berk Hess wrote:
> Hi,
>
> I found a solution for your problem!
>
> The basic issue is that a time step for the long range forces of 5*2 fs
> = 10 fs is too long
> in a system where charges move fast (probably mainly the water hydrogens).
> With the old, non-reversible, scheme there must be  some cancellation of
> errors, which makes
> the integration stable (although still very inaccurate).
> With the new, reversible, integration scheme the hydrogens in the
> protein heat up by a factor 1.5
> with a T-coupling strength of  tau_t=0.1 (worse than the old scheme) and
> they become unstable.
> Changing the long time step to 8 fs (nstlist=4 with dt=2fs) reduces the
> energy drift by a factor of 4
> and makes protein simulations stable.
> With this 8 fs the energy drift, for my protein test system, with the
> new scheme is a factor 3 lower
> than with the old scheme.
>
> The proper solution is using PME. But you said you couldn't use that for
> whatever reason.
> I hope that that is not because you need pair forces. Note that the
> reaction field force is NOT a pair
> force, it is a collective effect of all charges in the cut-off sphere,
> very similar to PME.
>
> Cheers,
>
> Berk
>
> On 2011-11-03 18:52, Miguel Machuqueiro wrote:
>> 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