Previous Thread
Next Thread
Print Thread
#26204 01/13/11 11:36 AM
Joined: Aug 2004
Posts: 192
E
Edo Offline OP
Forum Member
OP Offline
Forum Member
E
Joined: Aug 2004
Posts: 192
Dear All,
I am running all-atom simulations of a protein in explicit solvent and periodic boundary conditions using a parallel version of Charmm (c35). Unfortunately, the simulations die after several hundred thousand integration time-steps due to SHAKE errors (see below).

One potential correlation I see is that these simulations tend to die soon after one of the protein residues leaves the primary box; that is, the simulation runs fine while all the protein atoms are inside the primary box, but die a few thousand integration steps after one of the residues leave it. This would seem to suggest I've set up the image calculation wrong, but I'm not sure how that could be as the water molecules wrap fine during the simulations. Does anyone see anything wrong with my CHARMM script, or have a suggestion?

Please also note that I am using a non-standard dye molecule parameter and topology file.

Thanks!

Here is my input charmm script, and below this a sample output:
prnlev 3 node 0
set bmlev 0
bomlev @bmlev


!! Set up temperature run
set param ../setup/par_all27_prot_lipid.inp
set top ../setup/robs_dyes/top_all27_prot_lipid.inp
!set psf pass_this_variable
!set strt pass_this_variable
!set out pass_this_variable
!set rand pass_this_variable
!set dist 61.2 ! dye restraint distance

!! Set up periodic boundaries
set segname A ! segment name used in wrapping images
!set XSIZ 74.0 ! in angstroms
set YSIZ @XSIZ ! assumes cubic box
set ZSIZ @XSIZ
calc pcutim 10 !cutim
calc pcutnb 10 !cutnb

! read parameter and topology files
open unit 10 read form name @top
read rtf unit 10 card
close unit 10
open unit 1 read card name "../setup/dyes/alexa_top.inp"
read rtf card unit 1 append

open unit 10 read form name @param
read param unit 10 card
close unit 10
open unit 1 read card name "../setup/dyes/alexa_prm.inp"
read param card unit 1 append

! Get psf
read psf card name @psf

open read unit 3 card name @strt
read coor unit 3 card
close unit 3

! setup Periodic Boundary Conditions
set 9 @XSIZ
open unit 1 read form name ../setup/tips.img
read image card unit 1
close unit 1

crystal defi cubi @xsiz @ysiz @zsiz 90. 90. 90.
crystal build cutoff @pcutim nope 0
image byres xcen 0.0 ycen 0.0 zcen 0.0 sele all end

! SETUP NONBOND, this time we will be using electrostatics with PME, since
! we've set up a crystal structure.
nbond inbfrq -1 imgfrq -1 atom vatom cdie eps 1.0 -
elec ewald pmew fftx 64 ffty 64 fftz 64 kappa .34 spline order 6 -
vdw vswitch cutnb 16 cutim 16 ctofnb 12. ctonnb 10. wmin 1.0

cons fix sele .not. resname TIP3 end
mini sd nstep 50

cons fix sele none end
cons fix sele segid A .and. resid 4:79 end
mini sd nstep 50

cons fix sele none end

MMFP
GEO sphere RCM distance -
harmonic symmetric force 10.0 droff @dist -
select resname A488 .and. (type C1 .or. type C2 .or. -
type C3 .or. type C4 .or. type C5 .or. type O6) end -
select resname A594 .and. (type C1 .or. type C2 .or. -
type C3 .or. type C4 .or. type C5 .or. type O6) end
END

shake bond param
open unit 12 write form name prod/@OUT_@dist.res
open unit 21 write unform name prod/@OUT_@dist.dcd
open unit 22 write form name prod/@OUT_@dist.mon

