# sample vmd data processing script for Timeline
# Barry Isralewitz
#
# create a .tml data file for reading by Timeline
# use from commandline for pre-processing
# then read data in later
 
#temporary -- get header file definition
package require timeline 2.0

set theFilename titin-strand-contacts.tml

set theMolid 0

# load the file
mol new titin.psf type psf first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all
mol addfile titin.dcd type dcd first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all





proc myBatchDataWrite {filename molid} {
  # Instructions: edit this procedure for custom data analysis routines
  # The HEADER information is read by timeline 
  #
  # 0 for residue based calcs
  # 1 for selection based calcs
  set usesFreeSelection 1
  
  # variable dataName
  # variable trajdat
  # variable trajdatNum
  # variable numFrames

  set dataTitle "# contacts"
  set firstFrame 0
  set lastFrame 96 
  set numFrames [expr $lastFrame - $firstFrame + 1]
  set labeledSels  [list {"strand A" "resid 4 to 7"} \
          {"strand A'" "resid 11 to 15"} \
          {"strand B" "resid 18 to 26"}  \
          {"strand C" "resid 32 to 36"}  \
          {"strand D" "resid 47 to 52"}  \
          {"strand E" "resid 55 to 61"}  \
          {"strand F" "resid 69 to 75"}  \
          {"strand G" "resid 78 to 88"}]

  set numSelectionGroups [llength $labeledSels]

  if {$filename == ""  } {
    die "usage: batchDataWrite FILENAME MOLID"
   } 

  set outDataFile [open $filename w]

  
  # Calculation takes place here 
  # Also, assigns the freeSelLabel and freeSelString for each selection group. 
  set protSel [atomselect $molid "protein"]
  set i 0
  foreach {selPair} $labeledSels {
    foreach {label seltext} $selPair {}
    puts "starting sel. label: >>$label<<  seltext: >$seltext<"
    set trajdat(freeSelLabel,$i) $label 
    set trajdat(freeSelString,$i) $seltext 
    set sel [atomselect $molid $seltext]
    for {set trajFrame $firstFrame} {$trajFrame <= $lastFrame} {incr  trajFrame} {
      $sel frame $trajFrame
      $protSel frame $trajFrame
        
      # Calculation...
      set d  [llength [lindex [measure contacts 4.0 $sel $protSel] 0]] 

      set trajdat($trajFrame,$i) $d
      puts "$seltext, $trajFrame $i $d"
    }
    $sel delete
    incr i
  }

  set defaultVal 0 

  ::timeline::writeDataFileHeader $outDataFile $molid $dataTitle $numFrames $numSelectionGroups $usesFreeSelection 
 
  for {set i 0} {$i<$numSelectionGroups} {incr i} {
    puts $outDataFile "freeSelLabel $trajdat(freeSelLabel,$i)"
    puts $outDataFile "freeSelString $trajdat(freeSelString,$i)"
    for {set trajFrame 0} {$trajFrame < $numFrames} {incr  trajFrame} {
      if {![info exists trajdat($trajFrame,$i)]} {
         set v $defaultVal
      } else {
         set v $trajdat($trajFrame,$i)
      }
      puts $outDataFile "$trajFrame $v" 
    }
  } 

  close $outDataFile
  return
}



myBatchDataWrite $theFilename $theMolid
