# 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
# vmd -dispdev text  -eofexit titin-per-res-rmsd.tcl

package require timeline 2.0

set outFilename titin-per-res-rmsd.tml

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

set myMolid 0

proc myRmsdBatchCalc {filename molid} {

  # 0 for per-residue calcs
  # 1 for per-selection calcs
  set usesFreeSelection 0
  # include units in title
  set dataTitle "res. RMSD (A)"
  set firstFrame 0
  set lastFrame 96
  # the sample molecule is a 1-chain, 1-segment molecule, so we will do
  # a simple loop over residue numbers
  set theChain "T"
  set theSeg "TIT"
  set firstRes 1
  set lastRes 89

  set numFrames [expr $lastFrame - $firstFrame + 1]
  set numSelectionGroups [expr $lastRes-$firstRes + 1]




  if {$filename == ""  } {
      die "usage: myRMSDBatchCalc FILENAME MOLID\n FILENAME cannot be an empty string."
  }

  set outDataFile [open $filename w]



  ::timeline::writeDataFileHeader $outDataFile $molid $dataTitle $numFrames \
     $numSelectionGroups $usesFreeSelection

  set chain $theChain
  set seg $theSeg
  for {set r $firstRes} {$r <= $lastRes} {incr r} {
    set sela [atomselect $molid "resid $r"]
    set selb [atomselect $molid "resid $r"]
    $sela frame 0
    for {set f $firstFrame} {$f<=$lastFrame} {incr f} {
      $selb frame $f
      display update
      set val [measure rmsd $sela $selb]
      set resid $r
      puts $outDataFile "$resid $chain $seg $f $val"
    }
  }
  close $outDataFile
  return
}

myRmsdBatchCalc $outFilename $myMolid