dyna cpt leap strt time 0.002 nstep 25000000 iseed @rand -
pcons pref 1.0 pgamma 0 pmass 0 hoover reft 310 tmass 50000 -
iprfrq 10000 ieqfrq 0 ntrfrq 10000 -
iunwri 12 iuncrd 21 iunvel -1 kunit 22 -
ihtfrq 0 teminc 0 nprint 5000 nsavc 5000 nsavv 0 ihbfrq 0 -
inbfrq -1 imgfrq -1 isvfrq 5000 -
iasors 0 iasvel 1 iscvel 0 ichecw 0 bycb - !use bycube image nonbond list generator
echeck 10000

open unit 1 write form name prod/@OUT_@dist.cor
write coor card unit 1
* Cubic box
* ubq
*

stop

========================================================
========================================================
Excerpt from LOG file:
========================================================

Processing passed argument "xsiz=74.576"
Parameter: xsiz <- "74.576"
Processing passed argument "rand=63731"
Parameter: rand <- "63731"
Processing passed argument "psf=../setup/init/1ubq_dyes_water_1-neutralized.psf"
Parameter: psf <- "../setup/init/1ubq_dyes_water_1-neutralized.psf"
Processing passed argument "strt=../setup/equil/1ubq_1_57.66.cor"
Parameter: strt <- "../setup/equil/1ubq_1_57.66.cor"
Processing passed argument "out=1ubq_1"
Parameter: out <- "1ubq_1"
Processing passed argument "dist=57.66"
Parameter: dist <- "57.66"
1
Chemistry at HARvard Macromolecular Mechanics
Copyright(c) 1984-2001 President and Fellows of Harvard College
All Rights Reserved

Maximum number of ATOMS: 360720, and RESidues: 120240
Current HEAP size: 10240000, and STACK size: 10000000

Processing passed argument "xsiz=74.576"
Parameter: XSIZ <- "74.576"
Processing passed argument "rand=63731"
Parameter: RAND <- "63731"
Processing passed argument "psf=../setup/init/1ubq_dyes_water_1-neutralized.psf"
Parameter: PSF <- "../SETUP/INIT/1UBQ_DYES_WATER_1-NEUTRALIZED.PSF"
Processing passed argument "strt=../setup/equil/1ubq_1_57.66.cor"
Parameter: STRT <- "../SETUP/EQUIL/1UBQ_1_57.66.COR"
Processing passed argument "out=1ubq_1"
Parameter: OUT <- "1UBQ_1"
Processing passed argument "dist=57.66"
Parameter: DIST <- "57.66"
RDTITL> * CHARMM
RDTITL> *

CHARMM> !======== BEGIN: Set parameters ===========
CHARMM> prnlev 3 node 0

CHARMM> set bmlev 0
Parameter: BMLEV <- "0"

CHARMM> bomlev @bmlev
Parameter: BMLEV -> "0"

CHARMM>

CHARMM>

CHARMM> !! Set up temperature run
CHARMM> set param ../setup/par_all27_prot_lipid.inp
Parameter: PARAM <- "../SETUP/PAR_ALL27_PROT_LIPID.INP"

CHARMM> set top ../setup/dyes/top_all27_prot_lipid.inp
Parameter: TOP <- "../SETUP/DYES/TOP_ALL27_PROT_LIPID.INP"

CHARMM> !set psf init/1ubq_dyes_water_1-neutralized.psf
CHARMM> !set strt init/1ubq_dyes_water_1-neutralized.cor
CHARMM> !set out 1ubq_dyes_water
CHARMM> !set rand 4977913
CHARMM> !set dist 61.2 ! dye restraint distance
CHARMM>

CHARMM> !! Set up periodic boundaries
CHARMM> set segname A ! segment name used in wrapping images
Parameter: SEGNAME <- "A"

CHARMM> !set XSIZ 74.0 ! in angstroms
CHARMM> set YSIZ @XSIZ ! assumes cubic box
Parameter: XSIZ -> "74.576"
Parameter: YSIZ <- "74.576"

CHARMM> set ZSIZ @XSIZ
Parameter: XSIZ -> "74.576"
Parameter: ZSIZ <- "74.576"

CHARMM> calc pcutim 10 !cutim
Parameter: PCUTIM <- "10"

