Re: [AMBER] Problem with running QM/MM - MD simulation using AM1/d-PhoT method in AMBER12

From: Marc van der Kamp <marcvanderkamp.gmail.com>
Date: Wed, 14 Nov 2012 20:46:59 +0000

Hi all, esp. Gustavo and Brian,

I have been helping Vojtech with his project. Let me point out first of all
that, as far as I'm aware, AM1/d-PhoT was developed in the group of Darrin
York especially for the same system that he is working on:
J Am Chem Soc. <http://www.ncbi.nlm.nih.gov/pubmed/18345664#> 2008 Apr
9;130(14):4680-91. Epub 2008 Mar 18.
Quantum mechanical/molecular mechanical simulation study of the mechanism
of hairpin ribozyme catalysis.
Nam K<http://www.ncbi.nlm.nih.gov/pubmed?term=Nam%20K%5BAuthor%5D&cauthor=true&cauthor_uid=18345664>
, Gao J<http://www.ncbi.nlm.nih.gov/pubmed?term=Gao%20J%5BAuthor%5D&cauthor=true&cauthor_uid=18345664>
, York DM<http://www.ncbi.nlm.nih.gov/pubmed?term=York%20DM%5BAuthor%5D&cauthor=true&cauthor_uid=18345664>

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2655239/
(and reference therein)

So the method *should* be perfectly applicable to this system, if anything.
Minimizations in AMBER do indeed give a reasonable potential energy
surface, but the reason to use a semi-empirical method is of course to do
some QM/MM MD.

Also, we have been using (almost) the same QM region as the above paper, as
well as slightly larger ones, and the problem still occurs.
It seems that in the context Nam & York have used the AM1/d-PhoT
parameters, it works fine. This was in CHARMM (with a spherical system with
stochastic boundary conditions), with the CHARMM27 parameters, *and* with
GHO-boundary atoms.
This is the reason that we also tried with the CHARMM27 parameters (using a
CHAMBER prmtop). However, the instability problems still persist. It would
be surprising to me if the GHO-boundary atoms vs. link-atoms would make the
crucial difference... But perhaps they are?
The boundary conditions are a difference as well of course, but again not
very likely that this causes the extreme QM-region instability. (Vojtech is
currently running an extra test with qm_ewald=0 and a largish qmcut value
to see if that makes a difference...)

Any ideas welcome!

All the best,
Marc


On 14 November 2012 19:37, Brian Radak <radak004.umn.edu> wrote:

