All, I need some help understanding the behavior of GROMACS 3.3. ( Apologies on the reference to the old version, but I published some data that now others are having trouble reproducing, which necessitates me using the version I originally ran the calculations in).
I am computing the hydration free energy of acetamide in TIP4P-Ew water using the FFAMBER ports, and I am having a problem in the component of the calculation where the A state has acetamide in water with no partial charges on it, and the B state has acetamide in water with all of its interactions with water turned off (but intramolecular interactions on) which is accomplished through some topology file trickery (modifications to the pairs section of the topology and explicit specification of intramolecular nonbonded interactions). Anyway, the expected value for this component calculation is something around -1.8 kcal/mol. My published result (in GROMACS 3.3) was -2.1 kcal/mol (+/- 0.05 or so). I tried to reproduce the result of -2.1 kcal/mol to confirm that the difference was legitimate, and found that I instead got -4.1 kcal/mol roughly. After some testing, I established that the difference was due to a single setting: In my original calculations, I had used PME, with rcoulomb = 0.9 nm and rlist = 1.0 nm. When I tried to reproduce these, I used PME, with rcoulomb=1.0 nm and rlist=1.0 nm. Can someone help me understand what exactly is going on here with this (seemingly small) change in the mdp files making such a large difference in the computed hydration free energies? I realize that in later (3.3.1 and after) versions of gromacs, the code insists that rcoulomb=rlist with pme, essentially (which still seems strange to me), but I still am not clear on why this is (and what exactly is going wrong in the calculations). For what it's worth, I still get the correct density of water with both sets of parameters. Yet this parameter change seems to cause a *systematic* difference of some sort -- for example, if I look at the plot of dV/dlambda values (http://www.dillgroup.ucsf.edu/~dmobley/uploads/novdw_comparison.pdf) the differences between the two cases are entirely within a specific range of lambda values, and are in a consistent direction. Yet I find no systematic differences in system volume or density at these lambda values. That's my question -- any help will be appreciated. The astute will notice I now also have a separate problem: If the "correct" input parameters to use are rcoulomb=rlist=1.0 nm, then the correct result is -4.1 kcal/mol, which is even further from the expected answer (-1.8 kcal/mol) than I started off with. This is the expected answer, in that it is the answer I arrive at using other simulation packages (DESMOND, for example) with the same input parameters, and the answer my collaborators (using several other packages) are getting with the same input parameters. Obviously this is a cause for some concern and suggests there may be some sort of an underlying problem here. On my list of things to do is to see whether I still get the same results in GROMACS 4. Thanks for your help, David Mobley _______________________________________________ gmx-users mailing list [email protected] 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 [email protected]. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
