[gmx-developers] Reaction Filed crash!

[Berk Hess]

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

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