Topic Options
#12641 - 12/05/06 08:06 PM ** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE
accelrys Offline
Forum Member

Registered: 06/21/06
Posts: 21
Dear all,
When I used a bigger ECHECK such as 999999999.0 to avoid the Energy Exceed error,the problem:"** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE" displayed.How can I deal with this ?Thanks for any help!
My input file is below:
get aa definitions and parameters
read rtf card name $$$$/top_all22_prot.inp
read para card name $$$$/par_all22_prot.inp
! name of the system, used in filenames
set FILE pro
if @N .gt. 0 goto get_psf

bomblev -2
! Get psf and coordinates, from gen-prot.inp
read psf card name $$$$/pro_wat.psf
read coor card name $$$$/min_pro_wat.crd

! setup Periodic Boundary Conditions
! Here we have to use the crystal facility
! Variables XSIZ,YSIZ, and ZSIZ have already been defined
crystal define orthorombic 89.0 90.0 67.0 90.0 90.0 90.0
crystal build cutoff 14.0

! also specify that water molecules should be recentered as molecules (residues) if they
! get out of primary box, but that the protein molecule is a whole segment
image byresidue sele resn tip3 end xcen 0.0 ycen 0.0 zcen 0.0
image bysegment sele segi pro end xcen 0.0 ycen 0.0 zcen 0.0
! We need to keep the water molecules rigid, and also want to
! keep covalent X-H bonds fixed, so use SHAKE
shake bonH para ! take reference values for bond lengths from param file

set l 0
!!!!! INITIAL DYNAMICS
cons harm force 0.0 sele ALL end
open unit 12 write form name $$$$/dyn_pro@l.res
open unit 21 write unform name $$$$/dyn_pro@l.cor
open unit 22 write form name $$$$/dyn_pro@l.mon

! time to spend on heating may vary.
! 100 to 300 K in a few ps is probably typical
! and then the time for equilibration will have to be judged
! from case to case (5-50ps, or more)
! Use extended system pressure algorithm (see pressure.doc)
! with simple temperature check using veolicty scaling.
! For more precise temperature control Hoover's method may be usd
dynamics start CPT time 0.002 nstep 25000 iseed 31415 -
firstt 48.0 finalt 298.0 teminc 10.0 ihtfrq 500 -
PCONst pmass 400.0 pgamma 20.0 tbath 300.0 PREFerence 1.0 -
ieqfrq 2000 ichecw 1 twindl -5.0 twindh +5.0 iasors 0 -
nprint 500 iprfrq 500 -
atom cdie fshift vshift cutnb 14.0 ctofnb 12.0 -
inbfrq -1 ihbfrq 0 imgfrq -1 - ! -1 for inbfrq/imgfrq means heuristic update
iunwrit 12 iuncrd 21 kunit 22 nsavc 100 !ie update nonbond list when any atom has move
! > 0.5 * (CUTNB-CTOFNB) since last update
write coor card name $$$$/dyn_pro@l.crd
* pro in 89 x 90 x 67 waterbox
* After initial 50ps CPT dynamics
* RMSD vs minimized: ?RMS
*

! read initial coordinates again and compare
read coor card comp name $$$$/min_pro_wat.crd
coor orie rms sele segid pro end


set N 1
goto RUN_DYN
!!!!!!!!!!!!!!!! RESTART AND CONTINUE DYNAMICS
! have to get a few things in place before restarting dynamics
label GET_PSF
read psf card name $$$$/pro_wat.psf
read coor card name $$$$/min_pro_wat.crd

! setup PBC, this is not saved in the psf
! get variables XSIZ,YSIZ, and ZSIZ
stream $$$$/pro_box.str
! at restart the actual boxsize is taken from the restart file
! we just need to setup the transformations with a reasonably similar box
crystal define orthorombic @XSIZ @YSIZ @ZSIZ 90.0 90.0 90.0
crystal build cutoff 14.0
image byresidue sele resn tip3 end
image bysegment sele segi pro end

shake bonH param

label RUN_DYN

set M @N
decr M by 1
open unit 11 read form name $$$$/dyn_pro@M.res
open unit 12 write form name $$$$/dyn_pro@N.res
open unit 21 write unform name $$$$/dyn_pro@N.cor
open unit 22 write form name $$$$/dyn_pro@N.mon

