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#> 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
--
Dr. Andreas W. Goetz
Assistant Project Scientist
San Diego Supercomputer Center
Tel : +1-858-822-4771
Email: agoetz.sdsc.edu
Web : www.awgoetz.de
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Nov 14 2012 - 13:30:02 PST
