[gmx-users] problem with free energy calculations

[David Mobley]

Michael,

> I calculated the free energy of solvation of water
> (basically its chemical potential). The literature value for TIP3P
> water, which I used, is about -6.3 kcal/mol -
> I get something like +2.5 (!) ...

I suggest doing the charging component separately. I *think* this is
more efficient than doing the steps together -- plus there are
separate literature values for the charging component (I think it's
about 8.3 kcal/mol if I remember correctly, but I can dig up the
references for you if you want). This would give you a better idea
where you're going wrong.


> my calculations are fairly standard, I use PME with
> conservative parameters, soft core potentials,
> do 21 windows, each one lasting 500 pico seconds, then
> after discarding the first 100 ps of each dgdl*xvg file
> i use the averages and do a numerical integration,
> the works ...

When I do the charges separately, I usually use a total of 21 windows
also -- 5 for charging, 16 for LJ. I suspect you'll actually get much
smoother dV/dlambda curves for both if you do it this way rather than
trying to both at the same time.

> I see that the fluctuations of my dU/dlambda values are fairly
> large ... can it be that the soft core parameters i use (sc_alpha=0.7,
> sc_power=1, sc_sigma=0.3) are inappropriate if, as i do here,
> i turn on/off the VdW and Coulomb interactions simultanioulsy ??

For LJ transformations, the most efficient choice of parameters is the
same but with sc_alpha=0.5 (based on my testing and that of Michael
Shirts). 0.7 is actually a LOT worse than 0.5.

More below.

> integrator               = md

I agree with others that you should be using sd (Langevin). This is
the only thermostat currently in GROMACS that actually gives you the
correct distribution of kinetic energies without other known defects
(Nose-Hoover has known defects, for example).  You can check out my
settings here, if you like:
http://www.dillgroup.ucsf.edu/group/wiki/index.php?title=Free_Energy:_Tutorial&printable=yes

> coulombtype              = PME
> rcoulomb                 = 1.2
> vdwtype                  = shift
> rvdw                     = 1.2

I'd use switched vdw. And I think you can get by with a smaller
rcoulomb and rvdw if you tweak your PME parameters a bit.

> pme_order                = 4
> ewald_rtol               = 1e-05
> ewald_geometry           = 3d
> epsilon_surface          = 0
> optimize_fft             = no
> fourier_nx               = 64
> fourier_ny               = 64
> fourier_nz               = 64

What does your fourier spacing end up being? We're using shorter
coulomb cutoffs with pme_order=6 and:
fourierspacing           = 0.1
; EWALD/PME/PPPM parameters =
pme_order                = 6
ewald_rtol               = 1e-06

This is actually based on rather extensive testing to see both (a)
what will give us low enough errors in the potential energy (relative
to a really long cutoff) and (b) what is adequate for free energy
calculations.

I really suspect the biggest issue may be your soft core settings and
the fact that you're doing both transformations at the same time. Are
you sure in CHARMM you're using the same soft core functional form?
sc-power=1 here gives you Michael Shirts' modified functional form,
which is better if with sc=alpha=0.5, but sc-alpha=0.7 may actually
make it worse. Also, doing charges and LJ at the same time could give
charge-charge overlaps that might cause problems...

I agree the thermostat choice is not ideal and suggest you switch, but
I doubt this is the main source of problems.

David