From: JC Gumbart (gumbart_at_ks.uiuc.edu)
Date: Tue Jul 17 2012 - 19:17:05 CDT
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