[gmx-developers] Reaction Filed crash!

[Berk Hess]

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