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

[Berk Hess]

Hi,

All the equations for free-energy Coulomb interaction in the manual are 
incorrect.
The correct linear Coulomb interpolation is:
V=1/(4 pi eps0 rij) ((1-lambda) qAi qAj + lamda qBi qBj)
dV/dl = 1/(4 pi eps0 rij) (-qAi qAj + qBi qBj)

So there is no explicit lambda dependence for dV/dl
(unless you use soft-core interactions).

Equation 4.120 gives the correct dV/dl for soft-core interactions
and also without soft-core when you set alpha=0.

We will correct the equations for the 4.0 manual.

Berk.

Leontyev Igor wrote:
> 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
>
>
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the www 
> interface or send it to gmx-developers-request at gromacs.org.
>
> This email was Anti Virus checked by Astaro Security Gateway. 
> http://www.astaro.com
>