CHARMM> calc pcutnb 10 !cutnb
Parameter: PCUTNB <- "10"

CHARMM>

CHARMM> ! read parameter and topology files

CHARMM> open unit 10 read form name @top
Parameter: TOP -> "../SETUP/DYES/TOP_ALL27_PROT_LIPID.INP"
VOPEN> Attempting to open::../setup/dyes/top_all27_prot_lipid.inp::
OPNLGU> Unit 10 opened for READONLY access to ../setup/dyes/top_all27_prot_lipid.inp

CHARMM> read rtf unit 10 card
MAINIO> Residue topology file being read from unit 10.
TITLE> *>>>>>> COMBINED CHARMM ALL-HYDROGEN TOPOLOGY FILE FOR <<<<<<<<<
TITLE> *>>>>>>>>> CHARMM22 PROTEINS AND CHARMM27 LIPIDS <<<<<<<<<<
TITLE> *FROM
TITLE> *>>>>>>>>CHARMM22 ALL-HYDROGEN TOPOLOGY FILE FOR PROTEINS <<<<<<
TITLE> *>>>>>>>>>>>>>>>>>>>>> AUGUST 1999 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
TITLE> *>>>>>>> DIRECT COMMENTS TO ALEXANDER D. MACKERELL JR. <<<<<<<<<
TITLE> *>>>>>> 410-706-7442 OR EMAIL: ALEX,MMIRIS.AB.UMD.EDU <<<<<<<<<
TITLE> *AND
TITLE> * \\\\\\\ CHARMM27 ALL-HYDROGEN LIPID TOPOLOGY FILE ///////
TITLE> * \\\\\\\\\\\\\\\\\\ DEVELOPMENTAL /////////////////////////
TITLE> * ALEXANDER D. MACKERELL JR.
TITLE> * AUGUST 1999
TITLE> * ALL COMMENTS TO ADM JR. EMAIL: ALEX,MMIRIS.AB.UMD.EDU
TITLE> * TELEPHONE: 410-706-7442
TITLE> *

CHARMM> close unit 10
VCLOSE: Closing unit 10 with status "KEEP"

CHARMM> open unit 1 read card name "../setup/dyes/alexa_top.inp"
VOPEN> Attempting to open::../setup/dyes/alexa_top.inp::
OPNLGU> Unit 1 opened for READONLY access to ../setup/dyes/alexa_top.inp

CHARMM> read rtf card unit 1 append
MAINIO> Residue topology file being read from unit 1.
TITLE> * ADDITIONAL TOPOLOGY FOR RHODAMINE 6G
TITLE> * REF: VAIANA ET AL., J. COMP. CHEM., V24 (5), 632-639 (2003).
TITLE> *
WARNING from DECODF -- Zero length string being converted to 0.
WARNING from DECODF -- Zero length string being converted to 0.
WARNING from DECODF -- Zero length string being converted to 0.
WARNING from DECODF -- Zero length string being converted to 0.

CHARMM>

CHARMM> open unit 10 read form name @param
Parameter: PARAM -> "../SETUP/PAR_ALL27_PROT_LIPID.INP"
VOPEN> Attempting to open::../setup/par_all27_prot_lipid.inp::
OPNLGU> Unit 10 opened for READONLY access to ../setup/par_all27_prot_lipid.inp

CHARMM> read param unit 10 card

