Re: pKa calculation with thermodynamic integration (TI)

From: Vermaas, Joshua (Joshua.Vermaas_at_nrel.gov)
Date: Thu May 11 2017 - 11:52:48 CDT

Adding on to Brian's response, in the initial structure, psfgen will need to guess coordinates for the new atoms, or you correct them as a post-processing step. It is not a big deal, since psfgen's guesses aren't outrageous, and they will be pinned together anyway, so the positions will correctly work themselves out in my experience. For the pinning force constant, I happened to use 20 kcal/molA^2 as the force constant, since that is what the papers recommended as the "best", but it did not really matter much from what I recall, so anything on the tens to hundreds of kcal/mol A^2 scale would probably be acceptable.

-Josh

On 05/11/2017 09:27 AM, Brian Radak wrote:

Note also that atom types change, which includes Lennard-Jones and bonded terms. Unlike the AMBER and GROMACS force fields, the CHARMM force field permits these to change for different protonation states -- NAMD's "dual topology" approach is likely the easiest way to accommodate this (see the recent GROMACS constant pH paper where they dubiously modify the force field to avoid this problem).

In general, I think adding a dummy atom makes sense so that all alchemical atoms at one endpoint map to one atom at the other endpoint. For dummy atoms, so long as the attachment only involves one bond, one angle, and one dihedral, then their is a (nearly) analytic shift to the relative free energy, otherwise this will cancel for relative free energy differences (such as for pKa shifts).

The restraint between atoms in different alchemical groups will make them nearly coincident, but not rigorously so. You formally introduce new kinetic energy terms corresponding to ideal gas particles at the endpoints (depending on how you handle alchemical bonds). There is actually an open debate over whether or not this should impact the kinetic virial in NAMD (it currently does), so I would recommend not using constant pressure for your TI/FEP calculations.

I have had very good luck with a force constant as strong as 100 kcal/mol-A^2, even when the alchemical region involves a holonomic constraint (which one would think causes a bit of trouble).

Brian

On 05/10/2017 07:36 PM, Sadegh Faramarzi Ganjabad wrote:
Josh,

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.

Thanks,
Sadegh

On Mon, May 8, 2017 at 6:01 PM, Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov<mailto:Joshua.Vermaas_at_nrel.gov>> wrote:
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<https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2F3.0.CO&data=02%7C01%7CJoshua.Vermaas%40nrel.gov%7C35ccd3287a3840ee9ce508d498823d85%7Ca0f29d7e28cd4f5484427885aee7c080%7C1%7C0%7C636301132527295245&sdata=6XGUFRPgBcOaHS4n%2FQ%2FZ4Ii55FZyeIQZn3mqLuhNflA%3D&reserved=0>;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
GROUP
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
BOND CD2 CG2 CG2 HG21 CG2 HG22 CG2 CB
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

GROUP
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 !
BOND OE2 HE2

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

--
Brian Radak
Postdoctoral Appointee
Leadership Computing Facility
Argonne National Laboratory
9700 South Cass Avenue, Bldg. 240
Argonne, IL 60439-4854
(630) 252-8643
brian.radak_at_anl.gov<mailto:brian.radak_at_anl.gov>

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