[gmx-users] free energy calculations

[Alessandra Villa]

> Dear Gromacs users/developers,
> 
> I want to calculate the free energy difference between the native 
> and the mutated form of a protein segment. Before starting these 
> calculations, I tried to simulate a simple change: isoleucine --> 
> valine, in water. I use Gromacs 3.0.5 and the forcefield G43a1.
> The change is:
> HOOC(NH2)CH-CH(CH3)-CH2-CH3 (state A) -->
> --> HOOC(NH2)CH-CH(CH3)2-Dummy (state B)  
> (The dummy atom has m=15.035, charge=0, LJ parameters C6=0, 
> C12=0).
> For the state B I kept the same parameters for the bonded 
> interactions as for the state A (bonds, angles, dihedrals).
> The MD parameters:
> dt=0.002 (ps), nsteps= 100000,
> constraints=all-bonds (lincs); T coupling=berendsen, T=300K;
> P coupling=berendsen, isotropic.
> I made simulations at different lambda: 0.1, 0.2 etc (delta-
> lambda=0), with or without soft-core interactios.
> The results: the dG/dlambda values vs. time have a bad 
> (asymetrical) distribution. There are many peaks towards negative 
> values; for example, if the mean value is 7.5 (kJ mol-1 nm-2 lambda-
> 1) with a st. dev. +/- 4.5, there are many peaks that reach -15, -18.

It is common to have high fluactuation in the values of 
dG/dlambda. 

If you are slowing growing from lambda 0.0 to 1.0 in 200 ps 
I guess that maybe your sampling time is too short. 

Anyway I suggest to use the softcore. 

You can also avaluate the ensamble average at a number of discrete lambda
points between 0.0 and 1.0 by separate simulations for each lambda points
and then integrate them numerical to have dG.

Alessandra