Re: pKa calculation with thermodynamic integration (TI)

From: Sadegh Faramarzi Ganjabad (
Date: Wed May 10 2017 - 19:36:16 CDT


Thanks for your detailed instructions. Correct me if I'm wrong, so all
terminal atoms of GLU disappear and new atoms appear at same positions but
with new partial charges (except the titratable hydrogen that vanishes
during the perturbation). Just few questions; 1) In the initial structure,
will CG and CG2, CD and CD2 , etc. positioned on top of each other having
same coordinates? 2) What is the best value for the harmonic constraint? I
saw 20 kcal/mol in two of papers you mentioned.


On Mon, May 8, 2017 at 6:01 PM, Vermaas, Joshua <>

> Hi Sadegh,
> The "dual topology" part of NAMD's free energy methods is actually a bit
> annoying in your case, since you need two copies of each atom that changes
> in some way, which includes the charges. What I did when I was doing this
> in my own research was create a patch to duplicate the glutamate atoms that
> changed (see below), calling them CG2 instead of CG, CD2 instead of CD,
> etc. Then you cause one set to disappear and the other set to appear by
> tagging them appropriately in your FEP or TI setup. Note that for faster
> convergence, you should "pin" together corresponding atoms (CG to CG2, CD
> to CD2, etc) with a weak harmonic restraint (I used the extrabonds feature
> of NAMD). See 10.1002/(SICI)1096-987X(199808)19:11<1278::AID-JCC7>3.0.CO;2-H
> or 10.1021/jp807701h. The reason you can get away with this extra restraint
> is that you are basically imposing a single-topology style on the parts
> where only the charges are changing, rather than actual atomtypes. If you
> don't, the nearly decoupled piece will explore super-unphysical
> configurations, and the convergence at the ends of your lambda space gets
> really slow. I made this mistake once (10.1021/acs.biochem.5b00033, you
> want the discussion on page 2107), and my results suffered because of it
> until I dug up those old papers and corrected it.
> In terms of comparing TI and FEP, I personally find TI much easier to deal
> with, since if it turns out you aren't converged yet, you just make the
> simulations for each window longer, while the AlchFEP plugin works best if
> the FEP is carried out as one long stepwise simulation, which means to get
> more sampling you need to redo what you have already done. The nuts and
> bolts of running TI is the same as FEP, you just change the alchType
> keyword.
> -Josh
> PRES FDP 0.00 ! patch for protonated glutamic acid, suitable for
> FEP simulations
> ATOM CG2 CT2 -0.21 !
> ATOM HG21 HA2 0.09 ! HG12 OE21
> ATOM HG22 HA2 0.09 ! | //
> ATOM CD2 CD 0.75 ! -CG2--CD2
> ATOM OE21 OB -0.55 ! | \
> ATOM OE22 OH1 -0.61 ! HG22 OE22-HE22
> ATOM HE22 H 0.44 !
> BOND OE22 HE22 CD2 OE22 CD2 OE21
> IMPR CD2 CG2 OE22 OE21
> !Original GLU charges and atom names.
> !ATOM CG CT2 -0.28
> !ATOM HG1 HA2 0.09
> !ATOM HG2 HA2 0.09
> !ATOM CD CC 0.62
> !ATOM OE1 OC -0.76
> !ATOM OE2 OC -0.76
> On 05/08/2017 02:04 PM, Sadegh Faramarzi Ganjabad wrote:
> Dear all,
> I am planning to calculate pKa of a Glu in a membrane protein with
> thermodynamic integration method. As you know, CHARMM 36 has parameters for
> both protonated and deprotonated Glu. However, there is no NAMD tutorial on
> TI. I had few questions about my system setup. I assume that the dual
> topology should be the same as free energy perturbation (FEP). And the
> following terminal group of Glu needs to be changed during the perturbation
> process
> ATOM CG CT2 -0.21 !
> ATOM HG1 HA2 0.09 ! HG1 OE1
> ATOM HG2 HA2 0.09 ! | //
> ATOM CD CD 0.75 ! -CG--CD
> ATOM OE1 OB -0.55 ! | \
> ATOM OE2 OH1 -0.61 ! HG2 OE2-HE2
> ATOM HE2 H 0.44 !
> Then, is this group supposed to vanish during the perturbation and a new
> group should appear at the same position, with HE2 omitted and updated
> partial charges for the rest of atoms? or only HE2 vanished and the charges
> of other atoms are updated without vanishing/appearing? for the latter I am
> not sure what the dual topology would look like (with keeping the atoms at
> same positions and only reassigning partial charges).
> Also, I am not sure what differences would be between FEP and TI in terms
> of simulation procedure. I know the theory of each method, but compared to
> FEP, there is no clear instructions about running TI as I said. In
> 'fep.tcl' script of FEP tutorial files there is a section for TI but I
> don't know how to use it. Any help is highly appreciated.
> Thanks,
> Sadegh

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2018 - 23:20:17 CST