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

From: Sadegh Faramarzi Ganjabad (safaramarziganjabad_at_mix.wvu.edu)
Date: Fri Oct 06 2017 - 14:19:56 CDT

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> 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> 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
>> 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
> brian.radak_at_anl.gov
>

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2018 - 23:20:38 CST