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

[Moore, Jonathan (J)]

This is a follow-up from my post in January regarding my attempts to reproduce some free energy of solvation calculations for modified cellulose oligomers.

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.

See my previous post here for the reference for the results I'm trying to reproduce:
http://www.gromacs.org/pipermail/gmx-users/2006-January/019523.html

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.

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

Thanks,
Jonathan

____________________________
Jonathan Moore, Ph.D.
Research and Engineering Sciences - New Products
Core R&D
The Dow Chemical Company
1702 Building, Office 4E
Midland, MI 48674  USA
Phone:  (989) 636-9765 
Fax: (989) 636-4019
E Mail: jmoore2 at dow.com