AMBER Archive (2007)

Subject: AMBER: box size bomb in QMMM calculation

From: M. L. Dodson (mldodson_at_houston.rr.com)
Date: Tue Jun 12 2007 - 14:46:00 CDT

Hi Ambers!

I have a problem, possibly conceptual.

I want to do QMMM and include in the QM region 5 waters close (within 5
A) to a key atom in the active site of an enzyme. I did 1.2ns of
completely standard classical MD, then determined to closest 5 waters to
the key atom (rectangular box of dimensions indicated below). But when
I used the last restart file of the classical MD run as the beginning
state of my QMMM run, I got the following error:

   .
   .
   .
| CHECK switch(x): max rel err = 0.2738E-14 at 2.422500
| CHECK d/dx switch(x): max rel err = 0.8314E-11 at 2.736960
  ---------------------------------------------------
| Local SIZE OF NONBOND LIST = 10918606
| TOTAL SIZE OF NONBOND LIST = 10918606
   ****************************************************
   ERROR: QM region + cutoff larger than box dimension:
   QM-MM Cutoff = 9.0000
    Coord Lower Upper Size Box(a,b,c)
      X 13.439 109.822 96.383 73.890
      Y 21.807 127.606 105.799 78.442
      Z 20.002 117.176 97.174 70.617
   ****************************************************
  SANDER BOMB in subroutine QM_CHECK_PERIODIC<qm_mm.f>
  QM region + cutoff larger than box
  cannot continue, need larger box.

See the enclosed .out file for the whole story (including the input).
Now the only thing I can see is that my method of determining the water
atoms to include in the QM region is wrong. Here is what I did:

1. imaged the restart file to the original cell, generating another
restart file (with ptraj):

center :1-${proteinl} mass origin
image origin
center :1-${protein1D} mass origin
image origin
center :1-${length} mass origin
image origin familiar

where proteinl, protein1D, and length are the number of residues in the
protein, protein + first DNA strand, and protein + both DNA strands

2. generated a pdb file from the restart file generated in 1 using ambpdb

3. generated a pdb file with the closest waters (using ptraj, and the
imaged restart file generated in step 1 as input):

closestwater 5 ":TT5\@C1'" oxygen noimage
(this is in a perl script, hence the \@ instead of just @)

4. Used grep to locate the Cartesian coordinates (oxygen) of these
closest waters in the pdb file generated with ambpdb in step 2.

5. Used the oxygen atom numbers identified in 4 + the succeeding two
numbers (Hydrogens) as the QM water atoms

This bombed as indicated above and in the attached output file.

Any comments would be appreciated.

My conceptual question is whether imaging renumbers the atoms? My
method assumes it does not. That is all I have been able to come up
with so far.

Thanks,
Bud Dodson

PS, this was run with a development version of the Amber10 sander
obtained from Ross. 

-- 
M. L. Dodson
Email:	mldodson-at-houston-dot-rr-dot-com
Phone:	eight_three_two-56_three-386_one

          -------------------------------------------------------
          Amber 9 SANDER 2006
          -------------------------------------------------------

| Run on 06/12/2007 at 14:06:20

File Assignments:
| MDIN: prod_md_1vas.in
| MDOUT: prod_md_1vas.out
|INPCRD: begin.restrt
| PARM: 1vas.prmtop
|RESTRT: prod_md_1vas.restrt
| REFC: refc
| MDVEL: prod_md_1vas.vel
| MDEN: mden
| MDCRD: prod_md_1vas.traj
|MDINFO: mdinfo
|INPDIP: inpdip
|RSTDIP: rstdip

 
 Here is the input file:
 
