[gmx-users] Free energy calculation problems (bug?)

[Michael Shirts]

Two thoughts that probably don't explain the problem, since they
should have similar behavior between versions:

* vdw-type                 = Cut-off

is almost never a good idea.  If you are doing this to be
reproducible, it's fine, but I'd avoid it in the future.

* You almost certainly want DispCorr = EnerPres, or else the density
will change with cutoff, and dispersion interactions beyond the cutoff
will be neglected (also changing the energy with cutoff).

Can you file a redmine (http://redmine.gromacs.org) with .top  and
.gro files as well?  One never knows where a problem might be coming
from (maybe only with certain top files?), so having those helps as
well.  It would be better to have the softcore version, since that is
likely to be better behaved anyway, and it has the same difference.


On Sat, Feb 8, 2014 at 12:00 AM, Олег Титов <titovoi at qsar.chem.msu.ru> wrote:
> Dear GROMACS users!
>
> I've encountered a problem in thermodynamic integration with GROMACS. I'm
> interested in calculating of solvation free energy differences between
> holbenzenes and benzene. To calculate this values I use thermodynamic
> integration approach (5 point formula).
>
> Half a year ago I did my calculation with GROMACS 4.6.1. The results looked
> good and were consistent with literature data. Recently, I've discovered
> that versions 4.6.0 and 4.6.1 had a bug in free energy code, so I decided
> to recalculate everything with 4.6.5 version of GROMACS. New results look
> not so good and, moreover, are not consistent with literature.
>
> I use GAFF forcefield with TIP3P water models. To obtain GROMACS topology
> and coordinate files I convert AMBER topologies (generated with tleap) with
> amb2gmx.pl script and this topologies look fine.
>
>  For phenylbromide molecule calculated solvation free difference should be
> near -0.1 - -0.2 kcal/mol (Experimental data is 0.59 kcal/mol).
> 4.6.1 version gave me -0.1+-0.13 kcal/mol - Fine!
> I've performed the same calculation with AMBER11  - -0.27 +-0.16 kcal/mol -
> also fine.
> 4.6.5 version for the same system gives me 1.1 - 1.2 kcal/mol
>
> I've played with some options (tried adding soft-core potentials and turned
> on/off couple-intramol option). The results are:
> no_sc no_intramol 1.27 +- 0.18
> no_sc intramol 1.20 +- 0.20
> sc no_intramol 1.35 +- 0.21
> sc intramol 1.22 +- 0.19
>
> So now I don't know what to do. Old buggy version produced good results
> consistent with literature and other MD packages and new fixed version
> produces different results. The dH/dl curve for two versions of GROMACS
> looks similar for lambdas 0-0.7 (mutation is PhBr -> PhH ) but becomes
> different at lambda = 0.9 So something changed there... Moreover different
> elecrostatic approaches for PhBr molecule are consistent with each other
> (dddG are fine), but new version of GROMACS gives me overestimated ddG as
> shown above.
>
> Can anyone give me some advice what to do with it? May be there is some bug
> in new version? May be I'm doing something wrong?
>
> Attaching sample .mdp (for solvated and vacuum system) and .top files.
>
> Thanks for the help,
> Oleg
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
>