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?

  • $\begingroup$ 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. $\endgroup$ – 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.