From: Vermaas, Joshua (Joshua.Vermaas_at_nrel.gov)
Date: Mon May 08 2017 - 17:01:05 CDT
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-unphysi
cal 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 : Mon Dec 31 2018 - 23:20:16 CST