VMD-L Mailing List

From: Ajasja Ljubetič (ajasja.ljubetic_at_gmail.com)
Date: Tue May 24 2011 - 09:48:15 CDT

Hi,

since bigdcd loads one frame from memory processes it and then deletes it
you only have one frame loaded.
Try changing the line
 set sel [atomselect $molid "name OW and within $cut of ($seltext)" frame
$i]
to
 set sel [atomselect $molid "name OW and within $cut of ($seltext)"]

and also delete
 $sel frame $i

In general you have to be change any selection where you are referencing the
frame number.

Best regards,
Ajasja

On Tue, May 24, 2011 at 15:39, maria goranovic <mariagoranovic_at_gmail.com>wrote:

> Hi
>
> I am trying to use Justin's script closewater.tcl to find the closest 1000
> water molecules to the protein in each frame of the trajectory. Although
> this works perfectly, my trajectory is too big for vmd to handle (16,000
> frames of 85,000 atoms, no velocities). So, I thought of combining
> bigdcd.tcl to closewater.tcl to get the data out. This needed some minor
> modifications in closewater.tcl. I have the following problem though.
>
> closewater.tcl reads in the entire trajectory, and calculates the number of
> waters using the following:
>
> set sel [atomselect $molid "name OH2 and within $cut of ($seltext)" frame
> $i]
>
> (where molid is 0)
>
> when I use bigdcd, I call closewater using only the frame number, and get
> the following error:
>
> bigdcd aborting at frame 3
> atomsel: frame 3 out of range for molecule 0
>
> How can this be sorted out?
>
> I am attaching the entire script that is being use below:
>
> ##############################################
> mol load gro temp.gro
> source mybigdcd.tcl
>
> # closewater.tcl
> # Justin Gullingsrud
> # jgulling_at_mccammon.ucsd.edu
> # 8 November 2004
>
> # This script processes a trajectory to create a new file containing just
> # a selection of atoms and the N closest water to that selection. The N
> # waters are recomputed for each timestep, and need not be the same waters
> # in each timestep (in fact, they probably will not be); thus it is in
> general
> # meaningless to analyze the dynamics of individual waters. However, it
> # may be useful for analyzing the distribution of waters around a
> relatively
> # static protein or DNA chain.
>
> # usage: closewater <molid> <selection text> <# waters> <filename prefix>
>
> proc closewater {i} {
>
> set molid 0
> set seltext "protein"
> set nwat 3000
> set prefix "close"
> set wat [atomselect $molid "name OW"]
> if { [$wat num] < $nwat } {
> error "Molecule contains ony [wat $num] water molecules."
> }
> # set nframes [molinfo $molid get numframes]
> # puts "the number of frames is $nframes"
> # for { set i 0 } { $i < $nframes } { incr i } {
> puts "now working in frame $i"
> set numinner 0
> set inner [list]
> set cut 1
> while {1} {
> set sel [atomselect $molid "name OW and within $cut of ($seltext)"
> frame $i]
> set outer [$sel list]
> $sel delete
> set numouter [llength $outer]
> if { $numouter < $nwat } {
> set inner $outer
> set numinner $numouter
> incr cut
> continue
> }
> break
> }
> puts "Found $numouter waters at cutoff $cut"
> # take all waters from inner, and just enough from outer
> catch { unset ohash }
> foreach ind $outer { set ohash($ind) 1 }
> foreach ind $inner { unset ohash($ind) }
> set outer [lrange [array names ohash] 0 [expr $nwat - $numinner - 1]]
> set watind [concat $inner $outer]
> set sel [atomselect $molid "($seltext) or same residue as (index
> $watind)"]
> $sel frame $i
> set j [format %06d $i]
> # $sel writepdb [format $prefix-%04d.pdb $i]
> $sel writepdb [format $prefix-$j.pdb]
> $sel delete
> # }
> }
>
> bigdcd closewater temp.xtc
> bigdcd_wait_till_done
>
> ##############################################
> --
> Maria G.
> Technical University of Denmark
> Copenhagen
>