VMD-L Mailing List
From: Maxim Belkin (mbelkin_at_ks.uiuc.edu)
Date: Sun Jan 25 2015 - 06:24:13 CST
- Next message: andrea spitaleri: "Re: optixRender error"
- Previous message: Antonija Tomić: "Hbond per residue"
- In reply to: Antonija Tomić: "Hbond per residue"
- Next in thread: Antonija Tomić: "Re: Hbond per residue"
- Reply: Antonija Tomić: "Re: Hbond per residue"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] [ attachment ]
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"
> }
> }
- Next message: andrea spitaleri: "Re: optixRender error"
- Previous message: Antonija Tomić: "Hbond per residue"
- In reply to: Antonija Tomić: "Hbond per residue"
- Next in thread: Antonija Tomić: "Re: Hbond per residue"
- Reply: Antonija Tomić: "Re: Hbond per residue"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] [ attachment ]