* This is a general template script to create a solvated * simulation "box" using periodic boundaries with a few specific volumes. * The script will generate, minimize (or alternatively read in a psf * file for a solute system already generated) a solute system in a: * 1) cubic volume w/ images or crystal * 2) truncated octahedral box w/ images * 3) truncated octahedral box w/ crystal (reoriented from images) *** * Required files: top_all22_prot.inp, par_all22_prot.inp ** octa.img, tip125.crd, tip216.crd * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !* IMPORTANT USAGE NOTE: ONLY SINGLE CHAIN PROTEINS CAN BE PROCESSED !* WITH GENERATION/MINIMAZATION STAGE BECAUSE OF PSF GENERATION PROCEDURE !* USED, I.E. ASSUMPTION THAT PSF CAN BE GENERATED BY READING SEQUENCE OF !* PROTEIN FROM PDB FILE. !* FOR MORE GENERAL SOLVATION PROBLEMS, I.E. MULTI-CHAIN PROTEINS, !* SIMPLY JUST USE THE PSF=TRUE KEY AND THE SOLVATION PART. !**_______________________________________________________________________** !** WRITTEN BY C.L. BROOKS, III ** !** OCTOBER, 1999 ** !** ** !** COPYRIGHT C.L. BROOKS, III, THE SCRIPPS RESEARCH INSTITUTE ** !**_______________________________________________________________________** !******** !* The usage requires that several variables either be set in the !* input script or on the CHARMM command line at execution. !* The needed variables are: !** Variable Default Function !*__________________________________________________________________________ !* toppar ($toppar) ==> path to parameter and topology files !** !* pdbfile (...) ==> name of pdb coordinate file containing !** coordinates and sequence information !** (for psf=true key this is the name of the !** coordinate file in pdb format (pdbfile.pdb) !** and the pre-generated psf file (pdbfile.psf)) !* psf (false) ==> key specifying use of pre-generated psf and !** pdb coordinate file for use w/ multi-segid !** solute systems !* seqname (pdbfile) ==> name of protein sequence (segid) !** !* shape (cubic) ==> variable specifying whether truncated octahedral or !** cubic volume used (octahedral/cubic) !* crystal (false) ==> logical variable true/false specifying whether !** to use crystal orientation of simulation volume !* cutoff (9) ==> value of cutoff to be used in establishing how !** far images are to be from volume "edges", often set !** at or near nonbonded list cutoff value !** to use crystal orientation of simulation volume !* solvate (false) ==> logical flag which determines whether you stop after !** initial minimization and psf generation or proceed !** to solvation component. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if @?pdbfile eq 0 goto iamdead if @?seqname eq 0 set seqname = @pdbfile if @?psf eq 0 set psf = false if @?toppar eq 0 set toppar = "$toppar" if @?shape eq 0 set shape = cubic if @?crystal eq 0 set crystal = false if @?cutoff eq 0 set cutoff = 9 if @?solvate eq 0 set solvate = false ! Read the parameter and topology files: goto readparams label fromreadparams if psf ne true goto buildprotein label frombuildprotein if psf eq true goto readpsfandcoor label fromreadpsfandcoor if solvate eq true goto solvate label fromsolvate stop !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! label readpsfandcoor open unit 1 read form name @pdbfile.psf read psf card unit 1 close unit 1 open unit 1 read form name @pdbfile.pdb read coor pdb unit 1 close unit 1 goto fromreadpsfandcoor !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! label readparams prnlev 0 node 0 wrnlev 0 node 0 open unit 1 read form name "@TOPPAR/top_all22_prot.inp" read rtf card unit 1 close unit 1 open unit 1 read form name "@TOPPAR/par_all22_prot.inp" read param card unit 1 close unit 1 prnlev 5 node 0 wrnlev 2 node 0 goto fromreadparams !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! label buildprotein !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SOLUTE GENERATION COMPONENT: ASSUMES SINGLE CHAIN SOLUTE IN PDB ! ! FILE GIVEN BY KEY PDBFILE. IF MORE THAN ONE CHAIN IS PRESENT ! ! THE GENERATION ASPECTS WILL BE ERRONEOUS. ! ! THIS COMPONENT OF THE SCRIPT: ! ! 1) GENERATES PROTEIN SEQUENCE ! ! 2) BUILDS ALL MISSING ATOMS ! ! 3) MINIMIZES PROTEIN SUBJECT TO SUCCESSIVELY REDUCED HARMONIC RESTRAINTS ! ! 4) WRITES THE MINIMIZED PDB FILE ! ! CHECKS ARE DONE ON NUMBER OF MISSING ATOMS TO GUESS WHETHER FILE WAS ! ! PROCESSED CORRECTLY. EXAMINE OUTPUT FILE FOR POTENTIAL PROBLEMS! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Copy original file for safe keeping then "fix" pdb file system "echo @PDBFILE.pdb > filename" system "cp `awk '{print tolower($0)}' filename` pdb_orig" ! Make "standard" changes to accomodate pdb -> top_all22/28 atom/residue renames. system "echo 'BEGIN {write = 0}' > awk.file" system "echo '$1 == "TER" {write += 1;}'>> awk.file" system "echo '{' >> awk.file" system "echo 'if ($4 == "HIS") sub("HIS","HSD");' >> awk.file" system "echo 'if ( write == 0 ) {print $0 "@PDBFILE" ;}' >> awk.file" system "echo '}'>> awk.file" system "echo '{if ( write == 1 ) {' >> awk.file" system "echo 'if ($1 == "TER") {print "END" ;}'>> awk.file" system "echo '}}'>> awk.file" ! Beyond column 72 we should have the segid, some NMR files don't do this. Add ! segid after column 72 using cut and awk. system "cut -b 1-72 pdb_orig | awk -f awk.file > pdb_new" ! Clean-up temporary files system "rm awk.file" system "mv pdb_new `(awk '{print tolower($0)}' filename)`" open unit 1 read form name @PDBFILE.pdb read sequ pdb unit 1 close unit 1 generate @seqname setup first nter last cter ! Check for disulphide bonds in SSBOND records system "echo "* Disuplhides for @pdbfile" > disu.strm" system "echo "* " >> disu.strm" system "awk '$1 == "SSBOND" {print "patch disu @pdbfile",$4," @pdbfile",$6}' pdb_orig >> disu.strm" stream "disu.strm" system "rm disu.strm" rename atom cd1 select resname ile .and. type cd end rename atom oxt select type ot1 end rename atom o select ires ?nres .and. type ot2 end open unit 1 read form name @PDBFILE.pdb read coor pdb unit 1 close unit 1 rename atom cd select resname ile .and. type cd1 end rename atom ot1 select type oxt end rename atom ot2 select ires ?nres .and. type o end ! Now move original pdb file back and fnish clean-up system "mv pdb_orig `awk '{print tolower($0)}' filename`" system "rm filename" define nonH select .not. hydrogen end set nonHatoms = ?nsel define missing select .not. hydrogen .and. .not. initialized end set missingatoms = ?nsel Calc fraction = (@nonHatoms - @missingatoms) / @nonHatoms if fraction le 0.8 goto toomanywrong print coor select missing end ! Build missing atoms if just a few ic param ic build hbuild define missing select .not. initialized end if nsel ne 0 goto skipthis system "banner Some atoms not built" print coor select missing end stop label skipthis coor copy compare coor orie coor stat faster on set k 30.0 label minimize cons harm force @k select .not. hydrogen end minimize conj nstep 500 tolenr 0.001 - cutnb 14.0 rdie eps 2.0 vswitch shift coor orie rms select type ca end coor orie rms select .not. hydrogen end decr k by 10.0 if k ge 10.0 goto minimize open unit 1 write form name @pdbfile.psf write psf card unit 1 * PSF for @pdbfile sequence * open unit 1 write form name @pdbfile_minimized.pdb write coor pdb unit 1 * Coordinates for @pdbfile after minimization w/ restraints * goto frombuildprotein !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! label iamdead system "banner Must specify pdb filename" stop label toomanywrong system "banner Too many missing atoms" print coor select missing end stop label solvate !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! PROTEIN SOLVATION COMPONENT: SOLVATES PROTEIN SYSTEM IN SOLVENT VOLUME ! ! GIVEN BY KEYS SHAPE AND CRYSTAL AND CUTOFF ! ! THIS COMPONENT OF THE SCRIPT: ! ! 1) DETERMINES THE OPTIMAL SIZE FOR SOLVATION VOLUME ! ! 2) CONSTRUCTS A CUBIC VOLUME OF APPROPRIATE SIZE FROM BUILDING BLOCKS ! ! 3) CONSTRUCTS REQUESTED VOLUME FROM THIS CUBE AND WRITES FILE ! ! 4) DELETES OVERLAPPING WATER MOLECULES ! ! 5) WRITES CHARMM COORDINATE FILE FOR SOLVATED SYSTEM ! ! CHECKS ARE DONE ON NUMBER OF MISSING ATOMS TO GUESS WHETHER FILE WAS ! ! PROCESSED CORRECTLY. EXAMINE OUTPUT FILE FOR POTENTIAL PROBLEMS! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! coor orie coor stat ! Orient the solute system at the origin and use the coor stat command ! to figure out the maximum dimension for solvation. Set min = 0 Let min = mini ?xmin ?ymin Let min = mini @min ?zmin Let min = Abs @min Set max = 0 Let max = maxi ?xmax ?ymax Let max = maxi @max ?zmax Let max = maxi @max @min Calc maxD = @max ! For a cubic volume the cube size should be about MAXD + CUTOFF calc L = 2 * ( @maxD + @cutoff ) ! For a truncated octahedron volume the size of the cubic box ! used to construct the solvent volume should be around ! 4 * (CUTOFF + MAXD) / Sqrt(3). if shape eq octahedral - calc L = @L * 2 / Sqrt(3) ! @L is the edgelength of a cubic box from which you want to ! cut the truncated octahedron to solvate the solute. Figure ! out whether its optimal to construct this from a basic ! water box of 125 or 216 molecules. goto bldbox label donebldbox ! NOTE ASSUMPTION THAT ONLY SOLVENT SEGIDS BEGIN WITH W define solvent select segid w* end Calc size = @size * @nbox if shape ne cubic set size = @L Set 7 = @size set 8 = @7 set 9 = @7 set theta = 109.4712206344907 set xtl = @7 Calc xtl = sqrt (3) * @xtl / 2 if shape ne cubic goto makeshape label doneshape goto writesolvent label donewrite label deleteoverlapping define solute select .not. segid w* end ! Remove water overlaps with solute system delete atom sele .byres. ( segid w* .and. type OH2 .and. - (( solute .and. .not. hydrogen ) - .around. 2.8 )) end define solvent select segid w* .and. .not. hydrogen end set nwat ?nsel if shape eq cubic open unit 1 write form name @pdbfile_solv-cube.crd if shape eq octahedral open unit 1 write form name @pdbfile_solv-octa.crd write coor card unit 1 * X-ray structure of solute system @pdbfile solvated in a @shape volume constructed * from a cube of edgelength @size A. Volume contains @nwat water molecules. * If volume shape is cubic, use dimensions @size. * If volume shape is octahedral use basic dimensions of @7 x @8 x @9 and * a tetrahedral angle with octahedral symmetry for crystal facility * or use image facility with dimensions just noted. * Crystal facility logical flag = @crystal, use: * set theta = 109.4712206344907 * Crystal define octa @xtl @xtl @xtl @theta @theta @theta * set theta = 90.0 if shape eq octahedral set theta = 109.4712206344907 open unit 1 write form name @pdbfile_dimens.stream write title unit 1 ** Stream file to establish dimensions for ** dynamics simulation of @pdbfile using PBCs ** * set toppar = @toppar * set pdbfile = @pdbfile * set shape = @shape * set crystal = @crystal * set cutoff = @cutoff * set nwat = @nwat * set size = @7 * set theta = @theta goto fromsolvate !########################################################################### !########################################################################### label writesolvent define solvent select segid w* .and. .not. hydrogen end set nwat ?nsel if shape eq cubic goto writecube if crystal eq true open unit 1 write form name xtl_oct@nwat.crd if crystal ne true open unit 1 write form name img_oct@nwat.crd write coor card unit 1 select segid w* end * Coordinates for a truncated octahedron constructed from * a cubic volume of water molecules in a parallelpiped of * dimension @7 Ax@8 Ax@9 A on a side. There are @nwat molecules * in this file. * If key for crystal = @crystal use: * set theta = 109.4712206344907 * Crystal define octa @xtl @xtl @xtl @theta @theta @theta * goto donewrite label writecube open unit 1 write form name cube@nwat.crd write coor card unit 1 select segid w* end * Coordinates for a cubic volume of water molecules * in a parallelpiped of dimension @7 Ax@8 Ax@9 A on a side. * There are @nwat molecules in this file. * goto donewrite !########################################################################### !########################################################################### !########################################################################### !########################################################################### label makeshape prnlev 0 node 0 label setupimage open read card unit 1 name "/home/brooks/bin/INPUT/octa.img" read image card unit 1 close unit 1 Image byres select solvent end coor copy comp UPDATE inbfrq 1 imgfrq 1 ctonnb 3 ctofnb 3 cutnb 3 cutim 3 !Clean up the waters that have moved due to the UPDATE command coor diff comp scalar wmain = xcomp define xout select prop abs xcomp .gt. 0.1 .and. type oh2 end scalar wmain = ycomp define yout select prop abs ycomp .gt. 0.1 .and. type oh2 end scalar wmain = zcomp define zout select prop abs zcomp .gt. 0.1 .and. type oh2 end prnlev 5 node 0 delete atom select .byres. ( xout .or. yout .or. zout ) end if crystal ne true goto doneshape Calc onethird = 1.0000 / 3.0000 Calc twothird = 2.0000 / 3.0000 format (f15.8) coor rotate matrix @twothird -@onethird -@twothird ! image -> crystal -@onethird @twothird -@twothird ! image -> crystal @twothird @twothird @onethird ! image -> crystal goto doneshape !########################################################################### !########################################################################### !########################################################################### !########################################################################### label bldbox ! Choose from two different boxsizes as basic building unit. ! minimize remainder! ! Edglength of cubic volume of 125 molecules Calc nboxsmall = int (@L / 15.5575) calc remaindersmall = @L / 15.5575 - @nboxsmall if remaindersmall gt 0.0 Calc nboxsmall = @nboxsmall + 1 Calc remaindersmall = @nboxsmall * 15.5575 - @L ! Edglength of cubic volume of 216 molecules Calc nboxbig = int (@L / 18.856) calc remainderbig = @L / 18.856 - @nboxbig if remainderbig gt 0.0 Calc nboxbig = @nboxbig + 1 Calc remainderbig = @nboxbig * 18.856 - @L set size 15.5575 set nwat 125 Set nbox = @nboxsmall if remainderbig gt @remaindersmall goto notbig Set nbox = @nboxbig set size 18.856 set nwat 216 label notbig calc n (@nbox - 1) / 2 ! Build a big cube to surround solute prnlev 0 node 0 set totres 0 set xhigh @n set xlow -@n set yhigh @n set ylow -@n set zhigh @n set zlow -@n label start set seg ?nseg set xm @xlow set ym @ylow set zm @zlow label xmov set xdir @size mult xdir by @xm label ymov set ydir @size mult ydir by @ym label zmov set zdir @size mult zdir by @zm incr seg by 1 read sequ tip3 @nwat generate w@seg setup noangle nodihe open unit 1 read form name "/home/brooks/bin/INPUT/tip@NWAT.crd" !from test/data direc. if seg ne 1 goto append read coor card unit 1 goto skipappend label append read coor card unit 1 append label skipappend close unit 1 coor trans xdir @xdir ydir @ydir zdir @zdir select segid w@seg end incr zm by 1.0 if zm le @zhigh goto zmov set zm @zlow incr ym by 1.0 if ym le @yhigh goto ymov set zm @zlow set ym @ylow incr xm by 1.0 if xm le @xhigh goto xmov prnlev 5 node 0 coor stats goto donebldbox !########################################################################### !###########################################################################