PARAMETER FILE BEING READ FROM UNIT 10
TITLE> *>>>>>> COMBINED CHARMM ALL-HYDROGEN PARAMETER FILE FOR <<<<<<<<<
TITLE> *>>>>>>>>> CHARMM22 PROTEINS AND CHARMM27 LIPIDS <<<<<<<<<<
TITLE> *FROM
TITLE> *>>>> CHARMM22 ALL-HYDROGEN PARAMETER FILE FOR PROTEINS <<<<<<<<<<
TITLE> *>>>>>>>>>>>>>>>>>>>>>> AUGUST 1999 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
TITLE> *>>>>>>> DIRECT COMMENTS TO ALEXANDER D. MACKERELL JR. <<<<<<<<<
TITLE> *>>>>>> 410-706-7442 OR EMAIL: ALEX,MMIRIS.AB.UMD.EDU <<<<<<<<<
TITLE> *AND
TITLE> * \\\\\\\ CHARMM27 ALL-HYDROGEN LIPID PARAMETER FILE ///////
TITLE> * \\\\\\\\\\\\\\\\\\ DEVELOPMENTAL /////////////////////////
TITLE> * ALEXANDER D. MACKERELL JR.
TITLE> * AUGUST 1999
TITLE> * ALL COMMENTS TO ADM JR. EMAIL:ALEX,MMIRIS.AB.UMD.EDU
TITLE> * TELEPHONE: 410-706-7442
TITLE> *
CHARMM> close unit 10
VCLOSE: Closing unit 10 with status "KEEP"

CHARMM> open unit 1 read card name "../setup/dyes/alexa_prm.inp"
OPNLGU> Unit already open. The old file will be closed first.
VCLOSE: Closing unit 1 with status "KEEP"
VOPEN> Attempting to open::../setup/dyes/alexa_prm.inp::
OPNLGU> Unit 1 opened for READONLY access to ../setup/dyes/alexa_prm.inp

CHARMM> read param card unit 1 append

PARAMETER FILE BEING READ FROM UNIT 1
TITLE> * ADDITIONAL PARAMETERS FOR RHODAMINE 6G
TITLE> * REF: VAIANA ET AL., J. COMP. CHEM., V24 (5), 632-639 (2003).
TITLE> *
PARRDR> ERROR: Repeated torsion periodicity: INDEX 321 CODE54395571 PERIODICITY 2 HP -CA -CA1 -OS
PARMIO> NONBOND, HBOND lists and IMAGE atoms cleared.

CHARMM>

CHARMM> ! Get psf

CHARMM> read psf card name @psf
Parameter: PSF -> "../SETUP/INIT/1UBQ_DYES_WATER_1-NEUTRALIZED.PSF"
VOPEN> Attempting to open::../setup/init/1ubq_dyes_water_1-neutralized.psf::
MAINIO> Protein structure file being read from unit 90.
TITLE> * 1CBN IN RHDO
TITLE> *
PSFSUM> PSF modified: NONBOND lists and IMAGE atoms cleared.
PSFSUM> Summary of the structure file counters :
Number of segments = 69 Number of residues = 13359
Number of atoms = 41240 Number of groups = 13688
Number of bonds = 41232 Number of angles = 15945
Number of dihedrals = 633 Number of impropers = 216
Number of cross-terms = 0
Number of HB acceptors = 13372 Number of HB donors = 152
Number of NB exclusions = 0 Total charge = 0.00000
VCLOSE: Closing unit 90 with status "KEEP"

CHARMM>

CHARMM> open read unit 3 card name @strt
Parameter: STRT -> "../SETUP/EQUIL/1UBQ_1_57.66.COR"
VOPEN> Attempting to open::../setup/equil/1ubq_1_57.66.cor::
OPNLGU> Unit 3 opened for READONLY access to ../setup/equil/1ubq_1_57.66.cor

CHARMM> read coor unit 3 card
SPATIAL COORDINATES BEING READ FROM UNIT 3
TITLE> * CUBIC BOX SIZE IS 74.576
TITLE> * UBQ
TITLE> *

CHARMM> close unit 3
VCLOSE: Closing unit 3 with status "KEEP"

CHARMM>

CHARMM> ! setup Periodic Boundary Conditions

CHARMM> set 9 @XSIZ
Parameter: XSIZ -> "74.576"
Parameter: 9 <- "74.576"

CHARMM> open unit 1 read form name ../setup/tips.img
OPNLGU> Unit already open. The old file will be closed first.
VCLOSE: Closing unit 1 with status "KEEP"
VOPEN> Attempting to open::../setup/tips.img::
OPNLGU> Unit 1 opened for READONLY access to ../setup/tips.img

