VMD-L Mailing List

From: crockett c.h. (chc2g16) (chc2g16_at_soton.ac.uk)
Date: Mon Dec 09 2019 - 17:28:46 CST

Hi all,

Sorry to have to ask a very similar question to last time however I have now come across the same issue for calculating the RMSD per residue (averaged over the time course of the simulation)- for those if you who may have not seen the previous thread I started, I am new to running and analysing MD simulations. I have run 105ns simulations for a 153 amino acid protein solvated in a 10 angstrom water box but have had great difficulty analysing since the files produced since it is too big for my laptop to handle analysis via the VMD gui.

Thus far I have tried to use the script provided in the tutorial and edit it to allow for bigdcd (I have edited out my additions with ‘!’ to make it clear where I’ve added things). Further I also just need to analyse from 5ns to 105ns (excluding the first 5ns for equilibration – at a 2fs time steps and writing to the dcd every 250 steps so I guess from frame 10,000?) which I am not sure how to even start without having visual access to the trajectory in VMD. I would go about it in this way, without using .tcl scripts however VMD crashes on my computer every time I try and load the trajectory- I have a laptop and no access to a desktop compute with more RAM. Below is what I have tried so far but with no luck (and this is over the whole time frame rather than from 10,000)

# % $Id: residue_rmsd.tcl,v 1.4 2006/03/06 23:56:46 timisgro Exp $

!proc rmsd_residue_over_time { frame } {
! global reference compare all num_steps fil ref comp rmsd($r) res_b ave
!}
!source bigdcd.tcl
!set mol [mol new 2gbu_wbn.psf type psf waitfor all]

proc rmsd_residue_over_time {{mol top} res} {

    # use frame 0 for the reference
    set reference [atomselect $mol "protein" frame 0]
    # the frame being compared
    set compare [atomselect $mol "protein"]
    #make a selection with all atoms
    set all [atomselect top all]
    #get the number of frames
    set num_steps [molinfo $mol get numframes]
    #open file for writing
    set fil [open residue_rmsd.dat w]

    foreach r $res {
                set rmsd($r) 0
    }

    #loop over all frames in the trajectory
    for {set frame 0} {$frame < $num_steps} {incr frame} {
                puts "Calculating rmsd for frame $frame ..."
                # get the correct frame
                $compare frame $frame
        $all frame $frame
                # compute the transformation
                set trans_mat [measure fit $compare $reference]
                # do the alignment
                $all move $trans_mat

                # compute the RMSD
                #loop through all residues
                foreach r $res {
                    set ref [atomselect $mol "protein and resid $r and noh" frame 0]
                    set comp [atomselect $mol "protein and resid $r and noh" frame $frame]
                    set rmsd($r) [expr $rmsd($r) + [measure rmsd $comp $ref]]
                    $comp delete
                    $ref delete
                }
    }

    set ave 0
                foreach r $res {
                    set rmsd($r) [expr $rmsd($r)/$num_steps]
                    # print the RMSD
                    puts "RMSD of residue $r is $rmsd($r)"
                    puts $fil " $r \t $rmsd($r)"
                    set res_b [atomselect $mol "resid $r"]
            $res_b set user $rmsd($r)
            $res_b delete
                    set ave [expr $ave + $rmsd($r)]
                }

    set ave [expr $ave/[llength $res]]
    puts " Average rmsd per residue: $ave"
    close $fil
}

!rmsd_residue_over_time top $sel_resid
!bigdcd rmsd_residue_over_time 2gbu_18.dcd
!bigdcd_wait
!quit

Any help would be much appreciated as the results from these simulations are due to be presented imminently!

Many thanks again,
Catherine

Sent from Mail<https://go.microsoft.com/fwlink/?LinkId=550986> for Windows 10