[gmx-developers] Reaction Filed crash!

[baptista]

 On Sun, 18 Dec 2011 16:29:30 +0100, David van der Spoel wrote:
> On 12/18/11 1:27 AM, baptista at itqb.unl.pt wrote:
>> On Sat, 17 Dec 2011, Berk Hess wrote:
>>
>>> I will try to clarify the situation a bit.
>>> Pre-4.5 there were 4 approximation/errors in this typical "Gromos" 
>>> setup
>>> and a complicating factor:
>>> 1. Non-buffered pair-lists reused for nstlist-1 steps: no real
>>> Hamiltonian
>>> 2. Reaction-field with epsilon_rf<infinity: discontinuous energy 
>>> and
>>> force at cut-off
>>> 3. Twin-range interactions (or multiple time stepping) with 5*2=10 
>>> fs
>>> for the forces
>>> beyond 0.9 nm, but the reorientation of water only allows time 
>>> steps
>>> up to 5 fs
>>> 4. A non-reversible multiple time stepping scheme
>>> 5. A thermostat with strong coupling (pre-4.0 also Berendsen: one
>>> other error)
>>>
>>> Factors 1-4 cause serious integration errors in the long-range
>>> electrostatics
>>> which all lead to serious heating of the system, but different for
>>> water and protein.
>>> Factor 5 couples off this heat, but with such strong heating you no
>>> longer know
>>> what happening.
>>>
>>> In version 4.5 I fixed factor 4. This takes out one source of error
>>> and heating.
>>> But as different errors affect different aspects/components
>>> differently, it's very
>>> difficult to tell what will happen.
>>> Factors 1-3 can be fixed at a high computational cost (although 4.6
>>> might help here).
>>>
>>> My hypothesis is that factor 4 caused massive heating everywhere 
>>> and
>>> combined
>>> with a strongly coupled thermostat, all effects a seriously 
>>> dampened,
>>> including those of 1-3.
>>> Fixing factor 4 gives less heating, thus less thermostat coupling 
>>> and
>>> possibly less
>>> uniform heating between water and protein (due to effects of 1-3 
>>> being
>>> damped less),
>>> causing hydrogens in the protein to overheat. But as there are so 
>>> many
>>> complicating factors,
>>> it's difficult to find out what's really happening.
>>> When I fixed factors 1-3, by different mdp settings, I get perfect
>>> energy conservation in 4.5,
>>> but that's at a high computational cost.
>>>
>>> Cheers,
>>>
>>> Berk
>>
>> Thanks for clarifying the different factors involved. I see what you
>> mean. The relation between these factors is intricate and difficult 
>> to
>> antecipate, and trying to fix one of them ended up making things
>> unstable when using some setups (such as ours).
>
> As I pointed out earlier, due to large inherent heating in your
> setup, the system becomes unstable, and one can delay the "fatal
> fluctuation" but not indefinitely.
>
>>
>> You mention that 4.6 might help with factors 1-3. Does that mean 
>> that
>> you are making further changes? Assuming that we decide to stick 
>> with
>> our current setup, would you then advise us to stay with 4.0 for now
>> and wait for 4.6 to be released?
>
> If Berk has implemented these things correctly (which he usually
> does) 4.6 will make it more efficient to run with perfect energy
> conservation, but the inherent problems in your setup remain.
>
> What most likely "saved" you in previous versions, is the twin-range
> scheme (1) where the forces from far-away molecules are kept constant
> even though atoms may  have  moved out of the cut-off or into the
> cut-off.  In addition, the RF with epsilon < infinity in the original
> implementation has a minimum in the energy slightly beyond the
> cut-off. This may also improve stability.
>
> In summary: energy conservation is your friend!

 I totally agree! But realistic symmetry and boundary conditions (in its 
 broadest sense) are also your friends, and I tend to regard them as 
 crucial aspects in any physical model. So, I tend to follow them and 
 hope that the heat bath takes care of the problems associated with the 
 energy drift. But, as Michael rightly noted, even when that manages to 
 avoid numerical unstabilities and crashes, it risks putting your system 
 in a steady state. Choosing between friends it's not easy...

>
>>
>> Cheers,
>> Antonio
>>
>>>
>>> On 12/17/11 22:12 , Shirts, Michael (mrs5pt) wrote:
>>>> A couple thoughts:
>>>>
>>>> 1. if something was stable before (even if it was not best
>>>> practices), then
>>>> I agree that it should be stable in the current code, unless the
>>>> previous
>>>> algorithm was wrong to begin with. I haven't looked into the 
>>>> details
>>>> best
>>>> to know what the problem is yet.
>>>>
>>>> 2. Massive heating damped with temperature control will drive you
>>>> away from
>>>> an equilibrium to a steady state where some properties become
>>>> incorrect even
>>>> if others are correct. So you do have to be very careful that the
>>>> properties you looked at before are the same ones. Where exactly 
>>>> the
>>>> failure point is is very hard to say very clear information about 
>>>> the
>>>> specific system.
>>>>
>>>> 3.
>>>>> We can also discuss if, as a general rule, using a continuum 
>>>>> reaction
>>>>> field is better or worse than using a lattice method such as PME,
>>>>> but that
>>>>> is a different issue. Like many other people (Wilfred van 
>>>>> Gunsterem,
>>>>> Arieh
>>>>> Warshel, Alan Mark, Philippe Hünenberger, etc), I'm really not
>>>>> convinced
>>>>> that lattice methods offer any real advantage, because I find the
>>>>> published evidences for either their benefits or their artifacs 
>>>>> to be
>>>>> rather weak or contradictory.
>>>> For HOMOGENEOUS systems, then I agree; the configurations sampled 
>>>> can be
>>>> achieved by cheaper methods than PME -- basically, the 
>>>> configurations
>>>> are
>>>> dominated by short-range effects, and beyond a certain range, all
>>>> that the
>>>> extra interactions do is affect the overall energy, which can be
>>>> calculated
>>>> pretty well by many methods (including reaction field).
>>>>
>>>> When lattice methods are needed is inhomogeneous systems, such as
>>>> membranes
>>>> and interfaces, because the long range order affects the 
>>>> configurations
>>>> sampled, and RF and other continuum methods simply can't handle 
>>>> those
>>>> well;
>>>> properties will depend on cutoffs, etc.
>>>>
>>>> Best,
>>>> ~~~~~~~~~~~~
>>>> Michael Shirts
>>>> Assistant Professor
>>>> Department of Chemical Engineering
>>>> University of Virginia
>>>> michael.shirts at virginia.edu
>>>> (434)-243-1821
>>>>
>>>>
>>>>
>>>
>>> --
>>> gmx-developers mailing list
>>> gmx-developers at gromacs.org
>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>> Please don't post (un)subscribe requests to the list. Use the www
>>> interface or send it to gmx-developers-request at gromacs.org.
>>>
>>
>>
>
>
> --
> David.
> 
> ________________________________________________________________________
> David van der Spoel, PhD, Professor of Biology
> Dept. of Cell and Molecular Biology, Uppsala University.
> Husargatan 3, Box 596,  	75124 Uppsala, Sweden
> phone:	46 18 471 4205		fax: 46 18 511 755
> spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se
> 
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++