[gmx-users] Free energy calculation of multiple solutes

[Jeff Chen]

Dear All,

I am trying to calculate the change in free energy that occurs when
charging multiple type B molecules in an A-B mixture. For example, I have a
mixture of 900 molecules of type A and 100 of type B. When coulomb
interactions are on, B is identical to A (something like TIP4P water), but
when coulomb interactions are off, each B is simply a Lennard-Jones
particle (in the topology file, I've labelled B differently and called the
molecule SOLB, but all interaction information is identical to A). To do
the free energy calculations, I perform different simulations at each
lambda point. For example, for a simulation with lambda=0, I have the
following in my .mdp file:

free_energy = yes
init_lambda = 0.0
delta_lambda = 0
foreign_lambda= 0.05
couple-lambda0=vdw-q
couple-lambda1=vdw
couple-moltype=SOLB

and when lamba=1,

free_energy = yes
init_lambda = 1.0
delta_lambda = 0
foreign_lambda= 0.95
couple-lambda0=vdw-q
couple-lambda1=vdw
couple-moltype=SOLB

I expected that when lambda=0, B molecules would behave identically to A
molecules with the coulombic energy fully turned on (couple-lambda0=vdw-q),
and when lambda=1, B molecules would behave like LJ particles, with no
coulombic interactions, either with each other, or with A molecules.

However, when I look at the coulombic energy for B-B interactions (I used
tpbconv -o to make a new *.tpr for only B molecules, then used mdrun -rerun
with a trajectory containing only information about B molecules), I find
that both the short-range and reciprocal coulombic energies are NOT zero as
I had expected for lambda = 1, and are in fact exactly equal to the
negative of the respective coulombic energies obtained using the same
trajectory data when lambda=0.
I did verify that coulombic interactions between A and B molecules are zero
when lambda=1, as expected, so there would be no problem if I only had a
single solute molecule.

Obviously, I've misunderstood something. I would greatly appreciate any
help understanding my error/misunderstanding.

Thanks in advance.
Jeff