Fabricio,
So do you mean that setting lambda = 1.0 or 0.0 with the TI code gives
different results than the QM/MM code on its own? That would be concerning
indeed. I'm pretty sure that, all else equal, running an NVE simulation
(ntt = 0) or setting ig to an explicit value should result in identical
trajectories for both of those circumstances, even with multiple
processors. That is the most simplest and most obvious test that I can
think of (nstlim = 10, ntpr = 1, etc. should be adequate).
As for the alternate transformation method, I am not sure how much this is
done in the literature, but the idea would be to have the proton
"disappear" via an MM transformation with an appropriate force field (in
principle this can be arbitrary). Then, one would perform two additional
transformations from the MM endpoints where one endpoint is the MM
description and the other is the analogous QM/MM description. This should
in principle amount to a smaller perturbation of the bonded and
electrostatic terms (no guarantees!) as well as avoid the "problem" of
having weird bonded terms on the QM atoms.
I have no idea if the latter suggestion would solve your instability
issues, but if no other solutions present themselves it's probably worth a
try.
Regards,
Brian
On Sat, Feb 2, 2013 at 8:33 PM, Fabrício Bracht <bracht.iq.ufrj.br> wrote:
> 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
>
--
================================ 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
Received on Sun Feb 03 2013 - 09:30:03 PST
