From: J. Nathan Scott (scottjn_at_chemistry.montana.edu)
Date: Wed Feb 23 2011 - 12:10:38 CST

Hi Jérôme,

Thank you so much for your help, your example solved my issue
perfectly and now the script is working beautifully and updating the
water colors at every step.

Here is the working script (which I have not yet had time to comment
or turn into a procedure, but could still be useful to someone else
trying to accomplish per-atom per-timestep water coloring) in case it
might be of some use to another VMD user. The visualizations in this
script are rather crude, but could be modified to be much nicer. Thank
you again so much for sharing your expertise Jérôme!

Best Wishes,
Nathan

###################################################################
set i 0; #i will be timestep/filename indicator, j=i/2 will also be
the current frame

mol new {E:\\scottjn\\Data\\1stn\\wt\\1stn_wt_0.g96} type g96 waitfor all
mol addfile {E:\\scottjn\\Data\\1stn\\wt\\1stn_wt.xtc} type xtc
waitfor all first 0 last 500
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
mol representation CPK 0.350000 1.000000 8.000000 6.000000
mol selection {waters within 5 of resid 140}
mol addrep $mol_ID
mol colupdate 0 $mol_ID 1
mol selupdate 0 $mol_ID 1
mol color User

mol representation CPK 0.350000 1.000000
mol color Name
mol selection {resid 140}
mol addrep $mol_ID
mol colupdate 1 $mol_ID 1
mol selupdate 1 $mol_ID 1
mol modcolor 1 0 ColorID 7

while {$i < 1002} {
        set j [expr $i / 2]; #j is the current VMD frame number
        if {$i < 10} {
                    set fp [open "E:\\scottjn\\Data\\1stn\\wt\\wat\\water_step_0000$i.wat" r]
                    incr i 2
        } elseif {$i < 100 && $i >= 10 } {
                    set fp [open "E:\\scottjn\\Data\\1stn\\wt\\wat\\water_step_000$i.wat" r]
                    incr i 2
        } elseif {$i < 1000 && $i >= 100 } {
                    set fp [open "E:\\scottjn\\Data\\1stn\\wt\\wat\\water_step_00$i.wat" r]
                    incr i 2
        } elseif {$i < 10000 && $i >= 1000 } {
                    set fp [open "E:\\scottjn\\Data\\1stn\\wt\\wat\\water_step_0$i.wat" r]
                    incr i 2
        } elseif {$i >= 10000 } {
                    set fp [open "E:\\scottjn\\Data\\1stn\\wt\\wat\\water_step_$i.wat" r]
                    incr i 2
        }
        set file_data [read $fp]
        close $fp
        set data [split $file_data "\n"]
        foreach {one} $data {
                lappend shift [lindex $one 1]
        }

        set wat [atomselect $mol_ID waters frame $j]
        $wat frame $j
        set user_list {}
        for {set a 0} {$a < [expr [llength $shift] - 1]} {incr a} {
                set u [lindex $shift $a]
                lappend user_list $u $u $u
        }
        $wat set user $user_list
        $wat delete
        unset shift
}

color scale midpoint 0.379590
color scale min 0.000000
color scale max 1.000000
mol modcolor 0 0 User
mol scaleminmax 0 0 -12.932960 21.137920

