Re: convergence criterion of free energy perturbation (FEP) simulations for pKa calculation

From: Brian Radak (bradak_at_anl.gov)
Date: Fri Oct 06 2017 - 14:11:27 CDT

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.

HTH,
Brian

On 10/06/2017 03:04 PM, Sadegh Faramarzi Ganjabad wrote:
> Brian,
>
> 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?
>
> Best,
> Sadegh
>
> 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.
>
> HTH,
> Brian
>
> On 10/05/2017 02:36 PM, Sadegh Faramarzi Ganjabad wrote:
>
> Hello,
>
> 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?
>
> 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 <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