Just to informally close out this discussion:
The problem I encountered seems to be a function of the choice of box shape,
in as much as I repeated the same type of simulations (literally the same
input files, just different parm7's) in a cube as opposed to a rhombic
dodecahedron and have not yet observed any such errors (108 independent 500
ps simulations so far).
Without fully understanding the PBC code in sander or sqm, I will offer the
following hypothesis (if anyone can offer further insight, please do).
Since non-orthorhombic and non-truncated octahedral PBCs in AMBER are
treated as triclinic shapes, the image wrapping does not necessarily supply
the same advantages as in the "familiar" polyhedron forms. More
specifically, the largest inscribe-able sphere may be smaller in triclinic
coordinates and the QM region may see different images than if the
simulation were done in "familiar" coordinates.
Brian
On Thu, Sep 8, 2011 at 5:34 PM, Qiantao Wang <qiantao.wang.gmail.com> wrote:
> Hi Ross,
>
> The switching function is a simple degree five polynomial. We switch
> both energy and force. By introducing the switching function (say S),
> we can break the QM-MM interaction into two different parts; one is
> the normal QM-MM (can we say (uv,ss)) interaction), and the other is
> the charge-charge (Q,q) interaction. The EQMMM = (uv,ss)*S +
> (Q,q)*(1-S). If the distance is shorter than the switching distance,
> we set S=1, so EQMMM = (uv,ss); and if the distance is beyond the
> switching region, we set S=0, so EQMMM = (Q,q), which is equivalent to
> what we have in QM/MM Ewald. Then we do the analytical gradient to
> obtain the force. The lower and upper boundaries can be specified by
> the user. Currently the lower and upper boundaries are arbitrarily
> taken to be qm_cut - 2 and qm_cut in default.
>
> Best,
> Qiantao
>
>
> 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
> >
>
> _______________________________________________
> 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 Tue Sep 13 2011 - 12:30:03 PDT