CHARMM> read image card unit 1
TITLE> * IMAGE FILE FOR QUBIC TRANSFORMATION
TITLE> * BOX SIZE IS 74.576 ANGSTROMS ON A SIDE
TITLE> *

READING IMAGE FILE FROM UNIT 1


26 TRANSFORMATIONS HAVE BEEN READ

THERE ARE NO ROTATIONS FOR THIS TRANSFORMATION SET

CHARMM> close unit 1
VCLOSE: Closing unit 1 with status "KEEP"

CHARMM>

CHARMM> crystal defi cubi @xsiz @ysiz @zsiz 90. 90. 90.
Parameter: XSIZ -> "74.576"
Parameter: YSIZ -> "74.576"
Parameter: ZSIZ -> "74.576"
Crystal Parameters : Crystal Type = CUBI
A = 74.57600 B = 74.57600 C = 74.57600
Alpha = 90.00000 Beta = 90.00000 Gamma = 90.00000

CHARMM> crystal build cutoff @pcutim nope 0
Parameter: PCUTIM -> "10"

Range of Grid Search for Transformation 1 :
Lattice Vector A -2 TO 2
Lattice Vector B -2 TO 2
Lattice Vector C -2 TO 2


The number of transformations generated = 26


Number Symop A B C Distance

1 1 -1 -1 -1 4.3779
2 1 -1 0 -1 2.5800
3 1 -1 1 -1 4.3249
4 1 0 -1 -1 3.0833
5 1 0 0 -1 2.0565
6 1 0 1 -1 5.0243
7 1 -1 -1 0 2.4579
8 1 -1 0 0 2.3881
9 1 -1 1 0 4.1122
10 1 0 -1 0 2.9636
11 1 0 1 0 2.9636
12 1 -1 -1 1 7.6285
13 1 -1 0 1 5.7357
14 1 -1 1 1 5.7793
15 1 0 -1 1 5.0243
16 1 0 0 1 2.0565
17 1 0 1 1 3.0833
18 1 1 1 1 4.3779
19 1 1 0 1 2.5800
20 1 1 -1 1 4.3249
21 1 1 1 0 2.4579
22 1 1 0 0 2.3881
23 1 1 -1 0 4.1122
24 1 1 1 -1 7.6285
25 1 1 0 -1 5.7357
26 1 1 -1 -1 5.7793
THERE ARE NO ROTATIONS FOR THIS TRANSFORMATION SET
26 Transformations have been processed.


CHARMM> image byres xcen 0.0 ycen 0.0 zcen 0.0 sele all end
SELRPN> 41240 atoms have been selected out of 41240
IMAGE CENTERING ON FOR SOME ATOMS

.
.
.
.
.
DYNA> 365000 730.00000-108577.60454 25458.87899-134036.48353 310.63807
DYNA PROP> 16.37241-108519.96905 25627.44946 57.63549 15920.08962
DYNA INTERN> 0.00000 753.78265 93.12836 111.82497 57.83073
DYNA EXTERN> 13065.19899-120886.70083 0.00000 0.00000 0.00000
DYNA IMAGES> 307.45315 -8073.41164 0.00000 0.00000 0.00000
DYNA EWALD> 1278.51073-892069.63213 871325.50219 0.00000 0.00000
DYNA MMFP> 0.02930 0.00000 0.00000 0.00000 0.00000
DYNA PRESS> 7032.36558 -17645.75866 -1162.59197 -92.54486 414760.37337
DYNA XTLE> -153768.10957 -74.01043 -45189.90060 0.00000
---------- --------- --------- --------- --------- ---------
Crystal Parameters : Crystal Type = CUBI
DYNA A = 74.57600 B = 74.57600 C = 74.57600
DYNA Alpha = 90.00000 Beta = 90.00000 Gamma = 90.00000
DYNA PIXX = -19.99 PIYY = -48.50 PIZZ = -209.15
DYNA PIXY = -13.02 PIXZ = 46.72 PIYZ = -49.19
DYNA Gradient Norm = 282.89392

