Re: pKa calculation with thermodynamic integration (TI)

From: Sadegh Faramarzi Ganjabad (safaramarziganjabad_at_mix.wvu.edu)
Date: Wed May 17 2017 - 02:00:11 CDT

Brian,

Just wanted to make sure I get what you mean. Here is what I did; I applied
GLU2 (protonated) and GLUHD (deprotonated) patches from the library you
shared with me. Here are those patches

PRES GLU2 0.00 ! protonated GLU, proton @OE2, dummy proton @OE1
GROUP
ATOM CG CT2 -0.21
ATOM HG1 HA2 0.09 ! HG1 OE1..HE1
ATOM HG2 HA2 0.09 ! | //
ATOM CD CD 0.75 ! --CG--CD
ATOM OE1 OB -0.55 ! | \
ATOM HE1 HD 0.00 ! HG2 OE2-HE2
ATOM OE2 OH1 -0.61
ATOM HE2 H 0.44
BOND OE1 HE1 OE2 HE2
############################################################
PRES GLUHD -1.00
GROUP
ATOM CB CT2A -0.18
ATOM HB1 HA2 0.09
ATOM HB2 HA2 0.09
GROUP ! new atoms
ATOM CG1 CT2 -0.28
ATOM HG11 HA2 0.09 ! HG11 OE11..HE11
ATOM HG21 HA2 0.09 ! | //
ATOM CD1 CC 0.62 ! --CG1--CD1 (-)
ATOM OE11 OC -0.76 ! | \
ATOM HE11 HD 0.00 ! HG21 OE21..HE21
ATOM OE21 OC -0.76
ATOM HE21 HD 0.00
BOND CG1 CB CD1 CG1 OE21 CD1
BOND CG1 HG11 CG1 HG21
DOUBLE CD1 OE11
BOND OE11 HE11 OE21 HE21

Then I generated a psf/pdb with psfgen and here is what the dual Glu
residue looks like in the pdb file

ATOM 1333 N GLU P 90 -3.463 -2.808 7.915 1.00 0.00 PR

ATOM 1334 HN GLU P 90 -3.522 -2.934 6.928 1.00 0.00 PR

ATOM 1335 CA GLU P 90 -2.190 -2.281 8.337 1.00 0.00 PR

ATOM 1336 HA GLU P 90 -1.850 -2.879 9.169 1.00 0.00 PR

ATOM 1337 CB GLU P 90 -1.139 -2.375 7.236 1.00 0.00 PR

ATOM 1338 HB1 GLU P 90 -1.640 -1.943 6.343 1.00 0.00 PR

ATOM 1339 HB2 GLU P 90 -0.191 -1.808 7.355 1.00 0.00 PR

ATOM 1340 CG1 GLU P 90 -0.671 -3.792 6.841 1.00 0.00 PR

ATOM 1341 HG11 GLU P 90 -0.233 -4.452 7.620 1.00 0.00 PR

ATOM 1342 HG21 GLU P 90 -1.560 -4.346 6.469 1.00 0.00 PR

ATOM 1343 CD1 GLU P 90 0.313 -3.702 5.736 1.00 0.00 PR

ATOM 1344 OE11 GLU P 90 1.475 -3.716 5.860 1.00 0.00 PR

ATOM 1345 HE11 GLU P 90 1.925 -4.486 5.422 1.00 0.00 PR

ATOM 1346 OE21 GLU P 90 -0.301 -3.707 4.545 1.00 0.00 PR

ATOM 1347 HE21 GLU P 90 0.387 -3.549 3.894 1.00 0.00 PR

ATOM 1348 CG GLU P 90 -0.671 -3.792 6.841 1.00 0.00 PR

ATOM 1349 HG1 GLU P 90 -0.233 -4.452 7.620 1.00 0.00 PR

ATOM 1350 HG2 GLU P 90 -1.560 -4.346 6.469 1.00 0.00 PR

ATOM 1351 CD GLU P 90 0.313 -3.702 5.736 1.00 0.00 PR

ATOM 1352 OE1 GLU P 90 1.475 -3.716 5.860 1.00 0.00 PR

ATOM 1353 HE1 GLU P 90 1.925 -4.486 5.422 1.00 0.00 PR

ATOM 1354 OE2 GLU P 90 -0.301 -3.707 4.545 1.00 0.00 PR

ATOM 1355 HE2 GLU P 90 0.387 -3.549 3.894 1.00 0.00 PR

ATOM 1356 C GLU P 90 -2.249 -0.857 8.908 1.00 0.00 PR

Is that what I am supposed to do? if so, how can I tight (for example) CD
and CD1 atoms together by a harmonic restraint? at the beginning, there is
no distance between them.

Thanks,
Sadegh

On Thu, May 11, 2017 at 11:27 AM, Brian Radak <bradak_at_anl.gov> 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>
> 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;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
>

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