* 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 !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 if psf ne true goto buildprotein 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 solvate 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 = 1}' > awk.file" system "echo '$1 == "ENDMDL" {write = 0;}'>> awk.file" system "echo '$1 == "END" {print $0}'>> awk.file" system "echo '{if ( write == 1 ) {' >> awk.file" system "echo 'if ($4 == "HIS" && $5 == "6") sub("HIS","HSD");' >> awk.file" system "echo 'if ($4 == "HIS" && $5 == "29") sub("HIS","HSE");' >> awk.file" system "echo 'if ($1 == "ATOM") {print $0}' >> awk.file" system "echo '}}'>> awk.file" system "awk -f awk.file pdb_orig > 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 prnlev 6 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" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Patch for changing HSD to HSE for this particular case !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! patch hs2 @pdbfile 29 rename resn HSE sele resi 29 end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Done Patch !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 open unit 1 write form name @pdbfile.psf write psf card unit 1 open unit 1 write form name @pdbfile_NOTmin.pdb write coor pdb unit 1 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 open unit 1 write form name @pdbfile_minimized.pdb write coor pdb unit 1 if solvate eq true goto solvate stop 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 * close unit 1 if shape eq cubic open unit 1 write form name @pdbfile_solv-cube.pdb if shape eq octahedral open unit 1 write form name @pdbfile_solv-octa.pdb write coor pdb unit 1 close unit 1 if shape eq cubic open unit 1 write form name @pdbfile_solv-cube.psf if shape eq octahedral open unit 1 write form name @pdbfile_solv-octa.psf write psf card unit 1 close unit 1 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 stop !########################################################################### !########################################################################### 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 !########################################################################### !###########################################################################