VMD-L Mailing List

From: maria goranovic (mariagoranovic_at_gmail.com)
Date: Tue May 24 2011 - 08:39:08 CDT

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