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

[Berk Hess]

Hi,

There seem to be two different issues in this thread.
One is a problem of moving frozen atoms and the other is a discussion
on how important "exact" minimization is.

The first issue is a technical one.
In general one can have conflicting demands, one can freeze atoms
as well as require constraints. The question is what to do with frozen
atoms between which there are constraints that are not at the correct length
in the initial structure. I chose to use masses for frozen atoms that
are 10^24
higher then normal. In that way constraints between frozen atoms will
be correct, whereas frozen atoms should move 10^24 less than non-frozen
atoms connected by constraints. Therefore I would expect that your
water oxgyens should move very little.
By how much did they move?

The second question is on the importance of minimization.
For biological systems the relevant temperature is around 300 K.
Except for removing atomic clashes, or otherwise incorrect structures,
energy minimization will not help you to get more realistic structures,
since in reality you would have an esemble with a large spread.
This spread will in cost cases be much larger than the difference
between a non-minimized and a minimized structure.

Berk


Leontyev Igor wrote:
>>>> "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?
>
> _______________________________________________
> 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.