**From:** Sándor Kovács (*skovacs_at_wustl.edu*)

**Date:** Tue Jul 17 2012 - 19:35:35 CDT

**Next message:**felmerino_at_uchile.cl: "Umbrella Sampling plus restraints"**Previous message:**JC Gumbart: "Re: Instabilities in Minimization"**In reply to:**JC Gumbart: "Re: Instabilities in Minimization"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

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
*

*>>>>
*

*>>>
*

*>>
*

*>
*

**Next message:**felmerino_at_uchile.cl: "Umbrella Sampling plus restraints"**Previous message:**JC Gumbart: "Re: Instabilities in Minimization"**In reply to:**JC Gumbart: "Re: Instabilities in Minimization"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

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