Hi Ross,
To my knowledge the switching function is an atom based cubic polynomial on
the potential energy. The defaults suggested to me were to start at qm_cut
- 2 and end at qm_cut. That may not all be optimal, but it was shown to
enhance energy conservation, especially for systems with d-orbitals. If
Qiantao is watching the mailing list maybe he can clarify.
The latter case is almost certainly happening (or at least the former is
almost definitely not). The QM region (molecule really) does diffuse
considerably. However this never seemed to be a problem with, for example,
methanol (perhaps just fortuitous?).
There are bonds between QM atoms in my parm7 file and they are counted as a
single molecule there; I made great pains to make sure that that was true in
anticipation of the imaging issue. Was that not sufficient? On the other
hand, turning off imaging (iwrap = 0) does not seem to stop the error I am
seeing here (would you expect it to?). I am fairly certain (but am possibly
mistaken) that I have observed re-imagings of the QM region in other,
successful simulations.
Ross and Gustavo,
Attached are the input, output, and parm7 from a non-TI "standard" option
QM/MM simulation that crashed in 13 steps. The trajectory is too large to
send, so here is a link:
http://dl.dropbox.com/u/40984968/IS_NVT_0.00_1.nc
I have many others obviously.
Thanks for everyone's help on this,
Brian
On Thu, Sep 8, 2011 at 3:43 PM, Ross Walker <ross.rosswalker.co.uk> wrote:
> Hi Brian,
>
> I'd be interested to know the details of the switching function you are
> using for the QM region. Is this a 'force switching' function and if so how
> are you dealing with switching the force of the QMMM interaction.
>
> With regards to the problem you are seeing there are two likely causes of
> this that spring to mind.
>
> 1) You have solvent or some other mobile molecule in your QM region that is
> diffusing away from the rest of the QM region. The QM/MM cutoff is any MM
> atom which is within qmcut of ANY QM atom is included for ALL QM atoms.
> Thus
> if your QM region disperses as is the case if you have solvent or ions in
> your QM region then the effective QM region will grow in size until qmcut +
> max distance between any two QM atoms ends up bigger than half your box
> size.
>
> 2) Your QM region has diffused out of the central box. I am not sure how
> well this case was tested. If one of your QM atoms moves out the edge of
> the
> box it may not get imaged correctly within the QM region since the imaging
> is normally done as molecule based but the QM region is lacking bonds and
> thus each individual atom would be imaged. This could lead to the error you
> are seeing. The solution here is to reimage occasionally in ptraj centering
> on your QM region.
>
> In almost all cases showing the error you see it is because the QM region
> is
> diffusing in some way, the original code was designed for things like
> active
> sites of proteins etc which are essentially fixed over the typical
> simulation time scale. If your QM region can diffuse or dissolve as it were
> then you will likely have problems.
>
> All the best
> Ross
>
>
>
>
>
> > -----Original Message-----
> > From: Brian Radak [mailto:radak004.umn.edu]
> > Sent: Thursday, September 08, 2011 11:53 AM
> > To: AMBER Mailing List
> > Subject: Re: [AMBER] Fwd: false positive on too large QM region?
> >
> > That has been my feeling exactly whenever I encounter this problem (I'm
> > in
> > the dozens of systems now); nothing seems out of place at all. In fact
> > I am
> > wondering if a viable solution isn't just to comment out the SANDER
> > BOMB
> > call in the code and recompile (the Office Space approach to problem
> > solving, if you will).
> >
> > The "non-standard" options are well tested modifications to the sqm
> > module
> > done by our group and Case group. The first one just changes the hard
> > coded
> > AM1 parameters to the SRP I am using. Also, QM/MM MD won't always
> > conserve
> > energy (as well as one might like) unless qmmm_switch is used (we've
> > done ns
> > scale simulations to test that). They shouldn't affect the size extent
> > of
> > the system for any reason that I know of, but I'll give it a shot.
> >
> > I did run the system without TI options and got the identical crash (at
> > least in the lambda = 0 case the forces were all QM/MM anyway).
> >
> > Why would a Berendsen thermostat help here (other than hiding energy
> > conservation problems)? On the time scales of the error the dynamics
> > are
> > NVE, which would seem simpler to me.
> >
> > I feel like the problem must have *something* to do with imaging, but I
> > don't know when or why. Could there be something slightly naive in the
> > minimum image convention used in the QM region? Does it properly
> > account
> > for non-cubic boundary conditions? Perhaps this is my fault for not
> > knowing
> > more about the PBCs I chose...
> >
> > Brian
> >
> > On Thu, Sep 8, 2011 at 2:30 PM, Gustavo Seabra
> > <gustavo.seabra.gmail.com>wrote:
> >
> > > Hi Brian,
> > >
> > > I looked over your files, but nothing jumped out. However, I noticed
> > > that you are using some 'non-standard' options, such as
> > > "parameter_file = 'AM1D_HCNOPMGD.PAR'", and "qmmm_switch = 1", and
> > > some TI-related options. To pinpoint whats happening, I'd start by
> > > trimming the mdin file to the bare minimum default options and
> > > increase them progressively until something goes wrong. Also, you
> > > could try removing (as a test) the TI options, and changing the
> > > thermostat to Berendsen.
> > >
> > > Gustavo Seabra
> > > Professor Adjunto
> > > Departamento de Química Fundamental
> > > Universidade Federal de Pernambuco
> > > Fone: +55-81-2126-7417
> > >
> > >
> > >
> > > On Thu, Sep 8, 2011 at 11:42 AM, Brian Radak <radak004.umn.edu>
> > wrote:
> > > > Sorry, apparently the trajectory was too large. You probably
> > wanted the
> > > > output instead anyway.
> > > >
> > > > ---------- Forwarded message ----------
> > > > From: Brian Radak <radak004.umn.edu>
> > > > Date: Thu, Sep 8, 2011 at 10:32 AM
> > > > Subject: Re: [AMBER] false positive on too large QM region?
> > > > To: AMBER Mailing List <amber.ambermd.org>
> > > >
> > > >
> > > > Hi Gustavo,
> > > >
> > > > That was sensible advice; apparently I was missing 10 integration
> > steps
> > > > before the crash due to a large print step. Attached is a parm7
> > and a
> > > > NetCDF file with the 11 steps before the crash. Below is the input
> > file
> > > I
> > > > am using. Some probably very important clarifications:
> > > >
> > > > 1.) These are actually QM/MM-MM TI simulations; the error occurs in
> > ALL
> > > of
> > > > them in differing number of steps (they use the same starting
> > structure
> > > but
> > > > different random seeds)
> > > > 2.) The error occurs at lambda=0 (all QM) with and without coupling
> > to
> > > the
> > > > MM potential (i.e. I tried again with icfe = 0 in both inputs)
> > > > 3.) I have used this exact same protocol with other similar
> > molecules
> > > with
> > > > success, but have also had other systems where the error presents
> > itself
> > > > later in the simulation (say ~1ps in)
> > > > 4.) Suspecting an imaging problem, I tried iwrap = 0 and iwrap = 1,
> > the
> > > > problem seems to be identical in both
> > > >
> > > > I inspected the QM region pdb of one of the simulations and the
> > molecule
> > > is
> > > > intact.
> > > >
> > > > Brian
> > > >
> > > > QM/MM constant volume equilibration/production
> > > > &cntrl
> > > > nstlim = 500000, dt = 0.001,
> > > > ntpr = 1, ntave = 500000, ! frequency
> > of
> > > > printing energies and averages
> > > > ntwr = 500000, ntwx = 1, ioutfm = 1, ! frequency
> > of
> > > > restart and coord print in NetCDF format
> > > > ntb = 1, iwrap = 0, ! const
> > vol PBCs
> > > > with wrapping
> > > > ntc = 2, ntf = 2, tol = 1e-8, ! Shake
> > all
> > > > hydrogens, ignore all bonded terms
> > > > cut = 10, ! direct
> > space
> > > and
> > > > LJ MM cutoff
> > > > ntt = 2, temp0 = 310.0, tempi = 310.0, vrand = 2000, ! Andersen
> > > > thermostat
> > > > ig = 1315422706,
> > > > icfe = 1, klambda = 1,
> > > > clambda = 0.00,
> > > > ifqnt = 1,
> > > > noshakemask = ':1' ! Needs to
> > match
> > > > MM input for FEP-type calcs!
> > > > /
> > > > &ewald
> > > > ew_type = 0, order = 6,
> > > > nfft1 = 36, nfft2 = 36, nfft3 = 36
> > > > /
> > > > &qmmm
> > > > qmmask = ':1', qmcharge = 0,
> > > > qm_theory = 'AM1', parameter_file = 'AM1D_HCNOPMGD.PAR',
> > > > qm_ewald = 1, qmmm_switch = 1, qmcut = 10,
> > > > qmshake = 0,
> > > > writepdb = 1
> > > > /
> > > >
> > > >
> > > >
> > > > On Wed, Sep 7, 2011 at 8:19 PM, Gustavo Seabra
> > <gustavo.seabra.gmail.com
> > > >wrote:
> > > >
> > > >> Hi Brian,
> > > >>
> > > >> As you mentioned, yes, the coordinates are shifted for the QM
> > > calculation.
> > > >> In fact, the QM zone is put in the center of the box and the whole
> > > system is
> > > >> wrapped around it, but that only happens for the QM calculation,
> > without
> > > >> changing the actual coordinates.
> > > >>
> > > >> As to what is happening, it would be helpful to see your input and
> > > output
> > > >> files. Also, it may be a good idea to run the calculation printing
> > the
> > > >> coordinates at every step to follow your system. There is also an
> > option
> > > >> (can't recall now, but it is in the manual) that tells sander to
> > print
> > > the
> > > >> QM zone as a PDB file at every step. My guess is that somehow
> > sander is
> > > >> including in the QM zone some atom(s) that should not be there.
> > > >>
> > > >> Good luck,
> > > >> Gustavo.
> > > >> --
> > > >> Sent from my iPad.
> > > >>
> > > >> On 07/09/2011, at 18:37, Brian Radak <radak004.umn.edu> wrote:
> > > >>
> > > >> > I am running many QM/MM NVT simulations of small organic
> > molecules
> > > (QM)
> > > >> in
> > > >> > water (MM) using rhombic dodecahedral PBCs. The starting
> > structures
> > > come
> > > >> > from 1 ns NpT MM equilibrations with the volume reset to a
> > rounded
> > > value
> > > >> > (these are for TI calculations and I thought it sensible and not
> > > >> > unreasonable to set the volume constant amongst all lambda
> > values, the
> > > >> > change is usually less than 0.5%). Sometimes, but not always,
> > sander
> > > >> bombs
> > > >> > out with the following error, either after many steps of QM/MM
> > MD or
> > > >> > sometimes after no steps at all:
> > > >> >
> > > >> > ****************************************************
> > > >> > ERROR: QM region + cutoff larger than box dimension:
> > > >> > QM-MM Cutoff = 10.0000
> > > >> > Coord Lower Upper Size Radius of largest sphere
> > inside
> > > unit
> > > >> > cell
> > > >> > X -11.486 11.276 22.762 12.092
> > > >> > Y -12.143 12.058 24.200 12.092
> > > >> > Z -11.619 11.293 22.911 12.092
> > > >> > ****************************************************
> > > >> > SANDER BOMB in subroutine QM_CHECK_PERIODIC<qm_mm.f>
> > > >> > QM region + cutoff larger than box
> > > >> > cannot continue, need larger box.
> > > >> >
> > > >> > As a specific example this happens to isopropanol (generously ~6
> > A
> > > wide)
> > > >> in
> > > >> > a 34 A RHDO. I don't think I am wrong in concluding that there
> > is
> > > easily
> > > >> 12
> > > >> > A (qm_cut + skinnb) to the nearest RHDO face. In any event I
> > got no
> > > such
> > > >> > error during the MM simulation.
> > > >> >
> > > >> > I looked in the code and I don't fully understand how the table
> > of
> > > values
> > > >> is
> > > >> > arrived at. It would seem that the system is shifted by the
> > geometric
> > > >> > center of mass of the QM region (which is why "Lower" is
> > negative and
> > > >> > "Upper" is positive?), but I'm not exactly sure what the reason
> > for
> > > this
> > > >> is
> > > >> > (I'm sure it's a very good one, I just don't know what it is).
> > > Anyway,
> > > >> > somewhere in that shifting, the molecule seems to get split
> > across a
> > > >> > periodic wall and thus the "size" of the molecule is calculated
> > > >> > erroneously. There is also a comment in the code that something
> > bad
> > > >> (this?)
> > > >> > would happen if the QM region diffuses extensively.
> > > >> >
> > > >> > Can anyone explain what's going on here? Has anyone seen this
> > problem
> > > >> > before? I'm willing to give fixing this a try, but I think I
> > need a
> > > >> little
> > > >> > bit broader vision of the code as this extends across sqm, the
> > image
> > > >> code,
> > > >> > and the ewald code.
> > > >> >
> > > >> > Thanks!
> > > >> > Brian
> > > >> >
> > > >> > P.S. On a related note, 12.092 A seems too small for the largest
> > > sphere
> > > >> in
> > > >> > the "familiar" form of a RHDO. Is this a consequence of the
> > triclinic
> > > >> form
> > > >> > used in the simulation?
> > > >> >
> > > >> > --
> > > >> > ================================ 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 : Wright-
> > Rieman Hall
> > > >> 101
> > > >> > Graduate Program in Chemical Physics : 610 Taylor 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 : Wright-Rieman
> > Hall
> > > 101
> > > > Graduate Program in Chemical Physics : 610 Taylor 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.
> > > >
> > > >
> > > >
> > > >
> > > > --
> > > > ================================ 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 : Wright-Rieman
> > Hall
> > > 101
> > > > Graduate Program in Chemical Physics : 610 Taylor 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 : Wright-Rieman Hall
> > 101
> > Graduate Program in Chemical Physics : 610 Taylor 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 : Wright-Rieman Hall 101
Graduate Program in Chemical Physics : 610 Taylor 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 Sep 08 2011 - 14:00:03 PDT
