* This input file provides a general template for the production * of a solvated DNA box with Na+ counter ions. The script performs the * following tasks: * 1) Reads and processes the nucleic acid sequence from the input pdb file * 2) minimizes the DNA system with restraints to "regularize" the structure * OR reads the psf and coordinates for a DNA system * 3) Creates an excess of counter ions for establishment of optimal positioning * 4) Using the REPLICA algorithm and simple minimization, places the ** counter ions and eliminates the excess beyond those needed for ** charge neutrality. * 5) Optionally solvates the system * 6) writes a final coordinate set and PSF *** * Required files: top_all22_na.inp, par_all22_na.inp ** octa.img, tip125.crd, tip216.crd * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !* IMPORTANT USAGE NOTE: ONLY DOUBLE STRANDED DNA 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 !* DNA FROM PDB FILE. !**_______________________________________________________________________** !** WRITTEN BY C.L. BROOKS, III ** !** MAY, 2000 ** !** ** !** 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 !** !* 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 @?toppar eq 0 set toppar = "$toppar" if @?psf eq 0 set psf = false 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 ! Set-up DNA and minimize structure w/ restraints if psf ne true goto setupdna label fromsetupdna if psf eq true goto readpsfandcoor label fromreadpsfandcoor goto addcounterions label fromaddcounterions open unit 1 write form name coor/@pdbfile+NAmin.pdb write coor pdb unit 1 * Coordinates for @pdbfile after minimization w/ restraints * and addition of @TotNa Na+ counter ions to balance DNA charge * system "mkdir psf" open unit 1 write form name psf/@pdbfile+Namin.psf write psf card unit 1 * PSF for @pdbfile after minimization w/ restraints * and addition of @TotNa Na+ counter ions to balance DNA charge * if solvate eq true goto solvate label fromsolvate stop !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! label readpsfandcoor open unit 1 read form name psf/@pdbfile.psf read psf card unit 1 close unit 1 open unit 1 read form name coor/@pdbfile_minimized.pdb read coor pdb unit 1 close unit 1 goto fromreadpsfandcoor !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!PLACE COUNTER IONS AROUND DNA!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! label addcounterions ! What is the total charge to neutralize? scalar charge stats Calc TotNa = -?stot ! Use excess in search for minima in placement of ions. Calc Nreplicas = 3 * @TotNa prnlev 0 node 0 read sequ card unit 5 * Sequence for single "Fat" Sodium atom * 1 FSOD generate salt ! Make copies of segment and delete original ! Each replica segment does not interact with the others but "sees" the DNA replica s nreplicas @nreplicas select segid salt end delete atom select segid salt end ! Figure out dimenions of DNA for initial placement of Na+ ions coor orie select .not. segid s* end coor stat select .not. segid s* end Calc dx = ( ?Xmax - ?Xmin ) / ( @Nreplicas - 1 ) Set Dtheta = 15 Set XPos = ?Xmax Calc Ypos = ( ?Ymax - ?Ymin ) / 2 + 8 Set Theta = 0 set nrep = 1 label placereplica Scalar X set @Xpos select segid s@nrep end Scalar Y set @Ypos select segid s@nrep end Scalar Z set 0.0 select segid s@nrep end coor rotate xdir 1 ydir 0 zdir 0 phi @theta select segid s@nrep end Calc Xpos = @Xpos - @dx Calc Theta = @Theta + @DTheta incr nrep by 1 if nrep le @nreplicas goto placereplica open unit 1 write form name coor/@pdbfile_Na.pdb write coor pdb unit 1 * Coordinates for @pdbfile after addition of replica ions * coor copy compare ! Now minimize replicas in field of fixed DNA cons fix select .not. segid s* end mini sd nstep 1000 tolgrd 1.0e-5 cutnb 20 ctofnb 20 ctonnb 14 inbfrq 100 open unit 1 write form name coor/@pdbfile_NaM.pdb write coor pdb unit 1 * Coordinates for @pdbfile after addition of replica ions and minimization * ! Now find out how many ions moved coor diff compare scalar xcomp prod xcomp scalar ycomp prod ycomp scalar zcomp prod zcomp scalar xcomp store 1 scalar ycomp +store 1 scalar zcomp +store 1 scalar wmain recall 1 define moved select segid s* .and. property wmain .gt. 2 end set Nmoved = ?nsel if Nmoved le @TotNa goto dontdelete if Nmoved lt @nreplicas delete atom select segid s* .and. .not. moved end label dontdelete ! Now look for ions that me be too close to each other label findoverlaps define ions select segid s* end set nions = ?nsel set thision = ?selires label checkmind coor mind select ires @thision end select ( ions .and. .not. ires @thision ) end if ?mind gt 6 goto domoreions delete atom select ires @thision end goto findoverlaps label domoreions define ions select ions .and. .not. ires @thision end if ?nsel eq 1 goto ionsdone set thision = ?selires goto checkmind label ionsdone ! Check and see if we have excess or too few ions now scalar charge stat set Excession = ?stot if Excession eq 0 goto ionsbalance if Excession gt 0 goto removeions Calc Nadd = -@Excession read sequ FSOD @Nadd generate salt ! Figure out dimenions of DNA for initial placement of Na+ ions coor stat select .not. segid s* end Calc dx = ( ?Xmax - ?Xmin ) / ( @Nadd ) Set Dtheta = 360 / @Nadd Set XPos = ?Xmax Calc Ypos = ( ?Ymax - ?Ymin ) / 2 + 9 Set Theta = 0 set nion = 1 label placeion Scalar X set @Xpos select atom salt @nion sod end Scalar Y set @Ypos select atom salt @nion sod end Scalar Z set 0.0 select atom salt @nion sod end coor rotate xdir 1 ydir 0 zdir 0 phi @theta select atom salt @nion sod end Calc Xpos = @Xpos - @dx Calc Theta = @Theta + @DTheta incr nion by 1 if nion le @Nadd goto placeion goto ionsbalance label removeions system "date | awk '{split($4,s,":");print "* Title";print "*";print "set seed = " $3 s[1] s[2] s[3]}' > seed.str" stream seed.str system "rm seed.str" Calc seed = 2 * @seed + 1 random uniform iseed @seed scalar wcomp set 0 select ions end label loopions define ions select segid s* end if ?nsel eq @TotNa goto ionsbalance Calc thision = int ( ?rand * ?nsel ) + ?selires if thision gt ?nsel set thision = ?nsel delete atom select ires @thision end goto loopions label delthision label ionsbalance ! Now make the "fat" Na+ ions into real Na+ ions rename resname sod select resname fsod end open unit 1 write form name ions.crd write coor card unit 1 select segid s* end close unit 1 delete atom select segid s* end read sequ SOD @TotNa generate salt open unit 1 read form name ions.crd read coor card unit 1 select segid salt end close unit 1 system "rm ions.crd" prnlev 5 node 0 goto fromaddcounterions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! label readparams prnlev 0 node 0 wrnlev 0 node 0 open unit 1 read form name @TOPPAR/top_all22_na.inp read rtf card unit 1 close unit 1 read rtf card unit 5 append * This is an rtf patch to make a neutral C-terminus * appropriate for pH < ~3. 22 1 MASS 97 FSOD 22.989770 NA ! "Fat" Sodium Ion RESI FSOD 1.00 ! "Fat" sodium GROUP ATOM SOD FSOD 1.0 END open unit 1 read form name @TOPPAR/par_all22_na.inp read param card unit 1 close unit 1 read param card unit 5 append NONBOND FSOD 0.0 -0.01000 2.932 ! approximately Na + H2O radius w/ reduced well END prnlev 5 node 0 wrnlev 2 node 0 goto fromreadparams !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! label setupdna !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! 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 DNA SEQUENCE ! ! 2) BUILDS ALL MISSING ATOMS ! ! 3) WRITES THE PDB FILE ! ! CHECKS ARE DONE ON NUMBER OF MISSING ATOMS TO GUESS WHETHER FILE WAS ! ! PROCESSED CORRECTLY. EXAMINE OUTPUT FILE FOR POTENTIAL PROBLEMS! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE DOGENERATE label dogenerate prnlev 0 node 0 ! Copy original file for safe keeping then "fix" pdb file system "echo @PDBFILE.pdb > filename" system "cp coor/`awk '{print tolower($0)}' filename` pdb_orig" system "rm filename" ! 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 == "A") sub(substr($0,18,3),"ADE");' >> awk.file" system "echo 'if ($4 == "G") sub(substr($0,18,3),"GUA");' >> awk.file" system "echo 'if ($4 == "C") sub(substr($0,18,4),"CYT ");' >> awk.file" system "echo 'if ($4 == "T") sub(substr($0,18,3),"THY");' >> awk.file" system "echo 'if ($1 == "ATOM") sub("\\\*","z");' >> awk.file" system "echo '}' >> awk.file" system "echo '{if ( write == 0 ) {' >> awk.file" system "echo 'if ($1 == "ATOM") {print $0 "DNAA" >"dnaa.pdb";}' >> awk.file" system "echo '}}'>> awk.file" system "echo '{if ( write == 1 ) {' >> awk.file" system "echo 'if ($1 == "TER") {print "END" >"dnaa.pdb";}'>> awk.file" system "echo 'if ($1 == "ATOM") {print $0 "DNAB" >"dnab.pdb";}' >> awk.file" system "echo 'if ($1 == "END") {print "END" >"dnab.pdb";}'>> awk.file" system "echo '}}'>> awk.file" system "echo '{if ( write == 2 ) {' >> awk.file" system "echo 'if ($1 == "TER") {print "END" >"dnab.pdb";}'>> 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 -b1-72 pdb_orig | awk -f awk.file" system "mv DNAA.PDB dnaa.pdb" system "mv DNAB.PDB dnab.pdb" open unit 1 read form name dnaa.pdb read sequ pdb unit 1 close unit 1 generate dnaa setup first 5pho last 3pho ! Changed RNA to DNA w/ patch set nbase = 1 label dodnaa set pres = deo1 define thisresidue select ires @nbase end if ?selresn eq GUA set pres = deo2 if ?selresn eq ADE set pres = deo2 patch @pres dnaa @nbase setup warn increment nbase by 1 if nbase le ?nres goto dodnaa Calc nbase = ?nres + 1 set offset = ?nres rename atom AC1z select type AC1' end rename atom AC2z select type AC2' end rename atom AC3z select type AC3' end rename atom AC4z select type AC4' end rename atom AC5z select type AC5' end rename atom AH1z select type AH1' end rename atom AH2z select type AH2' end rename atom AH3z select type AH3' end rename atom AH4z select type AH4' end rename atom AH5z select type AH5' end rename atom AO2z select type AO2' end rename atom AO3z select type AO3' end rename atom AO4z select type AO4' end rename atom AO5z select type AO5' end rename atom C1z select type C1' end rename atom C2z select type C2' end rename atom C3z select type C3' end rename atom C4z select type C4' end rename atom C5z select type C5' end rename atom NC1z select type NC1' end rename atom NC2z select type NC2' end rename atom NC3z select type NC3' end rename atom NC4z select type NC4' end rename atom NC5z select type NC5' end rename atom NH1z select type NH1' end rename atom NH2z select type NH2' end rename atom NH3z select type NH3' end rename atom NH4z select type NH4' end rename atom NH5z select type NH5' end rename atom NO2z select type NO2' end rename atom NO3z select type NO3' end rename atom NO4z select type NO4' end rename atom NO5z select type NO5' end rename atom O2z select type O2' end rename atom O3z select type O3' end rename atom O4z select type O4' end rename atom O5z select type O5' end open unit 1 read form name dnaa.pdb read coor pdb unit 1 close unit 1 rename atom AC1' select type AC1z end rename atom AC2' select type AC2z end rename atom AC3' select type AC3z end rename atom AC4' select type AC4z end rename atom AC5' select type AC5z end rename atom AH1' select type AH1z end rename atom AH2' select type AH2z end rename atom AH3' select type AH3z end rename atom AH4' select type AH4z end rename atom AH5' select type AH5z end rename atom AO2' select type AO2z end rename atom AO3' select type AO3z end rename atom AO4' select type AO4z end rename atom AO5' select type AO5z end rename atom C1' select type C1z end rename atom C2' select type C2z end rename atom C3' select type C3z end rename atom C4' select type C4z end rename atom C5' select type C5z end rename atom NC1' select type NC1z end rename atom NC2' select type NC2z end rename atom NC3' select type NC3z end rename atom NC4' select type NC4z end rename atom NC5' select type NC5z end rename atom NH1' select type NH1z end rename atom NH2' select type NH2z end rename atom NH3' select type NH3z end rename atom NH4' select type NH4z end rename atom NH5' select type NH5z end rename atom NO2' select type NO2z end rename atom NO3' select type NO3z end rename atom NO4' select type NO4z end rename atom NO5' select type NO5z end rename atom O2' select type O2z end rename atom O3' select type O3z end rename atom O4' select type O4z end rename atom O5' select type O5z end open unit 1 read form name dnab.pdb read sequ pdb unit 1 close unit 1 generate dnab setup first 5pho last 3pho ! Changed RNA to DNA w/ patch label dodnab set pres = deo1 define thisresidue select ires @nbase end if ?selresn eq GUA set pres = deo2 if ?selresn eq ADE set pres = deo2 patch @pres dnab @nbase setup warn increment nbase by 1 if nbase le ?nres goto dodnab rename atom AC1z select segid dnab .and. type AC1' end rename atom AC2z select segid dnab .and. type AC2' end rename atom AC3z select segid dnab .and. type AC3' end rename atom AC4z select segid dnab .and. type AC4' end rename atom AC5z select segid dnab .and. type AC5' end rename atom AH1z select segid dnab .and. type AH1' end rename atom AH2z select segid dnab .and. type AH2' end rename atom AH3z select segid dnab .and. type AH3' end rename atom AH4z select segid dnab .and. type AH4' end rename atom AH5z select segid dnab .and. type AH5' end rename atom AO2z select segid dnab .and. type AO2' end rename atom AO3z select segid dnab .and. type AO3' end rename atom AO4z select segid dnab .and. type AO4' end rename atom AO5z select segid dnab .and. type AO5' end rename atom C1z select segid dnab .and. type C1' end rename atom C2z select segid dnab .and. type C2' end rename atom C3z select segid dnab .and. type C3' end rename atom C4z select segid dnab .and. type C4' end rename atom C5z select segid dnab .and. type C5' end rename atom NC1z select segid dnab .and. type NC1' end rename atom NC2z select segid dnab .and. type NC2' end rename atom NC3z select segid dnab .and. type NC3' end rename atom NC4z select segid dnab .and. type NC4' end rename atom NC5z select segid dnab .and. type NC5' end rename atom NH1z select segid dnab .and. type NH1' end rename atom NH2z select segid dnab .and. type NH2' end rename atom NH3z select segid dnab .and. type NH3' end rename atom NH4z select segid dnab .and. type NH4' end rename atom NH5z select segid dnab .and. type NH5' end rename atom NO2z select segid dnab .and. type NO2' end rename atom NO3z select segid dnab .and. type NO3' end rename atom NO4z select segid dnab .and. type NO4' end rename atom NO5z select segid dnab .and. type NO5' end rename atom O2z select segid dnab .and. type O2' end rename atom O3z select segid dnab .and. type O3' end rename atom O4z select segid dnab .and. type O4' end rename atom O5z select segid dnab .and. type O5' end open unit 1 read form name dnab.pdb read coor pdb unit 1 offset @offset close unit 1 rename atom AC1' select segid dnab .and. type AC1z end rename atom AC2' select segid dnab .and. type AC2z end rename atom AC3' select segid dnab .and. type AC3z end rename atom AC4' select segid dnab .and. type AC4z end rename atom AC5' select segid dnab .and. type AC5z end rename atom AH1' select segid dnab .and. type AH1z end rename atom AH2' select segid dnab .and. type AH2z end rename atom AH3' select segid dnab .and. type AH3z end rename atom AH4' select segid dnab .and. type AH4z end rename atom AH5' select segid dnab .and. type AH5z end rename atom AO2' select segid dnab .and. type AO2z end rename atom AO3' select segid dnab .and. type AO3z end rename atom AO4' select segid dnab .and. type AO4z end rename atom AO5' select segid dnab .and. type AO5z end rename atom C1' select segid dnab .and. type C1z end rename atom C2' select segid dnab .and. type C2z end rename atom C3' select segid dnab .and. type C3z end rename atom C4' select segid dnab .and. type C4z end rename atom C5' select segid dnab .and. type C5z end rename atom NC1' select segid dnab .and. type NC1z end rename atom NC2' select segid dnab .and. type NC2z end rename atom NC3' select segid dnab .and. type NC3z end rename atom NC4' select segid dnab .and. type NC4z end rename atom NC5' select segid dnab .and. type NC5z end rename atom NH1' select segid dnab .and. type NH1z end rename atom NH2' select segid dnab .and. type NH2z end rename atom NH3' select segid dnab .and. type NH3z end rename atom NH4' select segid dnab .and. type NH4z end rename atom NH5' select segid dnab .and. type NH5z end rename atom NO2' select segid dnab .and. type NO2z end rename atom NO3' select segid dnab .and. type NO3z end rename atom NO4' select segid dnab .and. type NO4z end rename atom NO5' select segid dnab .and. type NO5z end rename atom O2' select segid dnab .and. type O2z end rename atom O3' select segid dnab .and. type O3z end rename atom O4' select segid dnab .and. type O4z end rename atom O5' select segid dnab .and. type O5z end 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 prnlev 5 node 0 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 prnlev 0 node 0 faster on set k 30.0 label minimize cons harm force @k select .not. ( hydrogen .or. missing ) end minimize conj nstep 500 tolenr 0.001 - cutnb 14.0 rdie eps 2.0 vswitch shift coor orie rms select .not. hydrogen end decr k by 10.0 if k ge 10.0 goto minimize cons harm force 0 select .not. ( hydrogen .or. missing ) end ! Clean-up files from generate system "rm awk.file pdb_orig dnaa.pdb dnab.pdb" goto fromsetupdna !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SOLVATION COMPONENT: SOLVATES 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! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! For DNA don't include counter ions coor orie select .not. type sod end coor stat select .not. type sod end ! 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 coor/@pdbfile_solv-cube.crd if shape eq octahedral open unit 1 write form name coor/@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 coor/xtl_oct@nwat.crd if crystal ne true open unit 1 write form name coor/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 coor/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 "data/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 "data/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 !########################################################################### !###########################################################################