Calculating RMSF in VMD

Question

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?

Answer 1

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.