VMD-L Mailing List

From: J. Nathan Scott (scottjn_at_chemistry.montana.edu)
Date: Wed Feb 16 2011 - 17:17:28 CST

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