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

[Leontyev Igor]

Dear Gerrit,
The fist issue from which I started investigation of this question was
indeed: "What interactions are excluded from the energy?"  To simplify the
question I considered the model system (two point charges, see original
message) with no unperturbed atoms and made sure that there is no excluded
energy groups (in mdp-file) and excluded interactions (in topology). Each
point charge is specified as separate molecule. Thus, the interactions  for
which "charges of both atoms involved in the interaction are scaled by
lambda," should be present in reported below energies. I don't use QM/MM,
the "QM-system" is just synonym of "solute" or "perturbed system."

Igor




>Explicit lambda dependence comes in if the charges of both atoms
>involved inthe interaction are scaled by lambda. What I suspect is
>that you have no interactions within the QM subsystem (i atoms), and
>there is only interaction between the perturbed i aotms and the
>environment atoms j.
>
>Do you use QM/MM, or is the name QM system just chosen like that?
>
>Best,
>
>Gerrit


>>On 24 Jun 2008, at 02:35, 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