* This is a general template script to carry-out an FEP simulation * to compute the free energy difference for a thermodynamic half-cycle. * The FEP simulation uses the TSM module of CHARMM. * NOTE THIS EXAMPLE RUNS THE SIMULATIONS IN VACUUM AS AN ILLUSTRATION *** * Required files: top_all22_prot.inp, par_all22_prot.inp * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !* 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. !**_______________________________________________________________________** !** 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 !* resnum (33) ==> residue number of mutation site !* lambdastep (0.1) ==> stepsize for FEP perturbations !* dynamics (false) ==> logical determining whether FEP !** simulation to be performed !* analysis (false) ==> logical determining whether analysis of FEP !** simulation to be performed !** if @?pdbfile eq 0 goto iamdead if @?psf eq 0 set psf = false if @?toppar eq 0 set toppar = "$toppar" if @?resnum eq 0 set resnum = 33 if @?lambdastep eq 0 set lambdastep = 0.1 if @?dynamics eq 0 set dynamics = false if @?analysis eq 0 set analysis = false if dynamics eq true goto skipcheck if analysis eq true goto doanalysis label skipcheck ! Read the parameter and topology files: goto readparams label fromreadparams ! Here we read in the specific patch to make a hybrid residue read rtf card unit 5 append * Residue patch to make hybrid tyr/phe residue * PRES Y2F 0.0 GROUP ATOM CD2 CA -0.115 ATOM HD2 HP 0.115 GROUP ATOM CE1 CA -0.115 ATOM HE1 HP 0.115 GROUP !New atoms for Phe part of hybrid ATOM CZ1 CA -0.115 cz oh hh !Note we explicitly exclude all atoms of "reactant" ATOM HZ HP 0.115 cz oh hh !from those of product. This is required in TSM! bond cz1 ce2 bond cz1 ce1 bond cz1 hz angle cd1 ce1 cz1 !It is also necessary to specify all angles, dihdrals angle cd2 ce2 cz1 !and impropers (this system has none) associated with angle ce1 cz1 hz !hybrid residue patch specifically. One can also angle ce2 cz1 hz !use autogenerate, but then its necessary to delete angle he1 ce1 cz1 !extra ones that autogenerate makes, this is simpler. angle he2 ce2 cz1 dihedral cg cd1 ce1 cz1 dihedral cg cd2 ce2 cz1 dihedral cd1 ce1 cz1 ce2 dihedral ce1 cz1 ce2 cd2 dihedral cd1 ce1 cz1 hz dihedral cd2 ce2 cz1 hz dihedral he1 ce1 cz1 hz dihedral he2 ce2 cz1 hz dihedral hd2 cd2 ce2 cz1 dihedral hd1 cd1 ce1 cz1 ! ! Take ICs from top_all22_prot to permit building of missing atoms if necessary. ! IC CG CD1 CE1 CZ1 1.4059 120.6300 -0.1200 119.9300 1.4004 IC CZ1 CD1 *CE1 HE1 1.4004 119.9300 -179.6900 120.0100 1.0808 IC CZ1 CD2 *CE2 HE2 1.4000 119.9600 -179.9300 119.8700 1.0811 IC CE1 CE2 *CZ1 HZ 1.4004 119.9800 179.5100 119.9700 1.0807 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END if psf ne true goto buildprotein label frombuildprotein if psf eq true goto readpsfandcoor label fromreadpsfandcoor ! Add patch to complete hybrid residue patch y2f @pdbfile @resnum setup ! Now build missing atoms, place cz1 and cz at same position scalar x stats select atom @pdbfile @resnum cz end scalar x set ?stot select atom @pdbfile @resnum cz1 end scalar y stats select atom @pdbfile @resnum cz end scalar y set ?stot select atom @pdbfile @resnum cz1 end scalar z stats select atom @pdbfile @resnum cz end scalar z set ?stot select atom @pdbfile @resnum cz1 end hbuild !Check to make sure no atoms are missing coor print select .not. initialized end if ?nsel ne 0 goto stillmissing !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!SET-UP FOR TSM PERTURBATIONS!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Define reactant and product atoms define reactant select ( atom @pdbfile @resnum cz - .or. atom @pdbfile @resnum oh - .or. atom @pdbfile @resnum hh ) end define product select ( atom @pdbfile @resnum cz1 - .or. atom @pdbfile @resnum hz ) end ! Now set-up the TSM commands for FEP calculations and clean-up the structure a bit TSM Reactant select reactant end Product select product end Dont reactant bond angle Dont product bond angle Lambda 0.5 END mini sd nstep 200 TSM Clear !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!SET-UP FOR TSM PERTURBATIONS!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!! AND DYNAMICS ON SYSTEM !!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Write out initial structure to make sure it looks ok open unit 1 write form name coor/@pdbfile_y@{resnum}f_0.pdb write coor pdb unit 1 * Initial coordinates for protein atoms in hybrid structure * with Y@resnum => F hybrid * ! re-Define reactant and product atoms in case we solvated define reactant select ( atom @pdbfile @resnum cz - .or. atom @pdbfile @resnum oh - .or. atom @pdbfile @resnum hh ) end define product select ( atom @pdbfile @resnum cz1 - .or. atom @pdbfile @resnum hz ) end system "mkdir res crd pert" set nlambda = 1 set status start label collectpert Calc lambda = @lambdastep * @nlambda if lambda ge 1 goto donepert ! Now set-up the TSM commands for FEP calculations TSM Reactant select reactant end Product select product end Dont product bond angle Dont reactant bond angle Lambda @lambda Save Unit 10 Freq 10 END open unit 10 write form name pert/@pdbfile_@nlambda.prt open unit 11 write unform name crd/@pdbfile_@nlambda.crd open unit 12 write form name res/@pdbfile_@nlambda.res if nlambda eq 1 goto doingfirstlambda Calc n = @nlambda - 1 open unit 13 read form name res/@pdbfile_@n.res Set status = restart label doingfirstlambda prnlev 5 node 0 ! Run dynamics dynamics leap @status timestep 0.002 nstep 5000 nprint 1000 iprfrq 1000 - firstt 298 finalt 298 twindl -5.0 twindh 5.0 - ichecw 1 teminc 30.0 ihtfrq 20 ieqfrq 200 - iasors 1 iasvel 1 iscvel 0 isvfrq 1000 - iunwri 12 nsavc 10 nsavv 0 iunvel 0 - iunread 13 iuncrd 11 - !{* Nonbond options *} inbfrq -1 ilbfrq 0 - eps 1.0 cutnb 16 ctofnb 14 ctonnb 10 vswi fshift TSM Clear incr nlambda by 1 goto collectpert label donepert ! Write out final structure to make sure it looks ok open unit 1 write form name coor/@pdbfile_y@{resnum}f_1.pdb write coor pdb unit 1 * Final coordinates for protein atoms in hybrid structure * with Y@resnum => F hybrid * if analysis ne true stop !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!ANALYSIS OF TSM PERTURBATIONS!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! label doanalysis set nlambda = 1 label openfile Calc lambda = @lambdastep * @nlambda Calc filenumber = @nlambda + 9 if lambda ge 1 goto doneopenfile open unit @filenumber read form name pert/@pdbfile_@nlambda.prt incr nlambda by 1 goto openfile label doneopenfile Calc pstack = int ( 1 / @lambdastep ) * 2 tsm post pstack @pstack plot set nlambda = 1 label dowindow Calc lambda = @lambdastep * @nlambda if lambda ge 1 goto donewindow Calc filenumber = @nlambda + 9 Calc lambdap = @lambdastep * ( @nlambda - 1 ) if lambdap ne 0 Calc lambdap = @lambdap + @lambdastep / 2 !perturb lambda = @lambda => lambda' = @lambdap proc first @filenumber lambda @lambdap temp 298 deltat 10 bins 10 ctemp Calc lambdap = @lambdastep * ( @nlambda + 1 ) if lambdap gt 1 Set lambdap = 1 if lambdap ne 1 Calc lambdap = @lambdap - @lambdastep / 2 !perturb lambda = @lambda => lambda' = @lambdap proc first @filenumber lambda @lambdap temp 298 deltat 10 bins 10 ctemp incr nlambda by 1 goto dowindow label donewindow end 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 coor/`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 @pdbfile 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 "rm pdb_orig `awk '{print tolower($0)}' filename` 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 label stillmissing 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 cons harm force 0 select .not. hydrogen end system "mkdir psf" open unit 1 write form name psf/@pdbfile.psf write psf card unit 1 * PSF for @pdbfile sequence * open unit 1 write form name coor/@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