> Hi all,
>
> Having used AM1/d-PhoT rather extensively, I have found it to be somewhat
> numerically sensitive to physically large (as in spread out) QM regions. I
> don't know if this is unique to the d-orbital implementation (which is
> different from the normal sp-orbital NDDO code) or a general characteristic
> of NDDO. I would give PM6 or MNDO/d a try, as those should also use
> d-orbitals for phosphorous (although I would have to double check that in
> the code to be sure).
>
> Regards,
> Brian
>
> P.S. The are virtually no assurances that AM1/d-PhoT is a good model for
> nucleobases, as these are entirely outside of the training set. Although
> the nitrogen parameters are untouched from AM1, the changing of the C, O,
> and H parameters certainly could have unbalanced some very important
> effects when it is present. That being said, I'm pretty sure that
> AM1/d-PhoT predicts superior proton affinities (wrt experiment) for small
> organic compounds compared to AM1, which might salvage some of the effects
> you are looking for.
>
> On Wed, Nov 14, 2012 at 1:41 PM, Gustavo Seabra <gustavo.seabra.gmail.com
> >wrote:
>
> > Hi,
> >
> > On Nov 14, 2012, at 2:15 PM, Vojtech wrote:
> >
> > > I have also tested another two SE methods (SCC-DFTB and AM1) and did
> not
> > > observe that problem.
> >
> > Ok, this *may* be an indication that the SE method you've chose does not
> > work very well for your system.
> >
> > > Furthermore, could you please explain me little bit more what do you
> > mean by
> > > 'parametrization of QM region'?
> >
> > Although its not very clear, sander *does* need parameters for the QM
> > region. First, because there are steps (as I mentioned in the other
> e-mail)
> > that are usually done with the system full-MM. Specifically, the initial
> > energy minimization and thermalization.
> >
> > Secondly, there are QM/MM interactions that use MM parameters: the vdW
> > interactions between QM-MM atoms utilizes the MM parameters for the QM
> > atoms and, if your QM-MM border crosses bonds, any bond, angle or
> dihedral
> > term involving at least one QM atom is calculated using the regular MM
> > equations and parameters. So, the QM region does need MM parameters.
> >
> > Finally, I had the impression that your QM region is not a standard
> > residue. So, for any of those parts, you would need to create the
> > appropriate parameters for it.
> >
> > Gustavo Seabra
> > Professor Adjunto
> > Departamento de Química Fundamental
> > Universidade Federal de Pernambuco
> > Recife - PE - Brasil
> > +55-81-2126-7450 x5023
> > > ---------- Forwarded message ----------
> > > From: Gustavo Seabra
> > > Date: 14 November 2012 12:40
> > > Subject: Re: [AMBER] Problem with running QM/MM - MD simulation using
> > AM1/d-
> > > PhoT method in AMBER12
> > > To: AMBER Mailing List
> > >
> > >
> > > Hi Vojtech,
> > >
> > > Just to eliminate possibilities: Have you tried a different QM
> > hamiltonian?
> > > Does the same thing happen?
> > >
> > > Anyways, how did you parametrize the QM region?
> > >
> > > What I usually do when running QM/MM calculations is to first energy-
> > > minimise the system, then slowly heat it to the desired temperature,
> all
> > > with full-MM. Only after the system is reasonably thermalised at the
> > final
> > > temperature I turn on the QM region.
> > >
> > > Cheers,
> > > Gustavo Seabra
> > >
> > >
> > >
> > >
> > >
> > > On Wed Nov 14 2012 08:05:29 Vojtech said:
> > >> Dear all,
> > >>
> > >> I would like to mention here one strange issue that I have recently
> > >> encountered in AMBER12.
> > >>
> > >>
> > >>
> > >> I try to run QM/MM MD simulation using recently implemented SE -
> AM1/d-
> > > PhoT
> > >> method (abbreviation AM1/d) as the qm_theory for the QM part of
> > molecule.
> > >> The system is RNA-enzym (Ribozyme) and contains just NA-nucleobases,
> Na+
> > >> counter-ions and rectangular box of TIP3P molecules (around 30.000
> atoms
> > > in
> > >> total). The QM region contains just one phosphate with the ribose-ring
> > at
> > >> one site, methoxy group on the other and 2 nucleobases within H-bond
> > >> distance from both non-bridging oxygens of that phosphate (typically
> 49
> > >> atoms). The link atoms are thereby placed instead of C-C bonds (or C-N
> > > bonds
> > >> in case of nucleobases).
> > >>
> > >> Problem is that my QM region suddenly 'explodes' during initial
> > (typically
> > >> 50 ps-long) heating phase (constant volume, example of script at the
> > end)
> > > as
> > >> soon as the temperature of the system reaches ~250 K. Firstly hydrogen
> > > atoms
> > >> and them the other QM atoms just fly apart. I have also tried
> different
> > >> modifications of my input script, charmm parameters and charges
> > > (converting
> > >> starting pdb to charmm psf file, then backwards using chamber module),
> > > more
> > >> sequentional pre-heating with some equilibration steps and even larger
> > QM
> > >> regions. All modifications resulted in same error during pre-heating
> or
> > > just
> > >> at the beginning of subsequent equilibration. Notably, initial QM
> > >> minimalization fails (SCF convergency) for the system with largest QM
> > > region
> > >> (100 atoms). I have also tried to heat the system up to 200 K only
> and
> > >> thereby it 'survises' the pre-heating phase. But during the
> > equilibration
> > >> phase (constant pressure) same issue occurs (also suddenly, there is
> no
> > > hint
> > >> of something going wrong / any possible trend in both QM/Total
> Energies
> > ..
> > >>
> > >>
> > >>
> > >>
> > >> I am very confused and will much appreciate any reply or suggestion.
> > >>
> > >> Thanks in advance,
> > >> Vojtech
> >
> > _______________________________________________
> > 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 Wed Nov 14 2012 - 13:00:02 PST
Custom Search