VMD-L Mailing List

From: Maxim Belkin (mbelkin_at_ks.uiuc.edu)
Date: Mon Jan 26 2015 - 07:30:47 CST

Scripts seems to be similar, except that I don’t see “close $outfile” after the loop in the excerpt which could result in an un-written buffer at the end of the file.

> On Jan 26, 2015, at 03:08, Antonija Tomić <piculast_at_gmail.com> wrote:
>
> 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 <mailto: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 <mailto: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"
>> }
>> }
>
>