#!/bin/bash ################################################################ # calc_inter_model.sh # # bash script for calculating interaction energies with AMBER8 # Stefan Henrich # 07.07.2005 calc_interaction.sh # 12.09.2005 calc_inter_model # 15.09.2005 # # # EML/AZ version ################################################################ environment="combine2go_linux" #environment settings #EML_linux #AZ_linux #combine2go_linux if [ "$environment" = "EML_linux" ]; then echo EML_linux COMBINEHOME="/home/henricsn/combine2go/programs/" sdfdir=$COMBINEHOME"sdf/" prgdir=$COMBINEHOME"programs/linux/" #directory for programs paradir=$COMBINEHOME"parameter/" #directory for parameter leaprc="leaprc.ff03" #in directory: /sw/mcm/app/amber/amber8_linux/dat/leap/cmd/ #ligdir=$COMBINEHOME"data/urokinase/pdb/ligand/" #directory with prep/parm files of ligand #directory=$COMBINEHOME"data/urokinase/prepared/" ligdir=$COMBINEHOME"data/urokinase/prepared/" elif [ "$environment" = "AZ_linux" ]; then echo combine2go_AZ COMBINEHOME="/home/tassfun/data/combine2go/" sdfdir=$COMBINEHOME"sdf/" prgdir=$COMBINEHOME"programs/linux/" #directory for programs paradir=$COMBINEHOME"parameter/" #directory for parameter leaprc="/home/kemifr/AMBER8/dat/leap/cmd/leaprc.ff02" #in directory: /sw/mcm/app/amber/amber8_linux/dat/leap/cmd/ #ligdir=$COMBINEHOME"data/urokinase/pdb/ligand/" #directory with prep/parm files of ligand #directory=$COMBINEHOME"data/urokinase/prepared/" ligdir=$COMBINEHOME"data/urokinase/prepared/" elif [ "$environment" = "combine2go_linux" ]; then echo combine2go_linux COMBINEHOME="/home/henricsn/combine2go/" sdfdir=$COMBINEHOME"sdf/" prgdir=$COMBINEHOME"programs/linux/" #directory for programs paradir=$COMBINEHOME"parameter/" #directory for parameter leaprc="leaprc.ff03" #in directory: /sw/mcm/app/amber/amber8_linux/dat/leap/cmd/ #ligdir=$COMBINEHOME"data/urokinase/pdb/ligand/" #directory with prep/parm files of ligand #directory=$COMBINEHOME"data/urokinase/prepared/" targetdir="/scratch/mcm/henricsn/tmp-aout/" ############## ### max. length of path and name is 80 characters !!!! sdf=$sdfdir"all_300805_3D.sdf" ligdir=$COMBINEHOME"data/urokinase/prepared/model/model5/" # ligdir=$COMBINEHOME"data/urokinase/pdb/gold/uPA-docking/temp19/min/" # ligdir="../" # pdbout=$COMBINEHOME"data/urokinase/anal/model/model5/" # pdbout=$COMBINEHOME"data/urokinase/pdb/gold/uPA-docking/temp19/min/anal/" pdbout="./" ############## else echo give environment settings fi ################################################################ ### !!!! fuer anal berechnung muss eigentlich .xyz und nicht .cor eingelesen werden, da .cor von tleap stammt !!! #getfile=$(find $ligdir -name "*.*.minH.xyz") getfile=$(find $ligdir -name "*.*.min.xyz") #getfile=$(find $ligdir -name "*.cor") prop_activity=uPA golpeoutdir="./" golpeout="output" golpeout="minH" golpedat=$golpeoutdir$golpeout".dat" golpenames=$golpeoutdir$golpeout".dat.VarNames" golpevartyp=$golpeoutdir$golpeout".dat.VarType" ### pymol script echo run /home/henricsn/Program/pymol/scripts/data2bfactor.py > color-vdw.pml echo run /home/henricsn/Program/pymol/scripts/color_b.py >> color-vdw.pml echo >> color-vdw.pml echo run /home/henricsn/Program/pymol/scripts/data2bfactor.py > color-elec.pml echo run /home/henricsn/Program/pymol/scripts/color_b.py >> color-elec.pml echo >> color-elec.pml echo run /home/henricsn/Program/pymol/scripts/data2bfactor.py > color.pml echo run /home/henricsn/Program/pymol/scripts/color_b.py >> color.pml echo >> color.pml ################################################################## ################################################################## # DO NOT CHANGE ANYTHING BELOW THIS LINE # UNLESS YOU KNOW WHAT YOU ARE DOING ################################################################## ### write golpe file ($golpedat) molnumber=1 ################################################################## ### prepare amber8 ################################################################## if [ "$environment" = "EML_linux" ]; then source /sw/mcm/app/prepare/prepare.bash prepare amber8 elif [ "$environment" = "AZ_linux" ]; then AMBERHOME=/home/kemifr/AMBER8// AMBEREXE=${AMBERHOME}/exe_fix elif [ "$environment" = "combine2go_linux" ]; then source /sw/mcm/app/prepare/prepare.bash prepare amber8 else echo give environment settings for amber fi ################################################################## FILES="$getfile" for file in $FILES do ### (./A14.0.cor) base=`basename $file ` #file without path (A14.0.cor) name=`expr "$base" : '\([0-9a-zA-Z][0-9a-zA-Z]*\)' ` #base name (A14) suffix="${base#*.}" #suffix (0.cor) counter=`expr "$suffix" : '\([0-9a-zA-Z][0-9a-zA-Z]*\)' ` #counter (0) echo $base $name $suffix $counter ###grep residue name and sequence number of ligand pdb file line=$(grep "C1" $ligdir$name"."$counter".minbres.pdb"|uniq -w6) residue=${line:17:3} lig_number=${line:20:6} # case $residue in # *[!0-9]*|"") activity=$(grepsdtag.linux -t residue $residue < $sdf | extract_prop_sdf -property_name $prop_activity); # echo FAILURE;; # *) id=$residue;activity=$(grepsdtag.linux -t ID $id < $sdf | extract_prop_sdf -property_name $prop_activity); # echo SUCCESS;; # esac # case $residue in ## *[!0-9]*|"") echo FAILURE;; ## *) echo SUCCESS;; # *[!0-9]*|"") activity=$(grepsdtag.linux -t residue $residue < $sdf | extract_prop_sdf -property_name $prop_activity); # id=$(grepsdtag.linux -t residue $residue < $sdf | extract_prop_sdf -property_name ID); # echo residue okay;; # *) activity=$(grepsdtag.linux -t residue $residue < $sdf | extract_prop_sdf -property_name $prop_activity); # echo FAILURE;; # esac echo $residue $name"."$counter seq number $lig_number ### coordinates, topology, output coordinates=$ligdir$name"."$counter.min.xyz # coordinates=$ligdir$name"."$counter.cor topology=$ligdir$name"."$counter.top analout=$targetdir$name"."$counter.outa ################################################################## ### anal module of amber ################################################################## while read line # For as many lines as the input file has... do neu="$line" neu="${neu//999/$lig_number}" echo "$neu" done < $paradir"anal.in" > $targetdir"anal.in" # reading template of minimize script and output into new file analin=$targetdir"anal.in" analin=${analin:=$prgdir"anal.in"} ################################################################## ### !!!! fuer anal berechnung muss eigentlich .xyz und nicht .cor eingelesen werden, da .cor von tleap stammt !!! if [ "$environment" = "EML_linux" ]; then echo anal -O -i $analin -o $targetdir$pdb$alt.outa -p $directory$pdb$alt.top -c $directory$pdb$alt.cor anal -O -i $analin -o $targetdir$pdb$alt.outa -p $directory$pdb$alt.top -c $directory$pdb$alt.cor elif [ "$environment" = "AZ_linux" ]; then echo $AMBEREXE/anal -O -i $analin -o $targetdir$pdb$alt.outa -p $directory$pdb$alt.top -c $directory$pdb$alt.cor $AMBEREXE/anal -O -i $analin -o $targetdir$pdb$alt.outa -p $directory$pdb$alt.top -c $directory$pdb$alt.cor elif [ "$environment" = "combine2go_linux" ]; then echo anal -O -i $analin -o $analout -p $topology -c $coordinates anal -O -i $analin -o $analout -p $topology -c $coordinates else echo tleap version? fi #tleap environment ################################################################## ### end of anal module of amber ################################################################## #exit 99 ################################################################## ### write golpe file ($golpedat) ################################################################## if [ $molnumber == 1 ]; then #title echo Golpe dat for COMBINE analysis, $(date +"%D %T") > $golpedat prot_components=$[lig_number-1] lig_components=1 #number of variables (x+y) #ligand*protein*(vdw+ele)+(sum_vdw+sum_ele)+(999999.999+999999.999)+activity+deltaG echo $[($lig_components*$prot_components)*2+2+2+2] >> $golpedat #number of objects echo $getfile|wc| awk '{ print $2 } ' >> $golpedat ################################################################## #write .dat.VarNames #Reste zaehlen echo $[($lig_components*$prot_components)*2+2+2+2] > $golpenames # echo $(awk '{if($1=="ATOM" && $3=="CA"){num++}}END{print (num*2)}' $file) > $golpenames #variable names as residue names and numbers awk '($1=="ATOM" && $3=="CA") {print "v"$4_$5}' $ligdir$name"."$counter".minbres.pdb" >> $golpenames awk '($1=="ATOM" && $3=="CA") {print "e"$4_$5}' $ligdir$name"."$counter".minbres.pdb" >> $golpenames echo "sum_vdw" >> $golpenames echo "sum_elec" >> $golpenames echo "dummy" >> $golpenames echo "dummy" >> $golpenames echo "activity" >> $golpenames echo "deltaG" >> $golpenames ################################################################## #write .dat.VarType echo 7 > $golpevartyp #number of blocks echo 1 >> $golpevartyp echo $prot_components >> $golpevartyp echo 2 >> $golpevartyp echo " vdw" >> $golpevartyp echo 1 >> $golpevartyp #x variable echo $[prot_components+1] >> $golpevartyp echo $[prot_components*2] >> $golpevartyp echo 2 >> $golpevartyp echo " elec" >> $golpevartyp echo 1 >> $golpevartyp #x variable echo $[prot_components*2+1] >> $golpevartyp echo $[prot_components*2+1] >> $golpevartyp echo 2 >> $golpevartyp echo " sum_vdw" >> $golpevartyp echo 3 >> $golpevartyp #not used echo $[prot_components*2+2] >> $golpevartyp echo $[prot_components*2+2] >> $golpevartyp echo 2 >> $golpevartyp echo " sum_elec" >> $golpevartyp echo 3 >> $golpevartyp #not used echo $[prot_components*2+3] >> $golpevartyp echo $[prot_components*2+4] >> $golpevartyp echo 2 >> $golpevartyp echo " dummy" >> $golpevartyp echo 3 >> $golpevartyp #not used echo $[prot_components*2+5] >> $golpevartyp echo $[prot_components*2+5] >> $golpevartyp echo 2 >> $golpevartyp echo " activity" >> $golpevartyp echo 3 >> $golpevartyp #y varible not used echo $[prot_components*2+6] >> $golpevartyp echo $[prot_components*2+6] >> $golpevartyp echo 2 >> $golpevartyp echo " deltaG" >> $golpevartyp echo 2 >> $golpevartyp #y varible ################################################################## fi #sequential number of molecule echo $molnumber >> $golpedat #name of molecule echo $residue"."$counter >> $golpedat #x-variables $prgdir"anal2golpe" $analout $prot_components $lig_components >> $golpedat echo anal2golpe $analout $prot_components $lig_components #y-variables activity=$(grepsdtag.linux -t residue $residue < $sdf | extract_prop_sdf -property_name $prop_activity) #dummy lines echo 999999.999 >> $golpedat echo 999999.999 >> $golpedat #activity echo $activity >> $golpedat #deltaG #kJ/mol echo "298*0.008314*l($activity*10^-6)" | bc -l >> $golpedat ################################################################## ### end of write golpe file ($golpedat) ################################################################## molnumber=$[molnumber+1] rm -f $renumb ################################################################## $prgdir"anal2golpe-vdw" $analout $prot_components $lig_components > $name"."$counter"-temp.vdw" $prgdir"anal2golpe-elec" $analout $prot_components $lig_components > $name"."$counter"-temp.elec" awk '($1=="ATOM" && $3=="CA") {print $5,$4}' $ligdir$name"."$counter".minbres.pdb" > $name"."$counter".seq" paste $name"."$counter".seq" $name"."$counter"-temp.vdw" > $name"."$counter".vdw" paste $name"."$counter".seq" $name"."$counter"-temp.elec" > $name"."$counter".elec" #pymol: #run /home/henricsn/Program/pymol/scripts/data2bfactor.py #run /home/henricsn/Program/pymol/scripts/color_b.py #data2b_res A04.0-vdw,/home/henricsn/combine2go/data/urokinase/anal/model/model5/seq echo load $ligdir$name"."$counter".minbres.pdb",$name"."$counter".vdw" >> color-vdw.pml echo data2b_res $name"."$counter".vdw",$pdbout$name"."$counter".vdw" >> color-vdw.pml #color_b(selection='model5', item='b', gradient='rwb', mode='hist',minimum='-1',maximum='1') echo color_b\(selection=\'\/$name"."$counter".vdw"\',item=\'b\',gradient=\'rwb\',mode=\'hist\',minimum=\'-1\',maximum=\'1\'\) >> color-vdw.pml echo save $pdbout$name"."$counter".vdw.pdb", \($name"."$counter".vdw"\) >> color-vdw.pml echo load $ligdir$name"."$counter".minbres.pdb",$name"."$counter".elec" >> color-elec.pml echo data2b_res $name"."$counter".elec",$pdbout$name"."$counter".elec" >> color-elec.pml echo color_b\(selection=\'\/$name"."$counter".elec"\',item=\'b\',gradient=\'rwb\',mode=\'hist\'\) >> color-elec.pml echo save $pdbout$name"."$counter".elec.pdb", \($name"."$counter".elec"\) >> color-elec.pml echo load $ligdir$name"."$counter".minbres.pdb",$name"."$counter >> color.pml echo data2q_res $name"."$counter,$pdbout$name"."$counter".vdw" >> color.pml echo data2b_res $name"."$counter,$pdbout$name"."$counter".elec" >> color.pml echo save $pdbout$name"."$counter".color.pdb", \($name"."$counter\) >> color.pml echo color_q\(selection=\'\/$name"."$counter\',item=\'b\',gradient=\'rwb\',mode=\'hist\',minimum=\'-1\',maximum=\'1\'\) >> color.pml echo color_b\(selection=\'\/$name"."$counter\',item=\'b\',gradient=\'rwb\',mode=\'hist\'\) >> color.pml #load_models /home/henricsn/combine2go/data/urokinase/anal/model/model5/*.elec.pdb,elec,discrete=1 #cmd.show("spheres" ,"resi 1:255") #sel ligand, (res 256) #cmd.show("sticks" ,"ligand") #util.cba(26,"ligand") ################################################################## #rm -f $name"."$counter"-temp.vdw" $name"."$counter"-temp.elec" cp -f $name"."$counter".seq" sequence.seq rm -f $name"."$counter".seq" done paste A??.?-temp.vdw > all.vdw paste A??.?-temp.elec > all.elec