VMD-L Mailing List

From: Peter Freddolino (petefred_at_ks.uiuc.edu)
Date: Mon Apr 23 2007 - 10:57:10 CDT

Hi Jianhui,

if you still run into problems like this, it is usually useful to add
some print/puts statements throughout the code to see what is actually
going on. My suspicion is that things aren't getting "stuck", but that
the program is indeed running a good chunk of your script.

I'm not necessarily saying you need to get all of your data into lists
and then do more analysis -- that's one possible approach, but you could
just as easily leave the loops the way you had them arranged. But if you
are memory constrained, realize that aside from the dcd itself, the
biggest memory hog in analysis like this is extra atomselect objects
floating around, so it makes much more sense to create and destroy them
as needed than to create all of them up front and keep them around.

Unfortunately I'm not aware of a tcl builtin for set intersection; you'd
probably need to write it yourself.

Peter

Jianhui Tian wrote:
> Hi Peter,
>
> When I run the script with vmd -dispdev text -e script, it get killed
> after the output "atomselect0", which should come from "set systm
> [atomselect top all]". Then the program will stuck for a while and
> "get killed". When you say "You can use a list or array to store the
> data and unpack it at the end as desired." do you mean I can get all
> the lists for all the frames and all the residues and then do the
> analysis? Because what I need to do is to process all the lists per
> frame. And also I can release all the dcd trajectory in memory with
> only lists remained, which should be more memory efficient. Am I right?
>
> Is there a simply command line to check whether two list have any
> element in common and to combine two lists have elements in common to
> one "unique" list? Thank you very much.
>
> Jianhui
>
>
> On 4/23/07, *Peter Freddolino* <petefred_at_ks.uiuc.edu
> <mailto:petefred_at_ks.uiuc.edu>> wrote:
>
> Hi Jianhui,
> while we'd need a few more details (such as where it dies) to
> figure out
> exactly what is going on, I can recommend a couple of changes that
> should make this easier and more memory efficient:
>
> -Use lsort -unique instead of bothering to define a new procedure for
> removing duplicate elements from a list
>
> -Write your main analysis as nested loops instead of first
> creating all
> the selections you could possibly need and iterating over them, eg:
>
> for {set i 1} {$i <= 1536} {incr i} {
> set thiswatsel [atomselect top "water and within 5.5 of resid $i"]
>
> for {set j 0} {$j < $numfrm} {incr j} {
> $thiswatsel frame $j
> $thiswatsel update
> set closeresids [lsort -unique -integer [$thiswatsel get resid]]
> .....
>
> }
>
> $thiswatsel delete
> }
>
> You can use a list or array to store the data and unpack it at the
> end
> as desired.
>
> Peter
>
>
> Jianhui Tian wrote:
> > Hi all,
> >
> > I have a water box of 1536 molecules. Now what I want to do is
> to make
> > atomselect of each of the molecule and then get the residue list of
> > other molecules within some distance of each water molecule. The
> > problem is that the tcl code will get killed. I dont know the
> reason.
> > Is it due to the memory limitation. I only have 512 memory in my PC
> > now. How can I do what I want to do here? I also posted the code I
> > write about this part.
> >
> >
> > #########################################
> > #########################################
> > proc simplifylist {slist} {
> > set newlist { }
> > foreach elem $slist {
> > if {[lsearch $newlist $elem] == -1} {
> > lappend newlist $elem
> > }
> > }
> > set newlist [lsort $newlist]
> > return $newlist
> > }
> >
> > ########################################
> > # Load the trajectoty files.
> > mol new prmtop type parm7
> > mol addfile dcd type dcd waitfor all
> >
> > set numfrm [molinfo top get numframes]
> > set systm [atomselect top all]
> >
> > #######################################
> > #Get all the atomselect.
> > for {set i 1} {$i <= 1536} {incr i} {
> > set atom($i) [atomselect top "water within 5.5 of resid $i"]
> > }
> >
> > for {set i 0} {$i < $numfrm} {incr i} {
> >
> > ######### Get all the connected list for each of the water
> molecule.
> > for {set j 1} {$j <= 1536} {incr j} {
> > $atom($j) frame $i
> > $atom($j) update
> > set list($j) [$atom($j) get resid]
> > set list($j) [simplifylist $list($j)]
> > puts "frame $i residue $j: $list($j)\n"
> > }
> > }
> >
> > exit
> >
> >
>
>