[gmx-users] Trouble loading User data with a Tcl script (only half the data seems to load)

[J. Nathan Scott]

Hello all,

I have written a script that does per-time-step coloring of the water
molecules, residues, and ions in my protein based on electric field values I
previously calculated. I received excellent help here a few months sorting
out some problems with the script, but there is one big one remaining that I
hope someone can help with.

The problem is that only half the frames (the first half) in my simulation
are getting colored, and I can't figure out why. This happens no matter how
many frames I load, whether the full trajectory or just a few hundred
frames. For what it's worth, the machine I'm running on has 12 GB of RAM, a
very nice video card, and a new Xeon quad core processor, so I don't think
this is an issue of system resources. At least I can see that there is still
tons of free RAM when I am experiencing this issue.

Also, I would sincerely appreciate any advice on speeding up one of my
loops. To color the water molecules and ions I can simply atomselect them
all, and then loop through the values from the electric field data file for
that timestep and assign the color value to a user list either 3 or 1 times
respectively to color each atom.

However, for the residues in the protein I am having to atomselect them one
at a time, get the number of atoms for that residue, and then do another
loop to build a list of the correct length to assign the user values to each
atom of the residue. Can anyone recommend a more efficient way of coloring
the residues? When the residue coloring part of the script is stripped out
it runs so much faster despite the fact that there are far fewer protein
atoms than there are water or ion atoms, so I know this method of assignment
is sluggish.

Please see the script below my signature, and thank you in advance for any
advice you can provide.

-- 
----------
J. Nathan Scott, Ph.D.
Postdoctoral Fellow
Department of Chemistry and Biochemistry
Montana State University




#######################################################################
set first 0;
set last 1000;
set mut wt;
set mut_ wt_;
set i $first; #i will be timestep/filename indicator,
set j 0;

mol new /data/1stn/xtc/1stn_$mut.gro type gro waitfor all
mol addfile /data/1stn/xtc/1stn_$mut.xtc type xtc waitfor all first $first
last $l$
set mol_ID top;
set n [ molinfo $mol_ID get numframes ];
animate goto 0
animate delete  beg 0 end 0 skip 0 0
mol delrep 0 $mol_ID

while {$i <= $last} {
        set fp [open "/data/1stn/$mut/pd5/pd5.stripped/1stn_$mut_$i.pd5" r]
        set file_data [read $fp]
        close $fp

        set data [split $file_data "\n"]
        foreach {one} $data {
                lappend efield [lindex $one 3]
        }

### Looping through residues in the for loop,
### selecting one residue at a time, getting its number of atoms,
### and then building a user value list for each of those
### atoms in another for loop.

        for {set k 6} {$k <= 141} {incr k} {
        set prot [atomselect $mol_ID "resid $k" frame $j]
        set num [$prot num]
        $prot frame $j
        set user_list {}
        set u [lindex $efield [expr $k - 6]]
        for {set b 1} {$b <= $num} {incr b} {
           lappend user_list $u
        }
        $prot set user $user_list
        $prot delete
        }

### Water loop works well

        set wat [atomselect $mol_ID waters frame $j]
        $wat frame $j
        set user_list {}
        for {set a 136} {$a < 10233} {incr a} {
                set u [lindex $efield $a]
                lappend user_list $u $u $u
        }
        $wat set user $user_list
        $wat delete

### Ion loop works well too

        set ions [atomselect $mol_ID ions frame $j]
        $ions frame $j
        set user_list {}
        for {set a 10233} {$a < 10243} {incr a} {
                set u [lindex $efield $a]
                lappend user_list $u
        }
        $ions set user $user_list
        $ions delete

        unset efield
        incr j 1
        incr i 2
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110504/2e712b12/attachment.html>