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
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Mar 14 2017 - 11:00:03 PDT
