[gmx-developers] Re:Re:Re:Re:Re: Energy Minimization

[Leontyev Igor]

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

>>  Leontyev Igor wrote
>> Anyhow, the terminology does not affect the problem: the
>> optimization of water-hydrogens with frozen positions of water-oxygens
>> results to the structure with shifted water-oxygens.

>Mark Abraham
> The terminology does affect the problem if your description of the
> existence of a problem is confounded by a frameshift artifact.

There is no frameshift. In my case all atoms are frozen (freezedim YYY)
except just few water HYDROGENS (TIP3P).  All atoms remained indeed on the
same positions after the minimization ("steep"), except only those few
"frozen" water OXYGENS which are bonded to the released HYDROGENS. Centering
the input and output structures by those shifted water OXYGENS makes
situation even worse.



>>> Tsjerk Wassenaar wrote:
>>> 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.
>>>

This is correct.
But the second L-BFGS run starting from the first L-BFGS endpoint makes no
EM steps ("Low-Memory BFGS Minimizer converged to machine precision in 0
steps.") The initial configuration is weird because hydrogen of OH- group is
residing almost on neighboring N-atom (charge<0). This configuration
corresponds to the REGION of singular energy minimum (-infinity). Just to
remind, the last has been obtained in the first L-BFGS run starting from the
reasonable configuration (endpoint of "steep" EM) and overcoming a high
barier.



> Chris
> Most of us are simulating things on the ns timescale that are
> experimentally probed on the ms timescale, agreed. The idea that em
> can possibly get you out of this problem when md can not seems
> ludicrous.You talk about em being able to "bring our system to
> another (lower energy) local minimum separated from the first one by a
> high barrier (>>kT)" although this defies even the most simple
> definitions of a "minimization".

EM deals with a potential energy profile while MD tends to minimize a free
energy. In general, these profiles are different. If EM's downhill pathway
goes in between two high barriers through some narrow valley (like a canyon)
of the potential energy, while the free energy profile does not have the
narrow exit at all, I don't see the reason why EM "principally" can not
result on better than MD minimum. Here I would prefer to question rather
than to make statements since I am not an expert in the minimization
algorithms.

> Chris
> I do, however, readily admit the possibility that I don't comprehend
> your main point because I have not seen your structure and attempted
> em. It seems likely that posting your structure, em attempt, output
> gro, etc. to some file-sharing utility and then directing us there
> would be useful.

tpr-file corresponding to "steep" minimization is already available via
bugzilla: (http://bugzilla.gromacs.org/show_bug.cgi?id=332). I am ready
to share all other i/o files for any MD/EM runs.
############################################################
                                     T H E    S Y S T E M
############################################################
     The system is an enzyme which pumps protons across a membrane and
against a proton gradient. The mechanism of the proton pumping  is quite
complex and currently not completely known. It is known; however, that the
energy necessary to move protons against the proton gradient is supplied by
chemical reactions in the catalytic site of the enzyme, while a dynamics and
pathways of the proton transfer is controlled by electron translocation
accompanied by opening/closing some proton gates (protonatable key
residues).
     To observe the proton gating process in MD we need to start from the
configuration which is sufficiently close to the point of the gate opening
or closing. Otherwise, 10-20ns MD run will not be enough to relax the system
to the point when the opening/closing become statistically favorable. A
possible way to deal with this problem is manual preparation of the start
configuration. Particularly, internal water molecules (which are typically
not seen in the crystal structure) forming proton entrance/exit channels,
are manually placed in appropriate positions. To avoid clashes, the prepared
structure then should be minimized keeping gates and water molecules in
desirable state. Thus, we come to minimization with constraints. The
minimized configuration and consequent MD depends dramatically on a way of
the constrained minimization, i.e. what is constrained/released and in what
order of the multistep procedure. Particularly, depending on what I
minimize first (water, gate or surrounding protein environment), the
manually added water might stay or be expelled from the water pocket next to
the gate. The water may even stay in the desirable position after
minimization but after 2 ns MD run it is gone again. Playing with the
constrained EM it is possible to find better (lower) minimum which results
to desirable dynamics when the water stays right next to the gate during
whole simulation. This is what I called "the system dynamics depends on the
minimization." Finer minimization algorithm with constraints provide more
flexibility to deal with such problems. However, as I figured out,  the
current implementation of the minimizers (even steepest descent) do not
always correctly treat constraints such as TIP3P water.



>Mark Abraham
> EM algorithms are unlikely to cross any energy barriers. They all work
> by looking at the value of the function and/or (approximations to)
> various derivatives, deducing a likely direction to move, following a
> heuristic to guess how far to move, and then trying again. In the limit
> of a continuous target function (true for MM), accurate energies and
> forces (true for MM), accurate second derivatives (often true for MM,
> but perhaps not for some EM algorithms), and a sufficiently conservative
> heuristic for choosing step sizes, then an EM algorithm cannot leave the
> local minimum. Thus there can exist a point when the system is in the
> region of a suitable minimum, and further EM won't affect the quality of
> the starting structure for MD on an MM force field. However if the
> assumptions don't hold, as they might for systems with frozen atoms,
> constraints, or an EM algorithm that starts with an approximation to the
> Hessian (matrix of second derivatives), or an EM stepsize that is too
> large, then EM might cross barriers.

Indeed, in the case of continuous function U(r) all EM algorithms by
construction goes downhill therefore they can not go above a barrier.
Although, a behavior of different algorithms in narrow valleys as well as in
the region of saddle or stationary point is not crystally clear for me,
such cases seem to be rather rare in real systems. In general, this
discussion is valuable for me and I would really like to thank all
participants.
    Nevertheless, the unconstrained minimization was not an initial issue of
the thread. All problems pointed out in my initial post are related to the
minimization with constrained/frozen atoms. According to what you are
saying, the constraints introduce a discontinuity into the "target function"
or singularity to Hessian matrix which makes numerical EM algorithms
unstable. Isn't it? But on a top of my mind I don't see a source of the
discontinuity if the constrained/frozen degrees of freedom would be excluded
from the phase space (i.e. considered as an external field) during the
search. Is the discontinuity fundamental problem or just an artifact of
implementation?



>>  Leontyev Igor wrote
>> "Second example" is an optimization of hydrogen positions with frozen
>> heavy atoms which typically is employed for continuum electrostatic pKa
>> calculations. The "Simple option" is a capability to carry out the "cg"
>> minimization with constraints.

> Mark Abraham
> I've forgotten the original context for this. The "second example" was
> outside the original context of the discussion, which I understood to be
> EM preparing for MD. The "simple option" is something the developers did
> not envisage being needed, for the reasons they describe in 3.10.2. If
> enough people want it, it might get added. The joy of free software is
> that if you want it, you can write it too!

I didn't expect that "the original context" of the thread will be shifted
toward the discussion that a finer minimization is not important for MD.
Even if the section 3.10.2 is true there are still some other applications,
e.g. the  "second example", where the finer minimization is an issue.
Anyhow, the reason why I started this topic is that to report some technical
problems with CONSTRAINED energy minimization. If problem exist (see at the
beginning) then it should be at least documented even in the free software,
isn't it?