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
Received on Thu Sep 08 2011 - 12:00:03 PDT
