[gmx-users] Relative alchemical free energy calculations: some confusion about options

[Ryan Muraglia]

Dear GROMACS users,

I am attempting to carry out relative alchemical free energy calculations
for amino acid mutations. Below I've reproduced an excerpt from my topology
file for a leucine to alanine mutation (generated with the pmx package):

    17           N      2    L2A      N     17  -0.415700    14.0100
    18           H      2    L2A      H     18   0.271900     1.0080
    19          CT      2    L2A     CA     19  -0.051800    12.0100
   CT   0.033700    12.0100
    20          H1      2    L2A     HA     20   0.092200     1.0080
   H1   0.082300     1.0080
    21          CT      2    L2A     CB     21  -0.110200    12.0100
   CT  -0.182500    12.0100
    22          HC      2    L2A    HB1     22   0.045700     1.0080
   HC   0.060300     1.0080
    23          HC      2    L2A    HB2     23   0.045700     1.0080
   HC   0.060300     1.0080
    24          CT      2    L2A     CG     24   0.353100    12.0100
   HC   0.060300     1.0080
    25          HC      2    L2A     HG     25  -0.036100     1.0080
 DUM_HC   0.000000     1.0000
    26          CT      2    L2A    CD1     26  -0.412100    12.0100
 DUM_CT   0.000000     1.0000
    27          HC      2    L2A   HD11     27   0.100000     1.0080
 DUM_HC   0.000000     1.0000
    28          HC      2    L2A   HD12     28   0.100000     1.0080
 DUM_HC   0.000000     1.0000
    29          HC      2    L2A   HD13     29   0.100000     1.0080
 DUM_HC   0.000000     1.0000
    30          CT      2    L2A    CD2     30  -0.412100    12.0100
 DUM_CT   0.000000     1.0000
    31          HC      2    L2A   HD21     31   0.100000     1.0080
 DUM_HC   0.000000     1.0000
    32          HC      2    L2A   HD22     32   0.100000     1.0080
 DUM_HC   0.000000     1.0000
    33          HC      2    L2A   HD23     33   0.100000     1.0080
 DUM_HC   0.000000     1.0000
    34           C      2    L2A      C     34   0.597300    12.0100
    35           O      2    L2A      O     35  -0.567900    16.0000

My plan is to follow the general process of a three step alchemical
transformation: decharging the changing atoms, swapping the Lennard-Jones
terms, then recharging the new atoms. Looking at the available free energy
options, I'm unsure of how to proceed, and I want to clarify the function
of the couple-lambda0 and 1 flags. Here is one scheme I've thought of.
Please let me know if it is appropriate, or if I've made an error somewhere.

Step 1: use an unmutated leucine topology (not reproduced here, but
essentially same as above, without B state specified and with L2A renamed
to LEU) to calculate the decharging step.
Parameters are: couple-lambda0 = vdw-q, couple-lambda1 = vdw, coul_lambdas
= 0.0 0.5 1.0 (normally I will use more lambdas, but for clarity, I'm only
specifying three here).
My intention here is for the 0th state to be a fully interacting leucine,
and the 2nd state to be a decharged leucine. Is this correct? Is specifying
couple-lambda1 = vdw sufficient, or do I have to specify a B state in the
topology with explicit 0 charges?

Step 2: use the mutated topology above with the following parameters:
couple-lambda0 = vdw, couple-lambda1 = vdw-q,
coul_lambdas = 0.0 0.0 0.0 0.5 1.0
vdw_lambdas  = 0.0 0.5 1.0 1.0 1.0
bonded and mass lambdas same as vdw.
Here my intention is for the 0th state to be the same uncharged leucine
from before, then 2nd state to be the uncharged alanine, then the 4th state
to be the fully interacting alanine.

Again, I want to confirm that my understanding of the free energy options
is appropriate, and this calculation scheme will accomplish what I
described above.

Thank you for reading, and for any advice you can offer!

-- 
Ryan Muraglia
rmuraglia at gmail.com