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
rmurag...@gmail.com
--
Gromacs Users mailing list
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a
mail to gmx-users-requ...@gromacs.org.
