[gmx-developers] What is dH/dl functional for Coulomb interaction?

[Leontyev Igor]

Dear gmx-developers,
I am confused by the fact that dhdl-values calculated with the same
molecular structure for different values of the Coulomb interaction coupling
parameter lambda (L) does not depend on the lambda. According to the Gromacs
3.3 Manual the L-dependence of Coulomb part of Hamiltonian has form given by
Eq.(4.103) such as dH(L)/dL should explicitly depend on lambda. However, my
tests show opposite. I have tried to calculate a charging free energy of
QM-system (solute) for two different systems:
1) QM-system (solute) in protein. The QM-system consist of metal center
alone with imidozole rings of ligating residues (35 atoms) and the protein
is represented by 11938 atoms. The QM-system has non-zero atomic charges in
both the initial and final states. Total charge difference of the system
between final and initial state is -1.
2) Solute alone (no solvent). The model QM-system is constituted by two
point charges separated by 1nm . Point charges have values (0,0) in the
initial state and (-1,1) in the final state.

The independence of dhdl-values on the parameter lambda (L) was observed for
both systems. For the model system (2 point charges), dhdl-values  can be
calculated manually. According to Eq. (4.103), H(L)=f/(1nm)*L*(-1)*L*(1),
thus:
      H(L)=-138.935*L*L                        (1)
and
      dH(L)/dL=-2*138.935*L
However, results obtained by Gromacs 3.3.2 are not described by the above
functional (1). The results are:
    L              Potential energy             dhdl
  0.00000         0.00000e+00     -1.38935e+02
  0.50000       -6.94677e+01      -1.38935e+02
  1.00000       -1.38935e+02      -1.38935e+02

It is interesting to notice that these results are precisely described by
the conventional functional form for Coulomb coupling:
           H(L)=HA+L*(HB-HA),                        (2)

where HB and HA are Hamiltonians in the final and initial states,
respectively. Thus, could someone confirm that the Coulomb coupling actually
implemented in Gromacs is given by functional (2) instead of fuctional given
by Eqs (4.103)-(4.106) of the Manual ?
Thanks,
Igor