1vas, run steered QMMM, dropping snapshots along the way
 &cntrl
  ntx = 5, irest = 1, ntrx = 1, ntxo = 1,
  ntpr = 100, ntwx = 400, ntwv = 400,
  ntwe = 0,
  ntf = 2, ntb = 1,
  cut = 9.0, nsnb = 25,
  nstlim = 100000,
  t = 0.0, dt = 0.0005,
  ntt = 1, tautp = 10.0,
  temp0 = 300.0, tempi = 300.0,
  ig = 71277,
  tautp = 0.5,
  vlimit = 15.0,
  ntc = 2,
  jar = 1,
  ifqnt = 1,
  nmropt = 1,
 /
 &qmmm
  iqmatoms = 1,2,3,4,5,6,7,8,9,10,11,12,13,355,356,357,358,359,360,361,
             2465,2466,2467,2468,2469,2470,2471,2472,2473,
             2474,2475,2476,2477,2478,2479,2480,2481,2482,2483,2484,2485,
             2486,2487,2488,2489,2490,2491,2492,2493,
             2501,2502,2503,2504,2505,2506,2507,2508,
             2509,2510,2511,2512,2513,2514,2515,2516,2517,2518,2519,2520,
             2521,2522,2523,2524,2525,4171,4172,4173,25972,25973,25974,
             36247,36248,36249,36508,36509,36510,38641,38642,38643,
  qmcut = 9.0,
  qmtheory = 7,
  dftb_doscc = 1,
  dftb_disper = 0,
  qmcharge = 0,
  peptide_corr = 0,
  qmshake = 0,
  writepdb = 1,
  adjust_q = 2,
 /
 &wt type='DUMPFREQ', istep1 = 100,
 /
 &wt type='END',
 /
DISANG=dist.RST
DUMPAVE=dist_vs_t
LISTIN=POUT
LISTOUT=POUT
                                                                               
                                                                               
                                                                               
                                                                               

--------------------------------------------------------------------------------
   1. RESOURCE USE:
--------------------------------------------------------------------------------

 getting new box info from bottom of inpcrd
| INFO: Old style inpcrd file read

| peek_ewald_inpcrd: Box info found
|Largest sphere to fit in unit cell has radius = 35.308
| New format PARM file being parsed.
| Version = 1.000 Date = 06/16/06 Time = 10:14:05
 NATOM = 40152 NTYPES = 20 NBONH = 38477 MBONA = 1746
 NTHETH = 3289 MTHETA = 2473 NPHIH = 6394 MPHIA = 5478
 NHPARM = 0 NPARM = 0 NNB = 66561 NRES = 12528
 NBONA = 1746 NTHETA = 2473 NPHIA = 5478 NUMBND = 74
 NUMANG = 167 NPTRA = 70 NATYP = 40 NPHB = 1
 IFBOX = 1 NMXRS = 33 IFCAP = 0 NEXTRA = 0
 NCOPY = 0

| Memory Use Allocated
| Real 2012369
| Hollerith 253442
| Integer 1120837
| Max Pairs 17814104
| nblistReal 481824
| nblist Int 1668305
| Total 100957 kbytes
| Duplicated 0 dihedrals
| Duplicated 0 dihedrals

     BOX TYPE: RECTILINEAR

--------------------------------------------------------------------------------
   2. CONTROL DATA FOR THE RUN
--------------------------------------------------------------------------------

                                                                                

General flags:
     imin = 0, nmropt = 1

Nature and format of input:
     ntx = 5, irest = 1, ntrx = 1

Nature and format of output:
     ntxo = 1, ntpr = 100, ntrx = 1, ntwr = 500
     iwrap = 0, ntwx = 400, ntwv = 400, ntwe = 0
     ioutfm = 0, ntwprt = 0, idecomp = 0, rbornstat= 0

Potential function:
     ntf = 2, ntb = 1, igb = 0, nsnb = 25
     ipol = 0, gbsa = 0, iesp = 0
     dielc = 1.00000, cut = 9.00000, intdiel = 1.00000
     scnb = 2.00000, scee = 1.20000

Frozen or restrained atoms:
     ibelly = 0, ntr = 0

Molecular dynamics:
     nstlim = 100000, nscm = 1000, nrespa = 1
     t = 0.00000, dt = 0.00050, vlimit = 15.00000

Berendsen (weak-coupling) temperature regulation:
     temp0 = 300.00000, tempi = 300.00000, tautp = 0.50000

SHAKE:
     ntc = 2, jfastw = 0
     tol = 0.00001

NMR refinement options:
     iscale = 0, noeskp = 1, ipnlty = 1, mxsub = 1
     scalm = 100.00000, pencut = 0.10000, tausw = 0.10000

Ewald parameters:
     verbose = 0, ew_type = 0, nbflag = 1, use_pme = 1
     vdwmeth = 1, eedmeth = 1, netfrc = 1
     Box X = 73.890 Box Y = 78.442 Box Z = 70.617
     Alpha = 90.000 Beta = 90.000 Gamma = 90.000
     NFFT1 = 80 NFFT2 = 80 NFFT3 = 72
     Cutoff= 9.000 Tol =0.100E-04
     Ewald Coefficient = 0.30768
     Interpolation order = 4

