Re: Instabilities in Minimization

From: Sándor Kovács (skovacs_at_wustl.edu)
Date: Tue Jul 17 2012 - 19:35:35 CDT

Indeed we have. The trouble is, several of the sites identified by
propKA have very similar pKas but wildly different locations (and
hence local environments), suggesting the need for a more robust
method of differentiating between them.

Can you expand on what you mean by "dynamical measures of stability"?
Are you just referring to a convergence in RMSD for the entire
complex? If not, are there good references available?

I've thought about FEP--the charge balance could be maintained by
having a positively charged solvent ion disappear--but mutating each
possible sidechain into its protonated variant to obtain an estimate
of deltaG seems like overkill, no?
Also, would the deltaG calculated by FEP on a single proton even
appear significant when the entire complex (~20,000 atoms, before
solvation) is considered? I'm concerned about "signal-to-noise" ratio...

On Jul 17, 2012, at 7:17 PM, JC Gumbart wrote:

> I feel like dynamical measures of stability followed by free
> energies (FEP maybe?) would be much more informative. Have you also
> looked into using pKa calculations (e.g., propKa) on the complex as
> a starting point?
>
>
> On Jul 17, 2012, at 6:36 PM, Sándor Kovács wrote:
>
>> Thanks Aron and JC for the quick responses!
>>
>> Our collaborators have experimental measurements of the charge on a
>> protein complex, and we believe that charge plays a role in how the
>> complex assembles, but no one knows exactly where that charge is
>> localized (i.e. which sidechains are protonated). Our goal is to
>> predict (assign) the most likely protonation sites for such a
>> complex. Would a modified equilibration + minimization scheme as
>> Aron suggests below be sufficient?
>>
>> Cheers,
>> Sándor
>>
>>
>> On Jul 17, 2012, at 5:50 PM, JC Gumbart wrote:
>>
>>> At the very least, roundoff errors could be causing the deviations.
>>>
>>> I agree with Aron, and will further in saying that I would never
>>> use energies from a minimization as a metric for such complex
>>> systems. Can you explain more about what you really are trying to
>>> do?
>>>
>>>
>>> On Jul 17, 2012, at 5:34 PM, Aron Broom wrote:
>>>
>>>> I can't answer the part about the minimizer, but as a quick note
>>>> concerning your mentioning looking at differently protonated
>>>> sites: if you are thinking to use the POTENTIAL energy to decide
>>>> which protonation configuration is best, I think that will not
>>>> work. You should use the ELECT instead, as the full potential
>>>> will be taking into account bond parameters (dihedrals and such)
>>>> which may not be relevant to the question of which site is most
>>>> favourable for protonation, but which will change depending on
>>>> which sidechain has the proton attached. By using just the
>>>> electrostatic potential, you will overcome this problem and have
>>>> a better sense of where it is most favourable to place a given
>>>> charge, rather than where it is most favourable to make a fake
>>>> (since this is classical MD) covalent bond.
>>>>
>>>> Having said the above, I think minimization (at least this kind)
>>>> is a bad choice for comparing energies, because you ignore the
>>>> whole ensemble of states that are likely very close to the local
>>>> or global minimum, and which can be dramatically different
>>>> between two charge configurations. I would suggest
>>>> equilibrating, and then, once the system has stabilized (by RMSD
>>>> for instance), take a number of snapshots along another window of
>>>> equilibration and then compare them (minimizing each snapshot if
>>>> you wish).
>>>>
>>>> Just my two cents.
>>>>
>>>> ~Aron
>>>>
>>>> On Tue, Jul 17, 2012 at 6:14 PM, Sándor Kovács
>>>> <skovacs_at_wustl.edu> wrote:
>>>> Greetings experts,
>>>>
>>>> I had a question about the minimization algorithm as implemented
>>>> in NAMD.
>>>>
>>>> For a system of ~20,000 atoms, when starting from the same
>>>> initial conditions (READ: same atomic positions), the first 5000
>>>> steps of my minimizations are IDENTICAL (as tested with diff).
>>>> After this point, however, repeat trials begin to diverge from
>>>> one-another, and result in POTENTIAL energies differing by as
>>>> many as ±10 kcal/mol after 50K minimization steps.
>>>> These differences persist when the minimization is carried out
>>>> for 200K steps (I wanted to make sure each had completely
>>>> converged).
>>>>
>>>> e.g. for three separate trials after 5000 minimization steps, I
>>>> obtain identical results:
>>>>
>>>> [TRIAL 1]
>>>> TIMING: 5000 CPU: 184.51, 0.0350367/step Wall: 184.51,
>>>> 0.0350367/step, 1.89782 hours remaining, 106.558594 MB of memory
>>>> in use.
>>>> ETITLE: TS BOND ANGLE
>>>> DIHED IMPRP ELECT VDW
>>>> BOUNDARY MISC KINETIC
>>>> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
>>>>
>>>> ENERGY: 5000 1124.8303 4717.7323
>>>> 3713.3759 139.7911 -34583.7700
>>>> -6281.9927 0.0000 0.0000 0.0000
>>>> -31170.0331 0.0000 -31170.0331 -31170.0331
>>>> 0.0000
>>>>
>>>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 5000
>>>> LINE MINIMIZER BRACKET: DX 0.000499702 0.000999403 DU -0.0921908
>>>> 0.16983 DUDX -302.606 -66.3689 406.269
>>>>
>>>>
>>>> [TRIAL 2]
>>>> TIMING: 5000 CPU: 184.356, 0.0350742/step Wall: 184.356,
>>>> 0.0350742/step, 1.89985 hours remaining, 106.617188 MB of memory
>>>> in use.
>>>> ETITLE: TS BOND ANGLE
>>>> DIHED IMPRP ELECT VDW
>>>> BOUNDARY MISC KINETIC
>>>> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
>>>>
>>>> ENERGY: 5000 1124.8303 4717.7323
>>>> 3713.3759 139.7911 -34583.7701
>>>> -6281.9927 0.0000 0.0000 0.0000
>>>> -31170.0331 0.0000 -31170.0331 -31170.0331
>>>> 0.0000
>>>>
>>>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 5000
>>>>
>>>>
>>>> [TRIAL 3]
>>>> TIMING: 5000 CPU: 184.478, 0.035129/step Wall: 184.478,
>>>> 0.035129/step, 1.90282 hours remaining, 106.367188 MB of memory
>>>> in use.
>>>> ETITLE: TS BOND ANGLE
>>>> DIHED IMPRP ELECT VDW
>>>> BOUNDARY MISC KINETIC
>>>> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
>>>>
>>>> ENERGY: 5000 1124.8303 4717.7324
>>>> 3713.3759 139.7911 -34583.7702
>>>> -6281.9926 0.0000 0.0000 0.0000
>>>> -31170.0331 0.0000 -31170.0331 -31170.0331
>>>> 0.0000
>>>>
>>>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 5000
>>>> LINE MINIMIZER BRACKET: DX 0.000499698 0.000999396 DU -0.0921886
>>>> 0.169841 DUDX -302.606 -66.363 406.287
>>>>
>>>>
>>>> after 50K minimization steps, one sees a divergence in the tens
>>>> place:
>>>>
>>>> [TRIAL 1]
>>>> LINE MINIMIZER BRACKET: DX 0.000511803 9.74047e-05 DU
>>>> -0.000141415 5.12214e-06 DUDX -0.552612 3.6656e-08 0.105172
>>>> LINE MINIMIZER REDUCING GRADIENT FROM 1.15801 TO 0.00115801
>>>> TIMING: 50000 CPU: 1784.96, 0.0353003/step Wall: 1784.96,
>>>> 0.0353003/step, 1.47085 hours remaining, 106.535156 MB of memory
>>>> in use.
>>>> ETITLE: TS BOND ANGLE
>>>> DIHED IMPRP ELECT VDW
>>>> BOUNDARY MISC KINETIC
>>>> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
>>>>
>>>> ENERGY: 50000 1129.8369 4769.8503
>>>> 3772.9891 142.4910 -35065.8099
>>>> -6378.5568 0.0000 0.0000 0.0000
>>>> -31629.1995 0.0000 -31629.1995 -31629.1995
>>>> 0.0000
>>>>
>>>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 50000
>>>>
>>>> [TRIAL 2]
>>>> LINE MINIMIZER BRACKET: DX 0.000580501 0.001161 DU -0.000380669
>>>> 0.0010312 DUDX -1.1704 -0.141114 1.91753
>>>> TIMING: 50000 CPU: 1773.21, 0.0352761/step Wall: 1773.21,
>>>> 0.0352761/step, 1.46984 hours remaining, 106.558594 MB of memory
>>>> in use.
>>>> ETITLE: TS BOND ANGLE
>>>> DIHED IMPRP ELECT VDW
>>>> BOUNDARY MISC KINETIC
>>>> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
>>>>
>>>> ENERGY: 50000 1129.6438 4771.2127
>>>> 3764.2485 140.6916 -35041.1554
>>>> -6374.3260 0.0000 0.0000 0.0000
>>>> -31609.6847 0.0000 -31609.6847 -31609.6847
>>>> 0.0000
>>>>
>>>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 50000
>>>>
>>>> [TRIAL 3]
>>>> LINE MINIMIZER BRACKET: DX 0.000196034 0.000856693 DU -0.00200314
>>>> 0.0381771 DUDX -20.3922 -0.0399876 89.2622
>>>> LINE MINIMIZER REDUCING GRADIENT FROM 123.878 TO 0.123878
>>>> TIMING: 50000 CPU: 1836.9, 0.0352672/step Wall: 1836.9,
>>>> 0.0352672/step, 1.46947 hours remaining, 107.398438 MB of memory
>>>> in use.
>>>> ETITLE: TS BOND ANGLE
>>>> DIHED IMPRP ELECT VDW
>>>> BOUNDARY MISC KINETIC
>>>> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
>>>>
>>>> ENERGY: 50000 1128.4327 4767.2091
>>>> 3765.5526 139.9726 -35046.9333
>>>> -6368.0236 0.0000 0.0000 0.0000
>>>> -31613.7900 0.0000 -31613.7900 -31613.7900
>>>> 0.0000
>>>>
>>>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 50000
>>>>
>>>>
>>>> and after 200K minimization steps (to ensure convergence to a
>>>> local minima), the differences remain:
>>>>
>>>> [TRIAL 1]
>>>> LINE MINIMIZER BRACKET: DX 5.58613e-05 0.000330726 DU
>>>> -2.54659e-11 4.00178e-11 DUDX 2.62465e-08 3.68785e-08 9.98249e-08
>>>> TIMING: 200000 CPU: 7113.24, 0.0353183/step Wall: 7113.24,
>>>> 0.0353183/step, 0 hours remaining, 107.148438 MB of memory in use.
>>>> ETITLE: TS BOND ANGLE
>>>> DIHED IMPRP ELECT VDW
>>>> BOUNDARY MISC KINETIC
>>>> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
>>>>
>>>> ENERGY: 200000 1129.1383 4765.1845
>>>> 3784.9412 142.4877 -35093.0615
>>>> -6379.9375 0.0000 0.0000 0.0000
>>>> -31651.2473 0.0000 -31651.2473 -31651.2473
>>>> 0.0000
>>>>
>>>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 200000
>>>>
>>>> [TRIAL 2]
>>>> LINE MINIMIZER BRACKET: DX 0.000907904 0.000907904 DU
>>>> -7.80979e-07 5.97607e-07 DUDX -0.00161935 -0.000100987 0.00141738
>>>> TIMING: 200000 CPU: 7119.27, 0.0355447/step Wall: 7119.27,
>>>> 0.0355447/step, 0 hours remaining, 107.281250 MB of memory in use.
>>>> ETITLE: TS BOND ANGLE
>>>> DIHED IMPRP ELECT VDW
>>>> BOUNDARY MISC KINETIC
>>>> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
>>>>
>>>> ENERGY: 200000 1130.0652 4766.4965
>>>> 3776.1033 141.0543 -35082.0893
>>>> -6370.2935 0.0000 0.0000 0.0000
>>>> -31638.6636 0.0000 -31638.6636 -31638.6636
>>>> 0.0000
>>>>
>>>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 200000
>>>>
>>>> [TRIAL 3]
>>>> TIMING: 200000 CPU: 7281.36, 0.0350886/step Wall: 7281.36,
>>>> 0.0350886/step, 0 hours remaining, 107.398438 MB of memory in use.
>>>> ETITLE: TS BOND ANGLE
>>>> DIHED IMPRP ELECT VDW
>>>> BOUNDARY MISC KINETIC
>>>> TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
>>>>
>>>> ENERGY: 200000 1129.6536 4770.1929
>>>> 3775.7018 141.1868 -35087.2720
>>>> -6374.7239 0.0000 0.0000 0.0000
>>>> -31645.2607 0.0000 -31645.2607 -31645.2607
>>>> 0.0000
>>>>
>>>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 200000
>>>>
>>>>
>>>> I thought the minimization procedure was deterministic--i.e.
>>>> repeat trials should all converge to the same energy minimum
>>>> (ideally the global one, but I understand that they can be
>>>> 'stuck' in local minima).
>>>> Is there something going on (stochastic, perhaps) during the
>>>> minimization procedure which would account for this discrepancy?
>>>>
>>>> Alternatively, is there a better procedure I should be following
>>>> If I would like to obtain an estimate of the system's potential
>>>> energy for comparison to similar systems with, say, slightly
>>>> differently protonated sites?
>>>>
>>>> I can make input files available off-list if someone needs to
>>>> reproduce the above.
>>>>
>>>> Thanks,
>>>> Sándor
>>>> ________________________________________________________
>>>> Sándor Ádám Kovács, B.S. M.S.
>>>> President, Association of Graduate Engineering Students (AGES)
>>>> Ph.D. Candidate and Graduate Research Assistant
>>>> Department of Energy, Environmental and Chemical Engineering
>>>> Washington University in St. Louis -- St. Louis, MO 63130 -- U.S.A.
>>>>
>>>>
>>>>
>>>> --
>>>> Aron Broom M.Sc
>>>> PhD Student
>>>> Department of Chemistry
>>>> University of Waterloo
>>>>
>>>
>>
>

This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:22:16 CST