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.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Wed Jun 13 2007 - 06:07:34 PDT