QMMM options:
             ifqnt = True nquant = 89
              qmgb = 0 qmcharge = 0 adjust_q = 2
              spin = 1 qmcut = 9.0000 qmshake = 0
     lnk_atomic_no = 1 lnk_dis = 1.0900
          qmtheory = DFTB verbosity = 0
            qmqmdx = Analytical
      tight_p_conv = False (converge density to 0.05xSqrt[SCFCRT])
           scfconv = 0.100E-07 itrmax = 1000
      printcharges = False peptide_corr = False
    qmqmrij_incore = False qmmmrij_incore = False
  qmqm_erep_incore = False
       pseudo_diag = False
          qm_ewald = 1 qm_pme = True
            kmaxqx = 5 kmaxqy = 5 kmaxqz = 5 ksqmaxq = 27

--------------------------------------------------------------------------------
   3. ATOMIC COORDINATES AND VELOCITIES
--------------------------------------------------------------------------------

                                                                                
 begin time read from input coords = 1200.000 ps

           Begin reading energy term weight changes/NMR restraints
 WEIGHT CHANGES:
 DUMPFREQ 100 0 0.000000 0.000000 0 0
                         ** No weight changes given **

 RESTRAINTS:
 Requested file redirections:
  DISANG = dist.RST
  DUMPAVE = dist_vs_t
  LISTIN = POUT
  LISTOUT = POUT
 Restraints will be read from file: dist.RST
Here are comments from the DISANG input file:

 jar option running
******
 C1' ( 2472)-N1 ( 2477) NSTEP1= 0 NSTEP2=100000
R1 = -98.500 R2 = 1.500 R3 = 1.500 R4 = 101.500 RK2 =5000.000 RK3 = 5000.000
R1A= -98.500 R2A= 3.700 R3A= 3.700 R4A= 101.500 RK2A= 0.000 RK3A= 0.000
 Rcurr: 1.499 Rcurr-(R2+R3)/2: 0.001 MIN(Rcurr-R2,Rcurr-R3): 0.001
                       Number of restraints read = 1

                  Done reading weight changes/NMR restraints

 Number of triangulated 3-point waters found: 12341

     Sum of charges from parm topology file = -0.00000393
     Forcing neutrality...
QMMM: ADJUSTING CHARGES
QMMM: ----------------------------------------------------------------------
QMMM: adjust_q = 2
QMMM: Uniformally adjusting the charge of MM atoms to conserve total charge.
QMMM: qm_charge = 0
QMMM: QM atom RESP charge sum (inc MM link) = 3.183
QMMM: Adjusting each MM atom resp charge by = 0.000
QMMM: Sum of MM + QM region is now = 0.000
QMMM: ----------------------------------------------------------------------
| # of SOLUTE degrees of freedom (RNDFP): 82027.
| # of SOLVENT degrees of freedom (RNDFS): 0.
| NDFMIN = 82024. NUM_NOSHAKE = 0 CORRECTED RNDFP = 82024.
| TOTAL # of degrees of freedom (RNDF) = 82024.
 ---------------------------------------------------
 APPROXIMATING switch and d/dx switch using CUBIC SPLINE INTERPOLATION
 using 5000.0 points per unit in tabled values
 TESTING RELATIVE ERROR over r ranging from 0.0 to cutoff
| CHECK switch(x): max rel err = 0.2738E-14 at 2.422500
| CHECK d/dx switch(x): max rel err = 0.8314E-11 at 2.736960
 ---------------------------------------------------
| Local SIZE OF NONBOND LIST = 10918606
| TOTAL SIZE OF NONBOND LIST = 10918606
  ****************************************************
  ERROR: QM region + cutoff larger than box dimension:
  QM-MM Cutoff = 9.0000
   Coord Lower Upper Size Box(a,b,c)
     X 13.439 109.822 96.383 73.890
     Y 21.807 127.606 105.799 78.442
     Z 20.002 117.176 97.174 70.617
  ****************************************************
 SANDER BOMB in subroutine QM_CHECK_PERIODIC<qm_mm.f>
 QM region + cutoff larger than box
 cannot continue, need larger box.

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu