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:28:56 CDT

Oh, I must have missed those settings in your original email -- my
apologies.

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
currently have.

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:

1) segment0.fep
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.

Brian

On 10/06/2017 03:19 PM, Sadegh Faramarzi Ganjabad wrote:
> Brian,
>
> 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.
>
> Thanks,
> Sadegh
>
> 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.
>
> 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
>> <https://maps.google.com/?q=9700+South+Cass+Avenue,+Bldg.+240&entry=gmail&source=g>
>> 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
> <https://maps.google.com/?q=9700+South+Cass+Avenue,+Bldg.+240&entry=gmail&source=g>
> 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 : Mon Dec 31 2018 - 23:20:38 CST