On Thu, Feb 17, 2011 at 3:10 AM, Jérôme Hénin <jhenin_at_ifr88.cnrs-mrs.fr> wrote:
> Hi Nathan,
>
> In your script, the key line is: set wat [atomselect $mol_ID "resid $k"]
> It creates a new atomselect object, which takes a little while and
> allocates memory. If you do it in a loop, you get poor performance,
> and (unless you deallocate the objects with "$wat delete") a memory
> leak, which is what eventually crashes VMD.
>
> The best way is to create a single selection with all atoms, and to
> set the user field from Tcl list, containing one value per atom.
>
> set wat [atomselect $mol_ID <all the requested atoms> frame 0]
>
> for { <loop i on frames> } {
>  $wat frame $i
>  set user_list {}
>  for { <loop on data> } {
>    set u <get data for one molecule>
>    lappend user_list $u $u $u  # one value per atom
>  }
>  $wat set user $user_list
> }
>
> $wat delete
>
>
> Cheers,
> Jerome
>
> On 17 February 2011 00:17, J. Nathan Scott
> <scottjn_at_chemistry.montana.edu> wrote:
>> Hello there,
>>
>> I have the following procedure I would like to implement in VMD, if
>> possible, and would very much appreciate any feedback on its
>> feasibility or methods. I have looked through the previous posts I
>> could find that focused on coloring in VMD, but none of them seemed
>> particularly on point for this issue.
>>
>> For a 40 ns Gromacs (20001 frames) simulation I would like to color
>> the water molecules of each frame by the electric field they generate
>> at a certain point in a protein of interest. This information has
>> already been calculated for each frame, and what I have right now (in
>> addition to the trajectory) is a series of ASCII files named
>> water_step_xxxxx.wat each with a single column containing the magnitude and
>> sign of the electric field. The water information in the ASCII files
>> is ordered the same as the real ordering of the water molecules in the
>> Gromacs topology file and trajectory.
>>
>> Is this something that is possible? For what it's worth, there are
>> approximately 10000 water molecules in the simulation cell. From the
>> reading I have been doing it seems like this might not run dynamically
>> in VMD, which would be OK as long as I could generate images showing
>> the waters colored by electric field magnitude and sign.
>>
>> My Tcl script thus far, which has been cobbled together from other
>> scripts and the bits of Tcl I've manged to pick up
>> looks like this:
>> [code]
>> # First water is resid 142 for the 1stn simulations
>> # Largest - shift = -12.93296 #WT
>> # Largest + shift =  21.13792 #WT
>>
>> mol new {C:\\wat\\1stn_wt_0.g96} type g96 waitfor all
>> mol addfile {C:\\wat\\1stn_wt_0-5ns.xtc} type xtc waitfor all
>> animate goto 0
>> animate delete  beg 0 end 0 skip 0 0
>> mol delrep 0 top
>> mol representation Lines 1.000000
>> mol color User
>> mol selection {water}
>> mol addrep top
>> mol colupdate 0 top 1
>>
>> mol representation CPK 0.350000 1.000000
>> mol color Name
>> mol selection {protein}
>> mol addrep top
>> mol colupdate 1 top 1
>>
>> set mol_ID top;
>> set n [ molinfo $mol_ID get numframes ];
>>
>> set numwaters 10097; #only correct for WT 1stn
>> set i 0; #i will be timestep/filename indicator, i/2 will also be the
>> current frame
>>
>> # This loop should do all of the file reading and data scaling and
>> color changing.
>> while {$i < 10} {
>>        set j [expr $i / 2]; #j is the current VMD frame number
>>        if {$i < 10} {
>>                set fp [open "C:\\wat\\water_step_0000$i.wat" r]
>>                incr i 2
>>        } elseif {$i < 100 && $i >= 10 } {
>>                set fp [open "C:\\wat\\water_step_000$i.wat" r]
>>                incr i 2
>>        } elseif {$i < 1000 && $i >= 100 } {
>>                set fp [open "C:\\wat\\water_step_00$i.wat" r]
>>                incr i 2
>>        } elseif {$i < 10000 && $i >= 1000 } {
>>                set fp [open "C:\\wat\\water_step_0$i.wat" r]
>>                incr i 2
>>        } elseif {$i >= 10000 } {
>>                set fp [open "C:\\wat\\water_step_$i.wat" r]
>>                incr i 2
>>        }
>>        set file_data [read $fp]
>>        close $fp
>>        set data [split $file_data "\n"]
>>
>> #All Trp shifts for file $i are in the list "data"
>>        set a 0;
>>        set k 142; #first water resid
>>
>>        while {$a <= [expr $numwaters - 1 ]} {
>>                set wat [atomselect $mol_ID "resid $k"]
>>                $wat frame [expr $i / 2]
>>
>>                if {[expr [lindex $data $a]] <= 0} {
>>                set [lindex $data $a] [expr [lindex $data $a]/-12.93296]
>>                $wat set user [lindex $data $a]
>>                }
>>                elseif {[expr [lindex $data $a]] > 0} {
>>                set [lindex $data $a] [expr [lindex $data $a]/21.13792]
>>                $wat set user [lindex $data $a]
>>                }
>>
>>                incr a
>>                incr k
>>        }
>> }
>> [/code]
>>
>> Is there perhaps a better way to do something like this? As I have it
>> so far VMD crashes after a minute or so and also, though the "user"
>> field of the water molecules is presumably changing, I'm not seeing
>> any change in the color of any of the water molecule's when I just
>> work with a single frame/.wat file and copy and paste the relevant
>> commands from the script I pasted above. Is there something else I
>> need to do to get the colors to update? Finally, is coloring by "user"
>> similar to the other coloring methods where the value should range
>> from 0 to 1? That is what I had assumed, and thus my dividing every
>> water's electric field value by the maximum positive or negative
>> across the whole trajectory, but I'm not certain this is correct.
>>
>> Thanks sincerely in advance for any tips or insight you can provide.
>> I've been beating my head against this for about a week now and could
>> definitely use a little expert advice!
>>
>> Best wishes,
>> Nathan
>>
>> ----------
>> J. Nathan Scott, Ph.D.
>> Postdoctoral Fellow
>> Department of Chemistry and Biochemistry
>> Montana State University
>>
>

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