Hi Brian. Thank you for the advice. Indeed I had already seen the article,
but I think it is time to read it again with more attention. The lambda
values I am using are 0.01592; 0.08198; 0.19331; 0.33787; 0.5; 0.66213;
0.80669; 0.91802; 0.98408. These are the values listed in amber12 manual
page #117.
The odd thing though is the difference I get in behaviour when I run the
simulation still with qmmm and using the same basic input parameters; i.e.,
using the same input file as for the TI calculation, but without the TI
part. All atoms stay bound, no matter what happens. At least for the
hydroxyl ligand. I was, indeed, expecting some problems to arise with the
water molecule, since an oxygen atom bound to two hydrogens and a zinc atom
is unlikely to stay that way for long.
You mentioned doing a MM simulation. Well, I already did that, but I am
not sure I understood what you meant by "perform MM to QM/MM
transformations on the endpoints as a kind of correction."
Thank you again
Fabrício Bracht
2013/2/2 Brian Radak <radak004.umn.edu>
> Hi Fabricio,
>
> Can you clarify what exactly you mean by "start with the TI simulation"?
> What values of lambda? Does this behavior occur for all of them?
>
> Some other things that might be informative: What topolog(y/ies) are you
> using? What bonded terms exist on the protons?
>
> I believe the best established QM/MM protocols for pKa type calculations
> involve some kind of restraints, but how you choose to apply them is still
> something of an art. The following paper describes a protocol in CHARMM,
> but I'm pretty sure the ideas transfer well to AMBER:
>
> Li and Qiang, J Phys Chem B 2003, 107, 14521.
>
> I attempted similar calculations a while back but had a very difficult time
> getting reproducible results. I suspect that there are some zero of energy
> issues that need to be handled when you have the same number of electrons
> but different orbitals. I am not sure what the solution is there. An
> alternative path is to perform the proton disappearing as an MM to MM
> mutation and then perform MM to QM/MM transformations on the endpoints as a
> kind of correction. This is more work but I suspect will converge faster in
> most situations.
>
> Regards,
> Brian
>
> On Sat, Feb 2, 2013 at 2:35 PM, Fabrício Bracht <bracht.iq.ufrj.br> wrote:
>
> > Hi all. I am doing a TI caltulation in order to obtain the free energy
> > values associated with the loss of a proton from a hydroxyl group bound
> to
> > a zinc atom coordinated with 3 other aminoacids. System 1 comprises of a
> > water molecule bound to the zinc atom, and system 2 is the hydroxyl
> group.
> > When I run qmmm simulations, using DFTB, of both systems separately, the
> > system stays intact, i.e. no bonds are broken and both the water and the
> > hydroxyl stay bound to the zinc atom ( as was expected). But as soon as I
> > star with the TI simulation, the hydroxyl group always drifts away.
> > Sometimes the oxygen atom stays and the hydrogen goes for a walk. This
> > happens only on TI simulations and the water system stays intact. I have
> > tried different sizes of qm region, different cutoffs, tried several
> > different starting structures ( minimized, equilibrated, thermalized etc)
> > and the result is almost always the same. Below is the input file for the
> > hydroxyl group system. I thought about restraining the hydroxyl group
> using
> > the nmropt module. Can I use the dvdl_norest flag to ignore those
> > restraints in the energy calculation? Please let me know if you need any
> > other detail from the simulation and/or if you need/want any file from
> the
> > inputs.
> > Thank you
> > Fabrício Bracht
> >
> > TI of hydroxyl group step 1
> > &cntrl
> > imin = 0,
> > irest = 1,
> > ntx = 7,
> > ntb = 2, pres0 = 1.0, ntp = 1, taup = 2.0,
> > cut = 8.0,
> > ntr = 0,
> > ntc = 2,
> > ntf = 1,
> > temp0 = 298.0,
> > ntt = 3,
> > gamma_ln = 1.0,
> > nstlim = 1000000, dt = 0.0005, ntave = 100,
> > ntpr = 100, ntwx = 100, ntwr = 100,
> > ig = 10703, ioutfm = 1, iwrap = 1,
> > icfe = 1, ifsc = 1, clambda = 0.01592,
> > ifqnt = 1, scmask = ':342', idecomp = 0, nmropt = 0,
> > /
> > &qmmm
> > qmmask=
> > ':200,255,258&!.CA,C,HA,O,N,HN,H|:341,342|(:202,204&!.CA,C,HA,O,N,H)'
> > dftb_3rd_order = 'PA'
> > qmcharge=-1,
> > qm_theory='DFTB',
> > qmshake=0,
> > qm_ewald=1, qm_pme=1
> > qmcut=9.0
> > writepdb=1
> > /
> > Receptor residues
> > RRES 1 7853
> > END
> > Printing
> > RES 1 342
> > END
> > END
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
> >
>
>
>
> --
> ================================ Current Address =======================
> Brian Radak : BioMaPS
> Institute for Quantitative Biology
> PhD candidate - York Research Group : Rutgers, The State
> University of New Jersey
> University of Minnesota - Twin Cities : Center for Integrative
> Proteomics Room 308
> Graduate Program in Chemical Physics : 174 Frelinghuysen Road,
> Department of Chemistry : Piscataway, NJ
> 08854-8066
> radak004.umn.edu :
> radakb.biomaps.rutgers.edu
> ====================================================================
> Sorry for the multiple e-mail addresses, just use the institute appropriate
> address.
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sat Feb 02 2013 - 18:00:02 PST