dynamics restart timestep 0.002 nstep 50000 -
nprint 500 iprfrq 500 -
PCONst pmass 400.0 pgamma 20.0 tbath 300.0 PREFerence 1.0 -
ieqfrq 2000 ichecw 1 twindl -5.0 twindh +5.0 iasors 0 ECHECK 999999999.0 - !keep checking T, may not be needed
atom cdie fshift vshift cutnb 14.0 ctofnb 12.0 -
inbfrq -1 ihbfrq 0 imgfrq -1 - ! heuristic update
iunread 11 iunwrit 12 iuncrd 21 kunit 22 nsavc 100

write coor card name $$$$/dyn_pro@N.crd
* pro in 89 x 90 x 67 waterbox
* After initial 50ps + @N * 100ps CPT dynamics
* RMSD vs minimized: ?RMS
*
! read initial coordinates again and compare
read coor card comp name $$$$/min_pro_wat.crd
coor orie rms sele segid pro end

incr N by 1
if N le 100 goto RUN_DYN
stop

and some of the output file below:
SELECTED IMAGES ATOMS BEING CENTERED ABOUT 0.000000 0.000000 0.000000
UPDECI: Nonbond update at step 10280
UPDECI: Image update at step 10281

SELECTED IMAGES ATOMS BEING CENTERED ABOUT 0.000000 0.000000 0.000000
UPDECI: Nonbond update at step 10281
** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE
NITER= 0 LL= 3232 I= 6589 J= 6590
TOLER= 1.2343210000 DIFF=****************** RRPR= -1568.9700791794
X(I) 18.0150504659 -5.4850176563 5.6001591727
XREF(I) 17.9362001935 -4.9615621696 5.8995189580
X(J) 1251.5293421441 -1698.1147226588 -683.2654093789
XREF(J) 18.2176307928 -4.0793560047 6.5133838458
*** LEVEL 1 WARNING *** BOMLEV IS 0
** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE
NITER= 0 LL= 3233 I= 6589 J= 6591
TOLER= 1.2343210000 DIFF= -7321.9722864279 RRPR= -22.8421121133
X(I) 65.7676462986 -71.0111836219 -21.0676462189
XREF(I) 17.9362001935 -4.9615621696 5.8995189580
X(J) 18.0285045086 -4.8610249328 4.7845771669
XREF(J) 17.9816622948 -4.8427256383 4.7958287496
*** LEVEL 1 WARNING *** BOMLEV IS 0
** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE
NITER= 0 LL= 3234 I= 6589 J= 6592
TOLER= 1.2343210000 DIFF= -6668.5658340981 RRPR= -70.8628637134
X(I) 63.9198490732 -68.4507672611 -20.0670068095
XREF(I) 17.9362001935 -4.9615621696 5.8995189580
X(J) 18.6586021448 -5.7545066787 6.2084302100
XREF(J) 18.6486602509 -5.7311542721 6.2661927402
*** LEVEL 1 WARNING *** BOMLEV IS 0
** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE
NITER= 1 LL= 3232 I= 6589 J= 6590
TOLER= 1.2343210000 DIFF=****************** RRPR= -789.1104010089
X(I) 62.1679904897 -66.0240779915 -19.0500032358
XREF(I) 17.9362001935 -4.9615621696 5.8995189580
X(J) 682.5249487439 -917.3262510183 -365.5005179526
XREF(J) 18.2176307928 -4.0793560047 6.5133838458
*** LEVEL 1 WARNING *** BOMLEV IS 0
** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE
NITER= 1 LL= 3233 I= 6589 J= 6591
TOLER= 1.2343210000 DIFF= -6814.9362462743 RRPR= -22.4873673048
X(I) 86.1836286414 -98.9802086081 -32.4620075475
XREF(I) 17.9362001935 -4.9615621696 5.8995189580
X(J) 40.0462549798 -35.3701131367 -7.1387164301
XREF(J) 17.9816622948 -4.8427256383 4.7958287496
*** LEVEL 1 WARNING *** BOMLEV IS 0
** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE
NITER= 1 LL= 3234 I= 6589 J= 6592
TOLER= 1.2343210000 DIFF= -6490.7059551502 RRPR= -70.1854324042
X(I) 84.3978520485 -96.5181391284 -31.4818533483
XREF(I) 17.9362001935 -4.9615621696 5.8995189580
X(J) 39.5331789759 -34.6701463787 -5.9098534445
XREF(J) 18.6486602509 -5.7311542721 6.2661927402
...
*** LEVEL 1 WARNING *** BOMLEV IS 0
** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE
NITER= 18 LL= 3232 I= 6589 J= 6590
TOLER= 1.2343210000 DIFF= -0.0000000003 RRPR= -0.7907238613
X(I) 100.7999433561 -118.9984067630 -40.6209485712
XREF(I) 17.9362001935 -4.9615621696 5.8995189580
X(J) 101.4213741017 -119.8514572810 -40.9680071289
XREF(J) 18.2176307928 -4.0793560047 6.5133838458
*** LEVEL 1 WARNING *** BOMLEV IS 0
BOMLEV HAS BEEN SATISFIED. TERMINATING.

