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

From: Brian Radak <radak004.umn.edu>
Date: Thu, 15 Nov 2012 12:30:47 -0500

Hopefully stupid questions: have you double checked the charge? Is the
multiplicity singlet? I can't think of any other reasons, from my
experience, that this would happen.

Brian

On Thu, Nov 15, 2012 at 12:01 PM, <vojtech.mlynsky.upol.cz> wrote:

> Dear all again,
>
>
>
> .. and sorry for prolonging ..
>
> I have finally run the sqm test for only QM-region system in geometry just
> before disruption during QM/MM-MD occurs .. and SCF convergence of Single-
> Point calculation fails (using default AMBER settings)! What is quite
> interesting for me that actual geometry of QM-region (attached at the end)
> is in fact reasonable and *very* similar to that starting one (I have
> posted
> and tested before..) ..
>
>
>
>
> With best regards,
>
> Vojtech.
>
>
> --
> ==============================
> Department of Physical Chemistry
> Palacky University Olomouc
> tr. 17. listopadu 12, 779 00 Olomouc
> Czech Republic
> ==============================
>
>
>
>
>
>
> QM-region pdb just before disruption ..
>
>
>
>
>
>
> REMARK
>
> ATOM 1 C QMR 1 1.749 1.728 -1.108
>
> ATOM 2 H QMR 1 1.972 1.196 -2.132
>
> ATOM 3 O QMR 1 3.047 1.621 -0.520
>
> ATOM 4 C QMR 1 3.176 0.344 0.224
>
> ATOM 5 H QMR 1 4.059 -0.221 -0.251
>
> ATOM 6 C QMR 1 0.863 0.817 -0.259
>
> ATOM 7 H QMR 1 0.577 1.377 0.689
>
> ATOM 8 C QMR 1 1.845 -0.324 0.035
>
> ATOM 9 H QMR 1 1.528 -0.842 0.998
>
> ATOM 10 O QMR 1 1.975 -1.219 -1.125
>
> ATOM 11 H QMR 1 1.152 -1.897 -1.173
>
> ATOM 12 O QMR 1 -0.367 0.438 -1.005
>
> ATOM 13 P QMR 1 -1.207 -0.805 -0.182
>
> ATOM 14 O QMR 1 -1.337 -0.315 1.244
>
> ATOM 15 O QMR 1 -0.496 -2.040 -0.682
>
> ATOM 16 O QMR 1 -2.685 -0.651 -0.775
>
> ATOM 17 C QMR 1 -3.743 -1.621 -0.402
>
> ATOM 18 H QMR 1 -3.869 -1.443 0.710
>
> ATOM 19 H QMR 1 -3.302 -2.687 -0.669
>
> ATOM 20 N QMR 1 1.552 0.211 7.611
>
> ATOM 21 C QMR 1 2.107 1.530 7.707
>
> ATOM 22 H QMR 1 2.704 1.906 8.624
>
> ATOM 23 N QMR 1 1.933 2.252 6.563
>
> ATOM 24 C QMR 1 1.083 1.433 5.777
>
> ATOM 25 C QMR 1 0.685 1.658 4.448
>
> ATOM 26 O QMR 1 0.999 2.687 3.707
>
> ATOM 27 N QMR 1 -0.059 0.638 3.820
>
> ATOM 28 H QMR 1 -0.341 0.690 2.831
>
> ATOM 29 C QMR 1 -0.209 -0.585 4.508
>
> ATOM 30 N QMR 1 -0.816 -1.556 3.827
>
> ATOM 31 H QMR 1 -1.020 -1.407 2.891
>
> ATOM 32 H QMR 1 -0.972 -2.547 4.132
>
> ATOM 33 N QMR 1 0.241 -0.891 5.782
>
> ATOM 34 C QMR 1 0.915 0.148 6.393
>
> ATOM 35 N QMR 1 -1.365 -0.970 -7.884
>
> ATOM 36 C QMR 1 -0.179 -1.678 -7.779
>
> ATOM 37 H QMR 1 0.407 -2.136 -8.624
>
> ATOM 38 N QMR 1 0.101 -1.969 -6.501
>
> ATOM 39 C QMR 1 -0.870 -1.303 -5.705
>
> ATOM 40 C QMR 1 -1.113 -1.439 -4.287
>
> ATOM 41 N QMR 1 -0.338 -2.081 -3.428
>
> ATOM 42 H QMR 1 -0.451 -2.052 -2.412
>
> ATOM 43 H QMR 1 0.649 -2.266 -3.723
>
> ATOM 44 N QMR 1 -2.324 -0.795 -3.854
>
> ATOM 45 H QMR 1 -2.652 -1.033 -2.901
>
> ATOM 46 C QMR 1 -3.206 -0.223 -4.796
>
> ATOM 47 H QMR 1 -4.059 0.230 -4.278
>
> ATOM 48 N QMR 1 -2.966 -0.107 -6.082
>
> ATOM 49 C QMR 1 -1.836 -0.727 -6.586
>
> ATOM 50 H QMR 1 1.465 2.778 -1.172
>
> ATOM 51 H QMR 1 3.417 0.507 1.274
>
> ATOM 52 H QMR 1 -4.676 -1.450 -0.938
>
> ATOM 53 H QMR 1 1.688 -0.587 8.341
>
> ATOM 54 H QMR 1 -1.900 -0.680 -8.788
>
> END
>
>
>
>
>
>
>
>
>
> ---------- Původní zpráva ----------
> Od: Andreas Goetz <agoetz.sdsc.edu>
> Datum: 14. 11. 2012
> Předmět: Re: [AMBER] Problem with running QM/MM - MD simulation using
> AM1/d-
> PhoT method in AMBER12
> "Dear all,
>
> It might be an SCF convergence problem in AMBER. Once the SCF fails to
> converge the forces will be wrong and it is not surprising that the system
> blows up.
>
> Vojtech and Mark, as an additional test, could you extract the QM region
> (including link atoms) and check whether you can converge a pure QM single
> point SCF with sqm? You can use 'writepdb = 1' in the &qmmm namelist, from
> which you can extract the coordinates and prepare an sqm input file.
> Thanks.
>
> All the best,
> Andy
>
> On Nov 14, 2012, at 12:46 PM, Marc van der Kamp wrote:
>
> > 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#
> (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
> (
> 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
> (
> 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/pubmed?term=York%20DM%5BAuthor%5D&cauthor=true&cauthor_uid=18345664
> )
> >
> >
> > http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2655239/
> (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
> (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
> (http://lists.ambermd.org/mailman/listinfo/amber)
> >>
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
> (http://lists.ambermd.org/mailman/listinfo/amber)
>
> --
> Dr. Andreas W. Goetz
> Assistant Project Scientist
> San Diego Supercomputer Center
> Tel : +1-858-822-4771
> Email: agoetz.sdsc.edu
> Web : www.awgoetz.de(http://www.awgoetz.de)
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
> (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 Thu Nov 15 2012 - 10:00:02 PST
Custom Search