VMD-L Mailing List

From: Antonija Tomić (piculast_at_gmail.com)
Date: Mon Jan 26 2015 - 03:08:07 CST

Dear Maxim,

 thank you very much for your suggestions. I modified my script according
to yours suggestions:

 set nf [molinfo top get numframes]

set all [atomselect top all]

set resids [lsort -unique [$all get resid]]

foreach startres $resids {

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

set prot [atomselect top "(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"

}

$aa delete

$prot delete

close $outfile

}

 But I am wondering why I'm still having different results for some
residues depending whether I execute the above script or just “for” part of
the above script (manually defining one residue, script shown below).

 Kind regards,

Antonija

 set outfile [open res100-hbond.dat w]

set mol [molinfo top]

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

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

set nf [molinfo $mol get numframes]

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"

}

2015-01-25 13:24 GMT+01:00 Maxim Belkin <mbelkin_at_ks.uiuc.edu>:

> My guess would be that you don’t close the file you open at the end of the
> loop over residues AND you don’t flush the output, so you just don’t see
> the last bit of frames.
>
> Here is a list of problem with your script:
>
> 1. You don’t have "$aa delete" and “$prot delete” at the end of the loop
> over residue numbers.
> 2. You don’t close the opened file (close $outfile) at the end of the loop
> over residue numbers.
> 2. Use more robust loop over residues:
> set all [atomselect top all]
> set resids [lsort -unique [$all get resid]]
> foreach startres $resids { ….
>
> Maxim
>
> On Jan 25, 2015, at 03:56, Antonija Tomić <piculast_at_gmail.com> wrote:
>
> 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"
>
> }
>
> }
>
>
>