[gmx-developers] Inconsistency in Coulomb 1-4 energy term scaled by lambda

[David Mobley]

Matt,

> I am running some sanity checks on a system where I scale the charges
>  on a peptide in water.  Specifically, I can modify the charges in one
>  of two ways:
>  A) with free-energy = yes, and having lambda scale the charges
>  B) with free-energy = no, and modifying the charges directly in the
>  itp file.  To do this, I multiply the OPLS-AA charges by the factor
>  lambda, and write the modified itp file prior to processing with
>  grompp.
>
>  In principle, these two approaches should yield the same energies --
>  am I correct in assuming this?

I am not sure; you'll probably need to check the manual on how the
pairs are generated. In particular, one thing that I think makes this
complicated is the fact that usually you have a pairs list in your
topology file that specifies 1-4 interactions, but this doesn't state
interaction parameters explicitly. Then GROMACS will go through and
"generate" (calculate) the pair interactions when you run grompp (see
the gen-pairs entry in the header in your force field itp file). It
all hinges on how these are calculated. If these are calculated from
the modified charges in the topology file (scenario B), then they may
not be the same as if they're calculated from the A and B states with
free energy on. What this means is that you would be taking a somewhat
different path through lambda space -- the endpoints are the same but
the path is somewhat different.

To put this another way, when you set up a topology file with
different A and B state charges and linear scaling, your path follows
V = (1-lambda)*V_0 + lambda*V_1. Note that this is *not* equivalent to
changing the charges: Intramolecular interactions scale as q^2, but
with this functional form they are only decreased linearly with
lambda.

On the other hand, if you modify the charges directly in the topology
file, then your Hamiltonian has a different functional form that is
not as easily written. But now you're turning off intramolecular
interactions with lambda^2  (q^2), and intermolecular interactions
(between your solute and the environment) with lambda (q). So
obviously your pathway is different.

If this is right, of course, then not only should you see significant
differences in the 1-4 energies, but also in the SR coulomb, and it
looks like you *are* seeing such differences, although you didn't
mention them.

You can of course check that this is not a problem by ensuring that
you get the same free energy using both pathways. Both are completely
fine, as far as I'm aware; they're just nonequivalent.

Best wishes,
David

>  In fact, I find that the coulomb 1-4 energy are different, depending
>  on the lambda value used.  For lambda = 0 and 1, approaches A and B
>  give me the same energies (as read in the log file).  For lambda =
>  0.5, most of the energy terms are identical (or very similar), but the
>  Coulomb 1-4 terms differ depending on whether A or B are used.  I list
>  the relevant lines of the two log files below.
>
>  SCALING TYPE A
>            Bond          Angle    Proper Dih. Ryckaert-Bell.
>  LJ-14
>       2.96803e+00    7.06877e+01    2.24532e+01    7.20134e+02
>  1.14116e+02
>        Coulomb-14        LJ (SR)   Coulomb (SR)   Coulomb (LR)
>  Potential
>       1.12253e+03    5.36139e+04   -3.49154e+05    1.12833e+03
>  -2.92359e+05
>       Kinetic En.   Total Energy    Temperature Pressure (bar)
>  dVpot/dlambda
>       5.30803e+04   -2.39278e+05    2.97292e+02    1.17644e+03
>  -2.01533e+03
>
>  SCALING TYPE B
>            Bond          Angle    Proper Dih. Ryckaert-Bell.
>  LJ-14
>       2.96803e+00    7.06877e+01    2.24532e+01    7.20134e+02    1.14116e+02
>        Coulomb-14        LJ (SR)   Coulomb (SR)   Coulomb (LR)      Potential
>       5.61266e+02    5.36139e+04   -3.48176e+05    1.12811e+03   -2.91942e+05
>       Kinetic En.   Total Energy    Temperature Pressure (bar)
>       5.30805e+04   -2.38862e+05    2.97294e+02    1.18713e+03
>
>  This issue is observed for both Reaction Field and Coulomb
>  electrostatics, and in versions 3.3.1, 3.3.2 and
>  gromacs-3.3.99_development_20071120.
>
>  I'll be happy to provide the relevant files.
>
>  Regards,
>
>  Matt
>
>  --
>  Matt Wyczalkowski
>  Doctoral Candidate, Biomedical Engineering
>  Pappu Lab: http://lima.wustl.edu
>  Washington University in St. Louis
>  _______________________________________________
>  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.
>