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