From: Ajasja Ljubetič (ajasja.ljubetic_at_gmail.com)
Date: Tue Nov 18 2014 - 08:24:59 CST

Probably not the answer you are looking for, but my advice would be to
write the number of waters vs timestep to a text file and do the analysis
in Python (or R or Mathematica or whatever specialized software of your
choice).
When I started doing that, life became much easier:)
Unless of course you need to display the histogram interactively etc..

Best regards,
Ajasja

On 18 November 2014 14:40, MannyEful E <mannyeful_at_gmail.com> wrote:

> 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)”;
>
> }
>
>