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

From: <vojtech.mlynsky.upol.cz>
Date: Thu, 15 Nov 2012 14:08:15 +0100 (CET)

Dear all,



.. thank you all for your suggestions and Mark for adding extra
clarifications.


Firstly, I have just completed another tests with qm_ewald=0 and larger
qmcut values but it resulted with same problem. I have also tested what
Gustavo suggested (slowly heating the system without qm_theory and then run
subsequent equilibration with qm region // also starting from my classical
MD restart file adding qmmm flag ..) but same issue occurred. That test
starting from my MD restart 'survived' for ~10 ps QM/MM dynamics, but
suddenly blows up after ..




And indeed, it looks like SCF convergence problem. I obtain these SCF
warnings in output files. As I have said, my system with the largest QM
region could not be minimized because SCF convergence fails. I have already
try to chance some parameters for that; namely turning pseudo_diag to 0,
using larger itrmax (more SCF steps than 1000), changing ndiis_attempts to
turn on earlier than on step 800 and also testing different diag_routines ..
. but I none of these helped for particular system (QM region consists of
100 atoms). It still appears that SCF settings are adjusted wisely (as it is
often stated in manuals).




The pure Single-Point sqm test with of QM region structure with all default
settings ends without any problem (pdb at the end of this mail if you want
to have a look). I should point out that it is the geometry of QM region
from starting structure. I will also try to check it for the geometry of QM
region just before disruption during pre-heating phase. It is still
interesting that it could be issue with SCF. It seems that I could
interpreted it wrongly as QM region blows first and then SCF convergence
(obviously) could not be achieved... I will also double-check that as well.
Anyway, thank you again (Andreas) for pointing on that.




With best regards,

Vojtech

 


==============================
Department of Physical Chemistry
Palacky University Olomouc
tr. 17. listopadu 12, 779 00 Olomouc
Czech Republic
==============================





QM region pdb:





REMARK

ATOM      1   C  QMR     1       1.730   1.989  -1.062

ATOM      2   H  QMR     1       1.759   1.561  -2.104

ATOM      3   O  QMR     1       3.087   1.852  -0.530

ATOM      4   C  QMR     1       3.224   0.589   0.193

ATOM      5   H  QMR     1       4.050   0.026  -0.326

ATOM      6   C  QMR     1       0.874   1.090  -0.155

ATOM      7   H  QMR     1       0.566   1.616   0.797

ATOM      8   C  QMR     1       1.855  -0.103   0.109

ATOM      9   H  QMR     1       1.612  -0.720   1.021

ATOM     10   O  QMR     1       1.940  -0.981  -1.033

ATOM     11   H  QMR     1       1.229  -1.701  -0.996

ATOM     12   O  QMR     1      -0.349   0.680  -0.836

ATOM     13   P  QMR     1      -1.162  -0.562  -0.093

ATOM     14   O  QMR     1      -1.436  -0.142   1.344

ATOM     15   O  QMR     1      -0.476  -1.856  -0.449

ATOM     16   O  QMR     1      -2.633  -0.544  -0.861

ATOM     17   C  QMR     1      -3.442  -1.753  -0.586

ATOM     18   H  QMR     1      -3.510  -1.882   0.530

ATOM     19   H  QMR     1      -2.935  -2.635  -1.073

ATOM     20   N  QMR     1       1.574   0.033   7.582

ATOM     21   C  QMR     1       2.101   1.339   7.698

ATOM     22   H  QMR     1       2.680   1.657   8.596

ATOM     23   N  QMR     1       1.852   2.099   6.625

ATOM     24   C  QMR     1       1.141   1.272   5.746

ATOM     25   C  QMR     1       0.650   1.514   4.432

ATOM     26   O  QMR     1       0.833   2.643   3.807

ATOM     27   N  QMR     1      -0.038   0.481   3.815

ATOM     28   H  QMR     1      -0.445   0.633   2.876

ATOM     29   C  QMR     1      -0.234  -0.769   4.477

ATOM     30   N  QMR     1      -0.995  -1.699   3.821

ATOM     31   H  QMR     1      -1.324  -1.511   2.870

ATOM     32   H  QMR     1      -1.073  -2.643   4.194

ATOM     33   N  QMR     1       0.259  -1.056   5.735

ATOM     34   C  QMR     1       0.948  -0.037   6.336

ATOM     35   N  QMR     1      -1.327  -0.874  -7.845

ATOM     36   C  QMR     1      -0.154  -1.651  -7.713

ATOM     37   H  QMR     1       0.430  -1.994  -8.596

ATOM     38   N  QMR     1       0.147  -1.939  -6.433

ATOM     39   C  QMR     1      -0.870  -1.354  -5.676

ATOM     40   C  QMR     1      -1.118  -1.365  -4.247

ATOM     41   N  QMR     1      -0.307  -1.963  -3.378

ATOM     42   H  QMR     1      -0.530  -2.086  -2.377

ATOM     43   H  QMR     1       0.581  -2.382  -3.693

ATOM     44   N  QMR     1      -2.298  -0.723  -3.853

ATOM     45   H  QMR     1      -2.627  -0.778  -2.875

ATOM     46   C  QMR     1      -3.117  -0.098  -4.779

ATOM     47   H  QMR     1      -4.050   0.387  -4.361

ATOM     48   N  QMR     1      -2.931  -0.020  -6.087

ATOM     49   C  QMR     1      -1.805  -0.676  -6.553

ATOM     50   H  QMR     1       1.481   3.049  -1.122

ATOM     51   H  QMR     1       3.490   0.787   1.231

ATOM     52   H  QMR     1      -4.436  -1.605  -1.009

ATOM     53   H  QMR     1       1.616  -0.764   8.325

ATOM     54   H  QMR     1      -1.800  -0.552  -8.773

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
Received on Thu Nov 15 2012 - 05:30:03 PST
Custom Search