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 _______________________________________________ gmx-users mailing list gmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
