VMD-L Mailing List

From: Antonija Tomić (piculast_at_gmail.com)
Date: Sun Jan 25 2015 - 03:56:19 CST

 Dear VMD Users,

I've written a script to calculate a number of hydrogen bonds per residue
per frame (see below text). As shown my protein consists of 723 amino acids
residues and the trajectory is composed of 2500 frames. When I execute
script for all 723 residues (from 1st till 723th residue, as show in
script below) the output files (res$startres-hbond.dat) contain information
only for 2397 frames per residue, but when I consider only ~50 residues
(for example from 350th till 400th) I get output for all 2500 frames per
residue. I must emphasize that I have checked if I have enough free disk
space and enough memory (16 Gb).

I would appreciate if you could tell me where the problem is,

Antonija

-------------------------------

set mol [molinfo top]

set nf [molinfo $mol get numframes]

for {set startres 1} {$startres <= 723} {incr startres} {

set aa [atomselect $mol "(mass 14.01 16 32.06) and (resid $startres)"]

set prot [atomselect $mol "(mass 14.01 16 32.06) and (all not resid
$startres)"]

set outfile [open res$startres-hbond.dat w]

    for {set i 0} {$i < $nf} {incr i} {

        $aa frame $i

        $prot frame $i

        $aa update

        $prot update

        set nhb1 [llength [lindex [measure hbonds 3.5 60 $aa $prot] 0]]

        set nhb2 [llength [lindex [measure hbonds 3.5 60 $prot $aa] 0]]

        puts $outfile "$i $nhb1 $nhb2"

    }

}