From: Brian Radak (bradak_at_anl.gov)
Date: Fri Oct 06 2017 - 14:28:56 CDT
Oh, I must have missed those settings in your original email -- my
Computing an absolute pKa is difficult, but you are on the right path.
In general you are best off computing a relative pKa by also computing
the free energy of a reference compound -- there is ample literature on
this subject. In both cases you would need to use the setup that you
If you only have a single output file, then I guess parseFEP cannot do
the cumulative calculation that I described. It would be easiest to run
multiple segments, say with outputs
segment0.fep, segment1.fep, segment2.fep
and then run parseFEP by loading:
2) segment0.fep + segment1.fep
3) segment0.fep + segment1.fep + segment2.fep
and see how the answer changes. Ideally you would have several such
points at a medium sized interval (100s of ps, say). Again, there is no
established method that solves the problem and this is just one
recommended best practice.
On 10/06/2017 03:19 PM, Sadegh Faramarzi Ganjabad wrote:
> I just simply want to calculate pKa of a GLU residue in a membrane
> protein interior from DeltaG. I already have extra bonds on.
> Also, alchVdwLambdaEnd is already 1.0, with alchBondDecouple on.
> Apparently, ParseFEP does not split the FEP output file.
> On Fri, Oct 6, 2017 at 3:11 PM, Brian Radak <bradak_at_anl.gov
> <mailto:bradak_at_anl.gov>> wrote:
> Hm - that might be a bug in the output; I'll have to check.
> By default the bonded terms are not scaled, since this is not
> usually what people want to do.
> If you want to compute a relative free energy, then in your case
> you should set:
> alchBondLambdaEnd 1.0
> alchBondDecouple on
> and then add extraBonds.
> I have no idea if that is what you want to compute -- why are you
> running this simulation?
> parseFEP should very happily take different amounts of input,
> although I'm not sure that it permits subsampling of the loaded
> files. It also assumes that the number of samples from all lambdas
> are equivalent.
> You would have to load some files, run the BAR calculation, then
> load additional files and run the BAR calculation again.
> On 10/06/2017 03:04 PM, Sadegh Faramarzi Ganjabad wrote:
>> Thanks replying.
>> 1) Regarding bonded terms, I thought they are already included;
>> in the log files I have these lines
>> FEP CURRENT LAMBDA VALUE 0
>> Info: FEP ELEC. ACTIVE FOR ANNIHILATED PARTICLES BETWEEN LAMBDA =
>> 0 AND LAMBDA = 0.5
>> Info: FEP ELEC. ACTIVE FOR EXNIHILATED PARTICLES BETWEEN LAMBDA =
>> 0.5 AND LAMBDA = 1
>> Info: FEP VDW ACTIVE FOR ANNIHILATED PARTICLES BETWEEN LAMBDA = 0
>> AND LAMBDA = 1
>> Info: FEP VDW ACTIVE FOR EXNIHILATED PARTICLES BETWEEN LAMBDA = 0
>> AND LAMBDA = 1
>> Info: FEP BOND ACTIVE FOR ANNIHILATED PARTICLES BETWEEN LAMBDA =
>> 0 AND LAMBDA = 1
>> Info: FEP BOND ACTIVE FOR EXNIHILATED PARTICLES BETWEEN LAMBDA =
>> 0 AND LAMBDA = 1
>> FEPTITLE: TS BOND2 ELECT2 VDW2
>> FEP: 0 22006.7270 -93162.0780 2413.1778
>> If not, how could I include them?
>> 2) How necessary it is to compute a difference in relative
>> solvation free energies?
>> 3) Is there a way I could do study long time behavior of the BAR
>> estimate of the cumulative data with the ParseFEP plugin? I
>> assume I should split the entire simulation to smaller fractions?
>> On Thu, Oct 5, 2017 at 4:30 PM, Brian Radak <bradak_at_anl.gov
>> <mailto:bradak_at_anl.gov>> wrote:
>> Does your transformation also include bonded terms? You will
>> not compute a relative ionization free energy if you do not
>> include those terms, although you can, in principle, compute
>> a difference in relative solvation free energies by doing the
>> corresponding calculation in the gas phase. I don't know your
>> reasons for running your simulation, but I feel that I should
>> point this out -- there are new keywords as of NAMD 2.11 for
>> doing this (see the user guide).
>> Agreement between the forward and backward FEP calculations
>> is only a necessary condition for FEP convergence -- it is
>> neither a sufficient condition for convergence nor a
>> meaningful condition for other estimators, such as BAR.
>> Another necessary, but still not sufficient, condition for
>> convergence might be found in the long time behavior of the
>> BAR estimate of the cumulative data (i.e. use the first 10
>> ns, then the first 20 ns, etc. and see if the answer changes
>> as the simulation gets longer).
>> The forward and backward construction loses all meaning when
>> analyzing with thermodynamic integration, for the same reason
>> as it does for BAR.
>> On 10/05/2017 02:36 PM, Sadegh Faramarzi Ganjabad wrote:
>> I am running a test FEP simulation to calculate free
>> energy of the ionization of an amino acid inside a
>> membrane protein. Here is the FEP part of my namd config.
>> source ../tools/fep.tcl
>> alch on
>> alchType FEP
>> alchFile all-ion-beta.pdb
>> alchCol B
>> alchOutFile forward-shift-pr.fepout
>> alchOutFreq 10
>> alchBondDecouple on
>> alchBondLambdaEnd 1.0
>> alchVdwLambdaEnd 1.0
>> alchElecLambdaStart 0.5 #soft potential
>> alchVdWShiftCoeff 6.0
>> #alchDecouple off
>> alchEquilSteps 100
>> set numSteps 500
>> runFEPmin 0.0 1.0 0.0625 10000 1000 310
>> Then I used ParseFEP of VMD. The difference between
>> DeltaG of the forward and backward transformation is
>> about 5 kcal/mol for lambda=0 to 0.5, and it become
>> almost 0 for lambda = 0.6 to 1. I should probably
>> increase simulation time for each lambda from 10,000 to
>> say 10,000,000 until the plots of forward and backward
>> lie on top of each other. My question is that is this the
>> only criterion for convergence of FEP calculation? How
>> can I make sure my simulation is long enough?
>> Also, if I use thermodynamic integration method for free
>> energy calculations, what would be convergence criteria
>> for that?
>> Brian Radak
>> Postdoctoral Appointee
>> Leadership Computing Facility
>> Argonne National Laboratory
>> 9700 South Cass Avenue, Bldg. 240
>> Argonne, IL 60439-4854
>> (630) 252-8643 <tel:%28630%29%20252-8643>
>> brian.radak_at_anl.gov <mailto:brian.radak_at_anl.gov>
> Brian Radak
> Postdoctoral Appointee
> Leadership Computing Facility
> Argonne National Laboratory
> 9700 South Cass Avenue, Bldg. 240
> Argonne, IL 60439-4854
> (630) 252-8643 <tel:%28630%29%20252-8643>
> brian.radak_at_anl.gov <mailto:brian.radak_at_anl.gov>
-- 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 : Sun Dec 31 2017 - 23:21:42 CST