[gmx-users] free energy calculations

[CRISTIAN DANCIULESCU]

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.

Could someone tell me what is wrong in my simulations or give me 
some hints concerning the dummy atoms in the free energy 
calculations (when create or deleate atoms during a simulation)?
Thank you in advance,
Cristian