Calculating RMSF in VMD

As a part of my coursework I have to calculate root mean square of fluctuation (RMSF) of a protein trajectory over the last 5 nanoseconds of simulation. I have used the inbuilt measure rmsf function of VMD. However, I want to write a script by myself to get a better understanding of both VMD and Tcl. This is the script I have written:

set outfile [open ~/rmsf.dat w]
# trajectory is at top
set frame_num [molinfo top get numframes]
# calculating minimum frame of 5ns window
set min_frame [expr {$frame_num - ($frame_num / 5) - 1}]
set full [atomselect 0 "protein"] # selecting reference pdb
set prot [atomselect top "protein"]
# create empty list to store rmsf values
set rmsf {}
foreach i [$full get index] { lappend rmsf 0 } #iterate over frames for {set j$min_frame} {$j <$frame_num} {incr j} {
$prot frame$j
$prot update # iterate over each atom foreach i [$full get index] {
lset rmsf $i [expr {[lindex$rmsf $i]} + {[vecdist [lindex [$full get {x y z}] $i] [lindex [$prot get {x y z}] $i]]}] } } # putting output foreach i [$full get index] {
puts $outfile [expr {sqrt([lindex$rmsf $i]/250)}] }$prot delete
$full delete close$outfile


My calculated rmsf values differ from the default calculated values. I guess this can be because I have used a different reference compared to default, or I calculated it wrong.

Which frame does VMD use as a reference in measure rmsf the first frame in selection or the first frame of the entire trajectory?

• Thanks for bringing this over from Bioinformatics.SE, and welcome to the site!! +1 and we hope to see much more of you here!! I'll just point out that we prefer only 1 question per post. So maybe pick the one you most want answered, and then ask the other two in separate questions. You can provide a link to this question in your other two questions, so that the script you showed can be reused. – Nike Dattani Dec 7 '20 at 3:27

The problem in the sript I had written was that my system was not able to handle that many variables at one go. So I made the calculations piecewise i.e. 40 atoms at a time and the results closely match the values one gets from the inbuilt function. In any case, calculating RMSF for just the C alpha atoms gives more information about the structural changes in the protein than measuring it for all atoms. I wrote the following script to do the same:

set outfile [open savefile.dat w]
set ref [atomselect 0 "protein"]
set sel [atomselect top "protein"]
for { set f 0 } { $f < [molinfo top get numframes] } { incr f } {$sel frame $f$sel move [measure fit $sel$ref]
}
$sel delete$ref delete
set sel [atomselect top "name CA"]
set val [measure rmsf $sel] for {set i 0} {$i < [llength $val]} { incr i} { puts$outfile [lindex $val$i]
}
$sel delete close$outfile


I also figured out the inbuilt measure rmsf function uses the average of the whole trajectory (that is why I aligned the trajectories in the script) as a reference rather than a specific frame. One can actually modify this by going into /src/Measure.C file in the source code if required.