# version 1.3
# extract distribution of residues represented by aname and rname around sel1
#
# usage: getResDistrib "sel1" "aname" "rname" min_distance max_distance step output_prefix <startframe <endframe>>
#
# script counts number of residues of given type within range (min-max) 
# of distances from sel1 structure
#
# rname - residue name, script counts number of aname atoms
# min / max - minimal and maximal analyzed distance
# step - distance step
# optional parameters: start and endframe in the trajectory, if not used, whole
# trajectory is analyzed
# frames are numbered from 1
#
# example:
# getResDistrib "protein" "S1" "SO4" 1.0 50.0 0.5
# analyze distribution of SO4 residues (count S1 atoms) around protein
# start at 1.0A up to 50.0A with 0.5A steps, 
# example2:
# fResPair "protein" "S1" "SO4" 1.0 50.0 0.5 1 10
# analogous, only perform analysis from frames 1 to 10 only
# output goes to output_prefix-distrib.dat file and has following format:
# distance	cumuative sum of residues within this distance of sel1

proc getResDistrib { sel1 aname rname min max step output_prefix args } {
  puts "Analyzing distribution of $rname ($aname representation)"
  puts "Distance limits $min A - $max A, step $step A"

  # get number of frames
  set numframes [molinfo top get numframes]
  puts "total number of frames in trajectory: $numframes"
  if {[llength $args] == 0} {
    set startframe 1
    set endframe $numframes
  }
  if {[llength $args] == 1} {
    set startframe [lindex $args 0]
    set endframe $numframes
    set stepframe 1
    if {$startframe <= 0 || $startframe > $numframes} {
      puts "illegal value of startframe, changing to first frame"
      set startframe 1
    }
  }
  if {[llength $args] >= 2} {
    set startframe [lindex $args 0]
    set endframe [lindex $args 1]
    if {$startframe <= 0 || $startframe > $numframes} {
      puts "illegal value of startframe, changing to first frame"
      set startframe 1
    }
    if {$endframe < $startframe || $endframe > $numframes} {
      puts "illegal value of endframe, changing to last frame"
      set endframe $numframes
    }
  }

  set totframes [expr ($endframe - $startframe + 1)]
  puts "analysis will be performed on $totframes frame(s) ($startframe to $endframe)"

  set file "$output_prefix-distrib.dat"
  set fw [open $file w]

  # now subtract 1 from all frame indexes - numbering starts with 0
  set startframe [expr $startframe - 1]
  set endframe [expr $endframe - 1]

  for { set dc $min } { $dc <= $max } { set dc [expr { $dc + $step }] } {
    set rescount 0.0
    puts "analyzing distance $dc A"
    set sel [atomselect top "(same residue as (resname $rname and within $dc of ($sel1) )) and name $aname"]
    set n [molinfo top get numframes]
    for {set i $startframe; set d 1} {$i <= $endframe} {incr i; incr d} {
      # show activity
      if { [expr $d % 50] == 0 }  {
        puts -nonewline "."
        if { [expr $d % 2000] == 0 } { puts " " }
        flush stdout
      }
      # analyze distribution
      $sel frame $i
      $sel update
      set rescount [ expr { $rescount + [$sel num]} ]
    }
    set rescount [ expr { $rescount / $totframes} ]
    puts $fw "$dc $rescount"
    flush $fw
    $sel delete    
    puts " "
  }
  close $fw
  puts "done"
}