tk/tcl的script

我写的用来计算一个二聚体蛋白之间相靠近的tk/tcl的script

有什么批评和建议请尽管提出来！：）

proc Inter_Dist_Cal {} {

set di "1Q2W.pdb"
puts -nonewline "Enter PDB file (press "d" as "$di"): "
flush stdout
set infile [gets stdin]
if {[string equal -nocase $infile "d"]} {set infile $di}

set PDBname [lindex [split $infile .] 0]

set do1 "($PDBname)All_Inter_Con.txt"
puts "All Inter Contact will be "$do1" "
flush stdout
set outfile1 $do1

set do2 "($PDBname)Res_Inter_Con.txt"
puts "Residual Inter Contact will be "$do2" "
flush stdout
set outfile2 $do2

set do3 "($PDBname)Num_Res_Inter_Con.txt"
puts "The number of Inter Contact for each Residue will be "$do3" "
flush stdout
set outfile3 $do3

set dd "10"
puts -nonewline "Enter max distance (press "d" as "$dd"): "
flush stdout
set indist [gets stdin]
if {[string equal -nocase $indist "d"]} {set indist $dd}

set in [open $infile r]
set n1 0
set n2 0
while {[gets $in line]>=0} {
#set line [lindex [split $line """] 0]
if {([lindex $line 0]==" ")||($line=="")||([lindex [split [lindex $line 0] ""] 0]=="#")}
{continue}
if {[lindex $line 0]!="ATOM"} {continue}
if {[lindex $line 4]=="A"} {
   set 1id($n1) [lindex $line 1]
   set 1atom($n1) [lindex $line 2]
   set 1resty($n1) [lindex $line 3]
   set 1resno($n1) [lindex $line 5]  
   set 1x($n1) [lindex $line 6]  
   set 1y($n1) [lindex $line 7]  
   set 1z($n1) [lindex $line 8]  
   incr n1
   }
if {[lindex $line 4]=="B"} {
   set 2id($n2) [lindex $line 1]
   set 2atom($n2) [lindex $line 2]
   set 2resty($n2) [lindex $line 3]
   set 2resno($n2) [lindex $line 5]  
   set 2x($n2) [lindex $line 6]  
   set 2y($n2) [lindex $line 7]  
   set 2z($n2) [lindex $line 8]  
   incr n2
   }  
   }
close $in

set out1 [open $outfile1 w]
puts $out1 "# [clock format [clock seconds]]  BY [lindex [exec who] end-4]"
puts $out1 ""
puts $out1 [format "# %4s %5s %5s - %6s %5s %5s %9s" Mol1 Res Atom Mol2 Res Atom Distance]
puts $out1 ""

for {set i 0} {$i<$n1} {incr i} {
    for {set j 0} {$j<$n2} {incr j} {
         if {[set dis [expr pow(pow($1x($i)-$2x($j),2)+pow($1y($i)-$2y($j),
2)+pow($1z($i)-$2z($j),2),0.5)]]<=$indist} {
        set fm [format "%6d %5s %5s - %6d %5s %5s %9.3f" $1resno($i) $1resty($i) $1atom($i)
$2resno($j) $2resty($j) $2atom($j) $dis]
        puts $out1 $fm
        set temp1 $1resno($i)
        if {[info exists ACon($temp1)]} {
    incr ACon($temp1)
    }  else {set ACon($temp1) 1}
        set temp2 [expr $1resno($i)*1000+$2resno($j)]
    if {[info exists  ABCon($temp2)]} {
       incr  ABCon($temp2)
               } else {set ABCon($temp2) 1}
     }
        }
    }
close $out1

set out2 [open $outfile2 w]
puts $out2 "# [clock format [clock seconds]]  BY [lindex [exec who] end-4]"
puts $out2 ""
puts $out2 [format "# %4s %5s  %6s %5s %9s" Mol1 Res Mol2 Res Contact]
puts $out2 ""

for {set i 0} {$i<$n1} {incr i} {
    for {set j 0} {$j<$n2} {incr j} {
         if {[info exists  ABCon([expr $1resno($i)*1000+$2resno($j)])]} {
      set fm [format "%6d %5s %6d %5s %9d" $1resno($i) $1resty($i) $2resno($j) $2resty($j)
$ABCon([expr $1resno($i)*1000+$2resno($j)])]
             puts $out2 $fm
      unset ABCon([expr $1resno($i)*1000+$2resno($j)])
 
     }
        }
     }
close $out2 
set out3 [open $outfile3 w]
puts $out3 "# [clock format [clock seconds]]  BY [lindex [exec who] end-4]"
puts $out3 ""
puts $out3 [format "# %4s %5s %9s" Mol1 Res Contact]
puts $out3 ""

for {set i 0} {$i<$n1} {incr i} {
       if {[info exists ACon($1resno($i))]} {
          set fm [format "%6d %5s %9d" $1resno($i) $1resty($i) $ACon($1resno($i))]
          puts $out3 $fm
   unset ACon($1resno($i))
   }
     }
close $out3   

}