Hi Brian,
There's nothing special about the Berendsen thermostat other than it
has been more thoroughly tested. I understand that the additions you
use have been well tested. Still, something is going wrong, and at the
moment we don't see what it is. Did you try printing the whole system
at every step until it breaks? Does at any point any atom seem out of
place or near the edge of the box, as Ross mentioned? Also, you
apparently have eliminated the TI options, so the only thing remaining
to test are the modifications.
Do you think you could send us the input files so we can try it out?
(if the files are too big, I'd recommend that you use dropbox to send
them: http://db.tt/h3uLcni)
Gustavo Seabra
Professor Adjunto
Departamento de Química Fundamental
Universidade Federal de Pernambuco
Fone: +55-81-2126-7417
On Thu, Sep 8, 2011 at 3:52 PM, Brian Radak <radak004.umn.edu> wrote:
> 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
Received on Thu Sep 08 2011 - 13:30:04 PDT
