Fabricio,
Unfortunately, I have no idea if the softcore potential code is compatible
with the QM/MM code. I believe the topology information is handled in a
slightly unusual way, which may cause trouble. Perhaps Thomas (or Ross?)
will see this and be able to comment. However I can't imagine it is
necessary to use the softcore code since the proton likely has no
Lennard-Jones term anyway. Even if that is not the case, I would guess that
the epsilon parameters are probably small enough that they can be scaled
away linearly, although this is pure speculation.
I would try using the normal TI code and setting lambda = 1 or 0 and
verifying that you obtain the same dynamics you did with the QM/MM code on
its own.
Regards,
Brian
On Mon, Feb 4, 2013 at 1:31 PM, Fabrício Bracht <bracht.iq.ufrj.br> wrote:
> What I meant is that I have already performed a MM ti simulation, but since
> I am using soft core potentials, I cannot set Lambda = 1 or lambda = 0.
> This means that I have not done such simulations (lambda = 0 or 1) with
> qmmm either. Setting restraints to the bonds between the zinc atom and the
> hydroxyl group has not worked either. The addition of these restraints
> result in very abnormal behaviour, such as very distorted geometry of the
> entire qm region, as well as SCF failing to converge. What I still don't
> understand is why such abnormal behaviour happens only with TI simulations.
> As I said before, if I simulate only the hydroxyl system (no multisander
> simulation with the HOH system) with qmmm, the simulation runs smoothly. Do
> you imagine I will have more success doing Ti simulation without softcore
> potentials?
> Fabrício Bracht
>
> 2013/2/3 Brian Radak <radak004.umn.edu>
>
> > 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
> >
> _______________________________________________
> 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 Mon Feb 04 2013 - 11:00:05 PST
