VMD-L Mailing List

From: MannyEful E (mannyeful_at_gmail.com)
Date: Tue Nov 18 2014 - 07:40:29 CST

I need some advice with VMD and tcl. I am counting the number of waters
around a shell of a polymer and then I want to create a histogram based on
these values.
I have currently obtained the value for each frame, but I’m stuck on how to
increment the histogram array ($hist($i)) correctly.
I can’t quite put my finger on what needs to be changed in the script below
(see arrows). Any ideas??

This part of the script will be important later to help me calculate the
orientational order parameter of surrounding water [q=1-3/8
SUM_OVER_6_ANGLES (cos (phi_ijk) + 1/3) S. Lee et al (2005)]

Thank you for your time.

# Variables for Atomselection
set min 0 ; # lower limit of the hydration shell
set max 3.5 ; # higher limit of the hydration shell
set atom_ref pva5; # VMD residue name of the reference atom

# Variables for Calculations
set frame_start 0;
set num_steps [molinfo top get numframes];
set fp [open "g_traj_hydration_vmd_3.txt" w];
set frame_end $num_steps;
set stride 1;
set suma 0;

# Variables for Histogram
set min_val 0;
set max_val 100;
set nbins 100;
set binWidth [expr (($max_val - $min_val)/$nbins)];

# initialises an array befor counting by setting all bars to zeroes
for {set i 0} {$i < $nbins} {incr i} {
set hist($i) 0
}

#script runs calculates the number of oxygen atoms within 3.5 of the polymer
for {set frame $frame_start} {$frame < $frame_end} {incr frame} {
 # goes to the frame and updates selections
puts "Frame $frame";
set oxygens [atomselect top "(name OW) and (within $max of resname
$atom_ref)"];
$oxygens frame $frame;
        $oxygens update;

# number of neighbouring waters, excludes itself
set num_oxygens [$oxygens num] ;
puts "number of oxygens is: $num_oxygens”;

# sum total number of waters that have been near the polymer
set suma [expr $suma + $num_oxygens];
puts "current total is now: $suma”;

# identifies which bin to increment
set whichBin [expr {int ($num_oxygens - $min_val)/$binWidth}];
puts "whichBin is: $whichBin”;

set $hist($num_oxygens) [expr $hist($num_oxygens) + 1]; <------

# testing print out during the loop. I get no change here <-------
puts "the whichBin Bin is incremented to: $hist($num_oxygens)";

}

# print out the histogram
for {set i 0} {$i < $nbins} {incr i} {
puts "$i $hist($i)”;

}