VMD-L Mailing List

From: Jérôme Hénin (jhenin_at_ifr88.cnrs-mrs.fr)
Date: Thu Feb 17 2011 - 04:10:04 CST

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
>