ABNORMAL TERMINATION
MAXIMUM STACK SPACE USED IS 3990934
STACK CURRENTLY IN USE IS 172110
MOST SEVERE WARNING WAS AT LEVEL -2

$$$$$ JOB ACCOUNTING INFORMATION $$$$$
ELAPSED TIME: 11.82 HOURS
CPU TIME: 11.52 HOURS

Top
#12642 - 12/05/06 08:48 PM Re: ** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE [Re: accelrys]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
Setting ECHECK to much larger than 9999. is usually wishful thinking.

The error is a symptom with many causes, including errors in model building and a parallel CPT bug in older versions. I've also noticed that most people who report the error seem to be using TIMESTEP 0.002; I generally use TIMESTEP 0.001 and have not seen this error in my own work for quite some time.

I should add that CPT PCONS using the ORTH lattice type is not recommended; the box can change shape, your macromolecule can rotate, and Bad Things Can Happen. Isotropic shapes (RHDO, OCTA, CUBIC) are best for a solvated large molecule or complex; RHDO is the closest to spherical.

Speaking of CPT, the HOOVER thermostat is the preferred method of T control for doing NPT dynamics with that code; using the simple 3-step Verlet (non-zero IEQFRQ, ICHECW 1) is crude and only intended for initial equilibration. (It also means you aren't really doing true NPT dynamics, either.)
_________________________
Rick Venable
computational chemist


Top
#12643 - 12/10/06 01:33 PM Re: ** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE [Re: rmv]
user28 Offline
Forum Member

Registered: 08/10/04
Posts: 63
A related question, Rick. If one is using a 1fs timestep, does one need to use shake at all? what is advisable?

Thanks.

Top
#12644 - 12/10/06 02:11 PM Re: ** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE [Re: user28]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
The TIP3P water model is rigid, so one should use SHAKE for that; we often use SHAKE BONH for all molecules for consistency.
_________________________
Rick Venable
computational chemist


Top
#12645 - 10/05/07 01:29 AM Re: ** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE [Re: rmv]
vbrahma Offline
Forum Member

Registered: 11/01/06
Posts: 11
Loc: India
I have 2 problems

firstly the * ERROR IN SHAKEA * DEVIATION IN SHAKE TOO LARGE

secondly ***** LEVEL -2 WARNING FROM <DYNAMC> *****
***** ENERGY CHANGE TOLERANCE EXCEEDED
******************************************
BOMLEV ( -2) IS REACHED - TERMINATING. WRNLEV IS 5

the Total charge = 6.00000, Is this the reason??
how can it be circumvented?

because apart from that according to me everything seems ok?
INBFRQ=-1 ,IHBFRQ=0, CTOFNB - CUTONNB difference is 2, minimization steps are 100 for SD and ABNR, SHAKE IS DEFAULT 0.1000E -09, TIMESTEP IS 1fs, I am using RHDO as the isotropic shape, also ECHECK IS 999.0
IS IT BECAUSE I AM USING ieqfrq 2000 ichecw 1 ????
or do i need to increase twindl -5.0 twindh +5.0???


Attachments
15788-run-cpt-md_old.out0.txt (465 downloads)

_________________________
Vijaya Brahma Research Scholar Institute of MIcrobial Technology India

Top
#12646 - 10/05/07 10:25 AM Re: ** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE [Re: vbrahma]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
A new post would be more appropriate, as the questions are not related to this message thread.

The large initial VDW and image VDW energies indicate severe packing problems, i.e. there are way too many bad atom overlaps from the model building.

Top
#12647 - 10/05/07 10:35 AM Re: ** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE [Re: vbrahma]
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4742
Loc: ~ 59N, 15E
It is difficult to read your long output file; it seems you are reading rtf,params etc several times, which looks strange.
Please post brief input file as well as relevant parts of the output file. As to the total charge, and how to neutralize it, you should consider how this works in reality, and also perhaps discuss with your supervisor.
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top

Moderator:  John Legato, lennart