[gmx-users] re: need help with free energy of solvation calculation

[David Mobley]

Jonathan,

> First, a question: I know that certain values of sc-alpha are preferred for different cases, but are there any values that definitely should not be used?  In my calculations, some hydrogens (with charges but no L-J) are mutated into united atom methyls (with charges and L-J), so the dgdl curve has a large, negative slope near lambda=0.  I had started with sc-alpha=1.51 and then tried sc-alpha=0.51 and then sc-alpha=0.2 as those changes reduced the slope of dgdl.  Is there any reason that I shouldn't continue to reduce sc-alpha to produce a better-behaved curve?  By the way, the results with these different sc-alpha values are all in good agreement with each other, but the uncertainty tends to be smaller with smaller sc-alpha.

No, it is fine to fiddle with this. In my experience if you use
sc-power=1.0 (default in 3.3) and sc-alpha=0.5 it is roughly optimal;
if you make sc-alpha any smaller with sc-power=1, you get back other
undesirable features. But it sounds like your strategy is on the right
track. If you are using sc-power=2 (as in the original Beutler work),
I'd suggest sc-power=1, because it seems to work better in terms of a
smoother curve (see also the work of M. Shirts on this). However, if
you are getting results in good agreement with one another, it is
reassuring (I was getting this in my tests as well).


> I've now run cases for 4 different molecules, and I only find excellent agreement between my results and the published results for 1 of the 4.  The other three may or may not be overlapping depending on the large error bars, but certainly don't agree well.  I'd be interested in any comments.  I can continue to refine my method (e.g., changing L-J and charges separately, as has been suggested on the list), but I fear I may be stuck with my results in my actual study being (hopefully) internally consistent with each other but maybe not in very good agreement with published results.

> To answer David M.'s question here:
> http://www.gromacs.org/pipermail/gmx-users/2006-January/019526.html
> The Yu et al. paper uses the GROMOS code with soft cores and single-step perturbation.
>
> They ran for a total of 5 ns.  Since I'm using multiple lambda values instead of single-step perturbation, so far I've only run about 1 ns for each case.  As far as some of the other validation that I've done, I did confirm that for a given cellobiose conformation I calculated exactly the same energy components as did some of the original developers of the force field.

Haven't read the Yu paper, but I am suspicious of single step
perturbation and am inclined to believe your results first. If you
like you can try our methane in water test as a simple test case to
compare with another published result done using a method more similar
 to yours. http://www.dillgroup.ucsf.edu/~jchodera/group/wiki/index.php/Free_Energy:_Tutorial

> The differences between Run1, Run2, and Run3 are in sc-alpha and number and spacing of lambda points.  There are all free energies of solvation relative to cellononaose at 300 K in SPC water.  For Yu et al., the first number is the average from the full 5 ns simulation and the number in () is for the last 2 ns.
>
> Here is the comparison:
>
> 2,3,6-methylcellononaose
> =============================
> JDM-Run1        487 ± 130 kJ/mol
> JDM-Run2        463 ± 26 kJ/mol
> Yu et al.       393 (329) kJ/mol
>
> CE236Me
> =============================
> JDM-Run1        113 ± 39 kJ/mol
> JDM-Run2        121 ± 15 kJ/mol
> Yu et al.       224 (165) kJ/mol
>
> CE23Me
> =============================
> JDM-Run1        -3.4 ± 23 kJ/mol
> JDM-Run2        16 ± 16 kJ/mol
> JDM-Run3        19 ± 10 kJ/mol
> Yu et al.       136 (136) kJ/mol
>
> CE6Me
> =============================
> JDM-Run1        70 ± 15 kJ/mol
> JDM-Run2        79 ± 9 kJ/mol
> Yu et al.       68 (50) kJ/mol


If their error bars are the number in parentheses, then it looks to me
like you're within their error in every case, which looks like pretty
good agreement to me. In my opinion published error estimates are
often underestimates of the true error, so I suspect this is as good
as things will get for  you. That's probably the case even if these
aren't underestimates of the true error.

Again, I would reiterate that you can probably do things more
efficiently by turning off the charges seperately for the vdw. Right
now you're trying to find a set of soft core parameters that
simultaneously works well for both, which no one has been able to do
so far. It's easier to pick a set that works well for decharging
(simple linear scaling works great for that, and you can usually get
by with only four or five lambda values) and then a set that works
well for vdW (usually sc-power=1.0 with sc-alpha=0.5).

David