Hi Jason,
Many thanks for your quick reply, that makes more sense to me now.
Thanks a lot for your help!
Eric
On 14 March 2017 at 17:53, Jason Swails <jason.swails.gmail.com> wrote:
>
>
> On Tue, Mar 14, 2017 at 10:40 AM, Eric Lang <eric.lang.bristol.ac.uk>
> wrote:
>
>> Hello,
>>
>> I am actually implementing the changes in titratable_residues.py only
>> now, and I think I am not 100% sure about what I should do.
>> Jason, you said that :
>>
>> *Knowing that the total reference energy is*
>>
>> *Ref Energy = TI - R*T*ln(10)*pKa*
>>
>> *If we want to solve for the "correct" TI energy, you can isolate that by
>> calculating*
>>
>> *TI = Ref Energy + R*T*ln(10)*pKa*
>>
>> *So in your case, take the Reference energy from the cpin file, add
>> 0.00199*300*log(10)*4.82, and that value is what the TI energy *should*
>> have been. You can then insert that value back into titratable_residues.py
>> where refene2.solvent_energies() is set for the GL4 residue.*
>>
>>
>> So in my cpin file, the Reference energy is STATENE right? So in my case
>> I have STATENE=0.000000,21.248624,21.248624,21.248624,21.248624,
>> So it means that the reference energy I should use is 21.248624?
>>
>
> ​Yes.
> ​
>
>> If we do *Ref Energy + R*T*ln(10)*pKa =* 21.248624+0.00199*300*log(10)*4.82= 27.8744
>> so the TI energy should have been 27.8744 to give me the correct pKa, is
>> that correct?
>>
>> Here is an extract of the glutamate entry in titrable_residues.py:
>>
>> # Glutamate
>> refene1 = _ReferenceEnergy(igb1=0, igb2=0, igb5=0, igb7=0, igb8=0)
>> refene1.solvent_energies()
>> refene1.dielc2_energies(igb1=0, igb2=0, igb5=0, igb7=0, igb8=0)
>> refene1.dielc2.solvent_energies()
>> refene2 = _ReferenceEnergy(igb1=3.89691326, igb2=8.4057785,
>> igb5=8.0855764,
>> igb7=5.305949, igb8=8.3591335)
>> refene2.solvent_energies(igb2=15.20019319)
>> refene2.dielc2_energies(igb2=3.455596, igb5=3.957270)
>> refene2.dielc2.solvent_energies()
>> refene2.set_pKa(4.4, deprotonated=False)
>>
>>
>>
>> So in order to "correct" the TI energy in titrable_residues.py, I simply
>> need to replace refene2.solvent_energies(igb2=15.20019319) with
>> refene2.solvent_energies(igb2=27.8744) and that should give me the
>> correct pKa when I re-run the simulation. Is that it? I feel like I am
>> missing something... Especially because it is a quite large difference.
>>
>
> ​Ack. I got my signs crossed, which is incredibly easy to do. If you
> subtract that factor instead of add it, you'll get
>
> 21.248624 - 0.00199*300*ln(10)*4.82 = ​14.6228
>
> Your intuition was right -- a small pKa shift should result in a
> correspondingly small change in the reference energy. So in this case the
> change is only 0.6 kcal/mol.
>
> If I want to do the same for Lys were I got a pKa of 10.38 and STATENE=-0.845507,0.000000,
>> I should do TI = -0.845507 + 0.00199*300*log(10)*10.38= 13.4233 and then
>> I would need to replace refene1.solvent_energies(igb2=-15.1417977) with
>> refene1.solvent_energies(igb2=13.4233)? That doesn't look correct.
>>
>> Shouldn't it be actually a subtraction in both cases, i.e. *Ref Energy -
>> R*T*ln(10)*pKa*, in which case it would give an energy of 14.6228
>> and -15.1143 for Glu and Lys respectively, which would be much more in line
>> with the original refene2.solvent_energies of15.20019319 and
>> -15.1417977 respectively?
>>
>
> ​Yes. There are so many differences taken here that it's easy for me to
> get confused. And if you try both + and -, the correct one to use is
> obvious (as you found out here), so I've gotten into the bad habit of doing
> both and just taking the correct one.​
>
> HTH,
> Jason
>
> --
> Jason M. Swails
>
--
Eric Lang
BrisSynBio Postdoctoral Research Associate Modelling
Centre for Computational Chemistry
School of Chemistry - University of Bristol
Bristol BS8 1TS - United Kingdom
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Mar 15 2017 - 02:00:02 PDT
