[gmx-users] Rugged dG/dlambda when turning off charges

[Robert Johnson]

Hello everyone,
As I mentioned in a previous post, I'm trying to compute the free
energy of binding of a DNA base on a carbon nanotube. My system
consists of a single DNA nucleotide (base, sugar, phosphate group) on
a rigid carbon nanotube in aqueous solution. A single Na+ counterion
is also included for charge neutralization. I'm proceeding with this
by computing the free energy associated with the following
transformations:

Water + Nanotube + DNA + Na+ -> Water + Nanotube + Nothing
Water + DNA + Na+ -> Water + Nothing

For each of these transformations I FIRST turn off the charges on the
base and counterion over a set of 5 lambda values (0.0, 0.25, 0.5,
0.75, 1.0). From literature and other posts on the GMX mailing list,
it seems like this small set of values should be adequate because
dG/dlambda tends to be pretty smooth for discharging the molecule.

Unfortunately, that doesn't seem to be the case in my system. Here is
a table of the dG/dlambda values I get from running 5 ns trajectories
of my system:

lambda     <dG/dlambda>
0.0           6211.105
0.25         1055.254
0.5           1230.675
0.75         1128.359
1.0           756.2904

If you plot this you will see a large drop from lambda=0 to 0.25. Then
dG/dlambda increase slightly from lambda=0.25 to 0.5 and then smoothly
decreases after that. This seems like a pretty strange result. Also,
it seems that something weird is going on around lambda=0. If I
compute dG/dlambda for lambda=0.05, I get something around 1200-1300,
which is again a very sharp decrease from the value of 6211.105 I get
at lambda=0. Thus, it doesn't seem like simply adding more lambda
values is going to help me.

Here are the values I'm using for the free energy calculation:
free_energy         =  yes
init_lambda         =  0.05
delta_lambda        =  0
sc_alpha            =  0.5
sc_power            =  1
sc_sigma            =  0.3

Also, here is an excerpt from my topology file:
[ atoms ]
1       P       1       DG      P       1       1.1659  30.9700
 P       0.0000  30.9700
2       O       1       DG      O1P     1       -0.7761 16.0000
 O       0.0000  16.0000
3       O       1       DG      O2P     1       -0.7761 16.0000
 O       0.0000  16.0000

I'm not perturbing any of the bonding interactions or VdW parameters -
I'm just turning off all charges on the molecule.

Does anyone have any ideas or suggestions about how this can be improved?

Thanks,
Bob Johnson