From: Brian Radak (bradak_at_anl.gov)
Date: Thu May 11 2017 - 10:27:25 CDT
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
> <http://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:17 CST