WRIDYN: RESTart file was written at step 365000
Using Image CUBE search
Using Image CUBE search
Using Image CUBE search
Using Image CUBE search
Using Image CUBE search
Using Image CUBE search
Using Image CUBE search
Using Image CUBE search
Using Image CUBE search
Using Image CUBE search
** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE
NITER= 0 LL= 1471 I= 1449 J= 1452
TOLER= 2.3409000000 DIFF= -32.8950602770 RRPR= -1.4882244243
X(I) -46.1839340705 -1.9950233006 -4.9728667772
XREF(I) -46.1782928205 -1.9868437177 -4.9781598609
X(J) -40.6814047937 -1.0985169582 -7.0111033223
XREF(J) -46.6044272086 -0.5337749513 -4.7592977057
*** LEVEL 1 WARNING *** BOMLEV IS 0
.
.
.
.
*** LEVEL 1 WARNING *** BOMLEV IS 0
***** ERROR IN SHAKEA ***** COORDINATE RESETTING WAS NOT ACCOMPLISHED IN 500 ITERATIONS

*** LEVEL -2 WARNING *** BOMLEV IS 0
BOMLEV HAS BEEN SATISFIED. TERMINATING.



/---------
/
/
/
! XXXX XXXX !
! XXXX XXXX !
! XXX XXX !
! X !
-- XXX /--
! ! XXX ! !
! ! ! !
! I I I I I !
! I I I I !
/
-- --
---/
XXX XXX
XXXX XXXX
XXXXX XXXXX
XXX XXX
XXX XXX
XXXXX
XXX XXX
XXX XXX
XXX XXX
XXXXX XXXXX
XXXX XXXX
XXX XXX


Execution terminated due to the detection of a fatal error.

ABNORMAL TERMINATION
MAXIMUM STACK SPACE USED IS 3591444
STACK CURRENTLY IN USE IS 123720
MOST SEVERE WARNING WAS AT LEVEL -2
HEAP PRINTOUT- HEAP SIZE 10240000
SPACE CURRENTLY IN USE IS 9696312
MAXIMUM SPACE USED IS 11928470
FREE LIST
PRINHP> ADDRESS: 1 LENGTH: 4049126 NEXT: 6826541
PRINHP> ADDRESS: 6826541 LENGTH: 1966784 NEXT: 5863960675889
PRINHP> ADDRESS: 5863960675889 LENGTH: 34 NEXT: 0

$$$$$ JOB ACCOUNTING INFORMATION $$$$$
ELAPSED TIME: 16.35 HOURS
CPU TIME: 16.30 HOURS
application called MPI_Abort(MPI_COMM_WORLD, 12) - process 0rank 0 in job 1 p1970_34506 caused collective abort of all ranks
exit status of rank 0: killed by signal 9
/usr/local/mvapich2-pgi/bin/mpdallexit
exit ( 0 )

Edo #26206 01/13/11 12:13 PM
Joined: Sep 2003
Posts: 4,863
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,863
Likes: 10
You cannot image proteins by residue; the unit by which imaging is done cannot be bonded to anything else:
IMAG BYSEG SELE my_protein END
IMAG BYRES SELE RESN TIP3 END

I also do not think that you should read images AND use the crystal code at the same time.

A cubic box is inefficient - try RHDO instead, and you may want to consider LOOKup (nbonds.doc) as well for speed.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
lennart #26209 01/13/11 02:13 PM
Joined: Aug 2004
Posts: 192
E
Edo Offline OP
Forum Member
OP Offline
Forum Member
E
Joined: Aug 2004
Posts: 192
Thanks Prof. Nilson, I forgot about the image wrapping by segment for the protein. I see a 20% speedup when I use the lookup table ("look sele resn TIP3 end"), thanks for the suggestion.
Ed


Moderated by  lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u2 Page Time: 0.009s Queries: 20 (0.005s) Memory: 0.7595 MB (Peak: 0.8526 MB) Data Comp: Off Server Time: 2023-01-27 08:25:09 UTC
Valid HTML 5 and Valid CSS