From: Radak, Brian K (bradak_at_anl.gov)
Date: Tue May 09 2017 - 10:41:52 CDT
I will just second that the protocol described by Josh is generally the best one for NAMD, if not for pKa calculations in general. This is more or less the route we are using to implement constant pH, which utilizes the exact same transformation.
In general, the FEP and TI codes are _identical_ in the way that interactions are computed. The only differences are in the on-the-fly analysis. If you want to join the current century and use a multi-state reweighting technique, then, in my opinion, the TI output is actually a better bet. The exception to this is the ParseFEP plugin for VMD, which has a lot of clever tools and simplifies the analysis for newcomers.
Some useful new features in 2.12 that may not be widely known:
1) You can now turn on alchemical scaling for bonded terms using alchBondLambdaEnd and alchBondDecouple - these are essential if you want to compute rigorously correct relative free energies as opposed to relative free energy differences (i.e. check your thermodynamic cycle!)
2) The vdW switching schemes are now updated and correct some inconsistent behaviors from previous versions - it may not be possible to compare exactly with older versions of the code, although I doubt the differences are visible within statistical error.
3) There are also a number of new psfgen features that may help in constructing your alchemical system. Namely that you can now alter the bfactor columns and velocities from the Tcl interface (use the new "psfset" command). This may only be in the latest CVS at this point.
4) I have a ready made library of "alchemical patches" for use with the C36 protein force field that I'd be happy to share, although I can't guarantee all of the new dummy atom parameters at this point and they might not exactly fit your needs. Contact me off list if this is something you are interested in.
Cheers,
Brian
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
________________________________________
From: owner-namd-l_at_ks.uiuc.edu [owner-namd-l_at_ks.uiuc.edu] on behalf of Vermaas, Joshua [Joshua.Vermaas_at_nrel.gov]
Sent: Monday, May 08, 2017 5:01 PM
To: namd-l_at_ks.uiuc.edu; Sadegh Faramarzi Ganjabad
Subject: Re: namd-l: pKa calculation with thermodynamic integration (TI)
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 decouple!
d 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
This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:21:16 CST