[gmx-users] Free Energy Calculations
Justin A. Lemkul
jalemkul at vt.edu
Sun Feb 17 01:20:10 CET 2008
Quoting David Mobley <dmobley at gmail.com>:
> Which gromacs version are you using?
Version 3.3.1 (forgot to mention that before, sorry!)
>
> I vaguely remember, that at one point I may have found a bug relating
> to atoms that initially started out with no charges and then had
> charges that were turned on; I think it might have been a problem with
> PME or something similar. This would have been a bug with 3.3.1 or
> CVS, prior to the release of 3.3.2, I think; you can probably find it
> on Bugzilla. I would think it would have been fixed for 3.3.2 if
> that's what you're using, but it might be worth digging up the
> bugzilla to see for sure. Let me know if you have any trouble finding
> it.
I seem to remember a related discussion across the list not too long ago, but it
entered a level of theory that was (and still likely is) beyond my level of
comprehension. However, I have not been able to find the bugzilla.
> OK. Can you explain where your starting configurations come from?
> Ordinarily, one would expect that "forward" and "reverse" only are
> meaningful when you're (a) doing nonequilibrium simulations, or (b)
> beginning simulations at lambda values from the ending configurations
> at preceding lambda values. You're not doing (a); are you doing (b)?
> That is, does the simulation at lambda number i+1 begin with the
> ending coordinates from the simulation at lambda number i?
I am running independent simulations at each lambda value. Is that not correct?
I haven't been able to determine this type of protocol directly from the
literature, so I suppose I am doing neither (a) nor (b). Should I be doing
things as you've suggested under (b)?
>
> Also, it might be useful to know how many (and which) lambda values
> you are using for each transformation. Your equilibration protocol
> might be useful to know as well.
I am running a protocol very similar (if not identical to) the tutorial you've
posted - steepest descents minimization, 10 ps NVT, 100 ps NPT, and 5 ns of
production.
In removing the LJ/Coulombic interactions, my lambda values are: 0, 0.005, 0.01,
0.02, 0.025, 0.05, 0.075, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8,
0.85, 0.9, 0.95, and 1.
In restoring these parameters, my lambda values are: 0, 0.05, 0.1, 0.15, 0.2,
0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.925, 0.95, 0.975, 0.980,
0.990, 0.995, and 1.
The lambda values are not evenly spaced, because I attempted these calculations
once before and found a large peak at lambda=0 that dominates the curve when my
molecules are solvated. By using tiny lambda intervals between 0 and 0.1, I was
able to define a more reasonable shape for this curve. I think it is due to the
strong hydrogen bonding between my molecules (between 5 and 8 -OH groups) and
the solvent. There is a paper that describes such an interaction and the
lambda values used to resolve it (I believe it's one of Berk's papers, sorry if
I've confused it!)
>
> I am not seeing anything wrong with your topologies, at least at the
> level of detail you've sent.
>
>
> I think if I were you, I'd begin by:
> (a) Checking bugzilla
> (b) Looking at the energy components of dV/dlambda at different lambda
> values for your forward and reverse coulomb transformations. In
> particular, the dV/dlambda components at lambda=0.2 in the forward
> direction, say, should be equal (within uncertainty) to those at
> lambda=0.8 in the reverse direction. Looking at the contributions
> component-by-component (LJ, PME, etc) will help you nail down whether
> there's a specific energy term which is causing most of the
> difference. Figure out where the big differences you're seeing are
> coming from. The differences are so big, even short test simulations
> should be able to figure this out. This may ultimately point to either
> (i) a bug, or (ii) some kind of problem with your topology that we've
> missed.
I looked into a few of these parameters; maybe some of this output will shed
some light on things?
***g_energy output from lambda=0.2 (turning off charges)***
Statistics over 2500001 steps [ 0.0000 thru 5000.0000 ps ], 8 data sets
Energy Average RMSD Fluct. Drift Tot-Drift
-------------------------------------------------------------------------------
LJ-14 101.306 8.05907 8.05869 -5.42024e-05
-0.271012
Coulomb-14 -19.7668 9.42063 9.42059 1.84376e-05 0.0921881
LJ (SR) -14.3683 3.74098 3.74094 -1.12996e-05
-0.056498
LJ (LR) -0.113404 0.0204458 0.0204391 -3.62549e-07
-0.00181275
Coulomb (SR) -53.4647 7.64063 7.6406 1.54684e-05 0.0773422
Coul. recip. 32.9191 3.73187 3.71678 -0.000232277
-1.16139
dVpot/dlambda 50.3906 6.61614 6.60645 0.000247962 1.23981
dEkin/dlambda -2.32581e-05 0 0 0
0
***g_energy output from lambda=0.8 (turning on charges)***
Statistics over 2500001 steps [ 0.0000 thru 5000.0000 ps ], 8 data sets
Energy Average RMSD Fluct. Drift Tot-Drift
-------------------------------------------------------------------------------
LJ-14 101.387 8.06986 8.06959 -4.6474e-05 -0.23237
Coulomb-14 -18.9631 9.49656 9.49653 1.4605e-05 0.073025
LJ (SR) -14.1166 3.80216 3.80207 1.77125e-05 0.0885627
LJ (LR) -0.11146 0.0205708 0.0205702 -1.075e-07
-0.000537501
Coulomb (SR) -53.8835 7.808 7.80761 5.42018e-05 0.271009
Coul. recip. 10.655 3.77966 3.77966 2.09831e-06 0.0104916
dVpot/dlambda -77.7396 6.61417 6.61293 8.86327e-05 0.443164
dEkin/dlambda -2.33049e-05 0 0 0
0
It looks like the biggest discrepancies are arising in the Coul. recip. term.
Would this suggest that I'm falling victim to the PME bug you alluded to
before?
Analysis of the corresponding .edr files from turning on/off LJ parameters
showed all of these values to be in good agreement, like I thought they might
be from the decent output I got earlier.
Thank you very much for your detailed reply!
-Justin
>
> Good luck!
>
> David
> _______________________________________________
> gmx-users mailing list gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
========================================
Justin A. Lemkul
Graduate Research Assistant
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul at vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/
========================================
More information about the gromacs.org_gmx-users
mailing list