[gmx-developers] Re: Energy Minimization

[Mark Abraham]

Igor Leontyev wrote:
>> Igor Leontyev wrote:
>>>>> The bug with minimization seems to be fixed now, see
>>>>> http://bugzilla.gromacs.org/show_bug.cgi?id=332. Thank you Berk. With
>>>>> corrected sources "csettle.c" and "minimize.c" the new 4.0.5 
>>>>> version is
>>>>> able
>>>>> to minimize my MD system at least on the same level as it was
>>>>> possible in
>>>>> 3.3.1 version.
>>>>>
>>>>> At this point I would like to mention 2 more limitation of the
>>>>> minimization
>>>>> procedure implemented in gromacs.
>>>>> 1) Since a dynamics of complex biological system depends on quality of
>>>>> the
>>>>> minimized structure, a fine minimization is very important.
>>>
>>>> Generally, I would think not. The process of assigning atomic 
>>>> velocities
>>>> randomly sampled from a suitable distribution "undoes" the finest part
>>>> of
>>>> the minimization. Then one needs to equilibrate, possibly under 
>>>> position
>>>> restraints, which also leaves the structure perturbed from the input
>>>> structure, but hopefully still in the same local minimum region. See
>>>> also
>>>> comments in manual section 3.10.
>>>
>>> It's true if the finer energy minimization do not exceed ~(kT)
>>> relatively to
>>> the "steep" method for each site of the complex protein+solvent system.
>>> But
>>> it's hard to ensure for each site. Moreover, there are applications of
>>> the
>>> minimized structures other than MD simulations, in which there is no
>>> atomic
>>> velocities.  For example, continuum electrostatic pKa calculation needs
>>> as
>>> an input the protein structure with optimized hydrogen positions.
>>
>> This example was excluded by your earlier description of "dynamics of a
>> complex biological system", but anyway for such a pKa calculation,
>> optimizing such protein hydrogens with respect to some (arbitrary?)
>> frozen configuration of waters doesn't seem likely to be better than
>> optimizing with respect to some flexible configuration of waters, or
>> even without waters. Averaging over several such flexible-waters optima
>> might be better still.
> 
> There is no rigorous prove in the manual section 3.10 that the finer
> energy minimization is not essential for MD of complex biological system.

There doesn't need to be. The point of doing MD is to see a system cross 
energy barriers (or, rarely, to see that it doesn't). Thus it is 
expected to move from the starting position, and at a certain point of 
EM, the structure is good enough, i.e. close to one in the ensemble of 
interest so that equilibration happens reasonably quickly and without 
failures of the integrator. There's also no proof that it is necessary, 
nor even an example ;-)

> Moreover, the second example also indicates that gromacs users meets
> unexpected difficulties even with such a basic feature of any MD code as a
> fine optimization of hydrogen positions. I even may agree, that rigid TIP3P
> water model physically is not better then flexible one but the use of TIP3P
> is conventional  (therefore does not require detailed explanation) and 
> there
> is no good reason to use something else unless it can be proven that it is
> essential. In other words, as a user, I would rather prefer to have the
> simple option implemented then to have the extensive discussion about why
> this option in some cases is not necessary.

Sorry, I'm not sure what the "second example" or "simple option" to 
which you refer are.

>>>>> Is there way in
>>>>> gmx to reFINE minimization of the system with many different
>>>>> constraints
>>>>> ?
>>>
>>>> See manual. You can make the waters flexible, or use L-BFGS for
>>>> single-processor minimization.
>>>
>>>
>>> I have been tried to use all described in manual options but result was
>>> not
>>> satisfactory.
>>> Flexible water model does not solve the problem for my system with
>>> constrained hbonds and other constraints or frozen atoms. L-bfgs results
>>> to
>>> the structure where some hydrogens are locked in the local minimum which
>>> is
>>> far away from the global one.
>>
>> Well maybe you need simulated annealing.
> 
> It's not the case for simulated annealing. It looks as a drawback of
> L-BFGS procedure because steep minimization does not result to the
> problem.

EM with an approximate Hessian with frozen atoms in an MM force field 
need not be a well-formed problem.

>>> I just ran l-bfgs minimization starting from
>>> the minimized by the "steep" method configuration and found that one of
>>> van-der-Waals-less hydrogen atoms in the final structure is transferred
>>> to
>>> the vicinity (~0.1A) of some negative nitrogen but far away (~3.0A) from
>>> its heavy atom . This local minimum is far away from global one due to
>>> huge
>>> bond-stretching penalty of the h-bond.
>>
>> It moved because it was relaxing some other more severe strain, created
>> either by the structure or the unbalanced model physics created by the
>> constraints.
> 
> This seems to be in contardiction with the statement that "steep"
> minimization produces good enough structure for MD.

There's no contradiction at all. It was probably in a position suitable 
for MD after steepest descent!

You did steepest-descent EM on a system with lots of frozen atoms, if I 
understand correctly. You then switched to L-BFGS with frozen atoms, 
which starts with an estimated Hessian. The Hessian updates then are 
based on calculated forces whose components are unbalanced because of 
the constraints. Such a procedure could well be numerically unstable. 
You may be proving that. Try a second L-BFGS from your first L-BFGS 
endpoint and see where it goes. If it's some different place, then you 
may be demonstrating that this minimizer can't deal with the 
(artificial) problem you have constructed.

>> What kind of non-water hydrogens should be modelled without
>>  vdW?
>>
> 
> The atom type HO (hydroxyl group) in AMBERFF  has no vdW.

OK, so if the L-BFGS problem is well-formed numerically, you may have a 
model physics that doesn't treat your problem well.

>>>>> 2) After steepest descent minimization  of the protein+solvent system,
>>>>> the
>>>>> positions of the frozen TIP3P water molecules are slightly shifted
>>>>> (~0.3
>>>>> A).
>>>>> The shift of frozen molecules is really bothering.
>>>
>>>> Are you centering the structures on the same frozen atom? Otherwise you
>>>> may just be seeing a shift of reference frame.
>>>
>>> "Frozen" (by freezegrps) means in respect to the initial configuration,
>>> isn't it?
>>
>> Dunno. Presumably so.
> 
> It was not a question, it was an answer regarding "shift of reference
> frame", which applicability in respect to frozen atoms is not clear for me.

"Frozen" might have meant with respect to the initial configuration, or 
as someone pointed out, that the atomic positions were merely never 
updated. So (for example) if the initial frame was translated for some 
reason, it could still be translated in the output. This could produce a 
deceptive frameshift. You can test for this by doing what I said - 
centering the input and output structures on the same frozen atom 
post-mdrun.

Mark