# Kipp <kippjohnson@uchicago.edu>
# for questions
#
# To run this file, from the TK console
# type "source RMSF.tcl" (without quotes)
# with a trajectory already loaded in VMD
#
# Script calculates RMSF of each 
# backbone CA in your molecule
#
# Set output name below (in quotes)
set outfile [open "rmsf.csv" w]

set sel [atomselect top "name CA"]

set nframes [molinfo top get numframes]
set nframes2 [expr $nframes - 1]

# Comment out below line if you do not want a header in output
puts $outfile "ResidueID, RMSF"

foreach ids [$sel get resid] { 
     set rmsf [measure rmsf [atomselect top "resid $ids"]] 
     puts $outfile "$ids , $rmsf" 
} 

close $outfile

# output will be file with two tab-seperated columns
# first column = residue number
# second column = RMSF
#
# the following R code will produce a quick plot
# Just input data file location
#
# par = mfrow(1,1)
# rmsf = read.table("<rmsf.dat-file-location>",header=TRUE, sep=",")
# plot(rmsf$V1, rmsf$V2, type="l", col="RED", xlab="Residue Number", ylab="RMSF (Å)")
# Kipp <kippjohnson@uchicago.edu>
# for questions
#
# To run this file, from the TK console
# type "source RMSF.tcl" (without quotes)
# with a trajectory already loaded in VMD
#
# Script calculates RMSF of each 
# backbone CA in your molecule
#
# Set output name below (in quotes)
set outfile [open "rmsf.csv" w]

set sel [atomselect top "name CA"]

set nframes [molinfo top get numframes]
set nframes2 [expr $nframes - 1]

# Comment out below line if you do not want a header in output
puts $outfile "Residue ID, RMSF"

foreach ids [$sel get resid] { 
     set rmsf [measure rmsf [atomselect top "resid $ids"]] 
     puts $outfile "$ids , $rmsf" 
} 

close $outfile

# output will be file with two tab-seperated columns
# first column = residue number
# second column = RMSF
#
# the following R code will produce a quick plot
# Just input data file location
#
# par = mfrow(1,1)
# rmsf = read.table("<rmsf.dat-file-location>",header=TRUE, sep=",")
# plot(rmsf$V1, rmsf$V2, type="l", col="RED", xlab="Residue Number", ylab="RMSF (Å)")
#