From: MannyEful E (mannyeful_at_gmail.com)
Date: Sat Apr 02 2016 - 12:28:06 CDT

Thank you John and Norman!

On Thu, Mar 31, 2016 at 4:38 PM, John Stone <johns_at_ks.uiuc.edu> wrote:

> Don't use the variable name "list" because "list" is a Tcl command.
> You could rename that to "lst" or something more indicative of what
> it actually contains and the problem might go away.
> In general you should avoid using Tcl command names as variables, as
> they can cause confusing syntax errors and the like.
> That seems to be a likely source of your problem.
>
> Cheers,
> John
>
> On Thu, Mar 31, 2016 at 01:51:48AM +0100, MannyEful E wrote:
> > Hello Everyone,A
> > I have seen a few requests for the average residence time of water
> > molecules around the protein and am currently working on such a
> script. I
> > get a few errors which I believe are the result of the incorrect
> syntax.
> > Could someone please point me in the right direction?
> > The script is attached below and the two lines responsible for the
> errors
> > are highlighted with arrows. (<-----------)
> > The error message:
> > A invalid character "." in expression "...st +])/[llength $list]."A
> > In sum, the script does the following:
> > 1. For each water molecules it counts the number of separate occasions
> > that it approaches the polymer, including the times that it leaves but
> > returns (A).A
> > 2. It also records how long they stay (B).
> > 3. An array of A/B is produced (C).
> > 4. Finally, the average and stdev of C is obtained and printed. A (The
> > errors crop up here when I try to output my results in double
> precision,
> > whilst using a proc function.)
> > #################### A The Residence Script: A A
> #####################
> > puts "Bash usage: vmd -dispdev text -gro file.gro -xtc file.xtc -e
> > residence_times_v2.tcl"
> > puts "Next, type: residence nth unit pol"
> > puts "e.g.: residence 1st 500 pvc "
> > #################################################################
> > # process used to calculate the average of values in an array
> > proc lavg list {
> > expr {double([join $list +])/[llength $list].}
> > }
> > A
> > A # process used to calculate the stdev of values in an array
> > proc stddev list {
> > A set m [lavg $list]
> > A expr {sqrt([mean2 $list]-$m*$m)}
> > A }
> > A
> > #################################################################
> > A
> > #proc residence {nth pol unit} {
> > # Variables for Atomselection
> > set nth 1st; # which repeat is it
> > set unit 5; # how many units is the polymer
> > set pol pva; # what is the polymer name
> > set min1 0.0; A A # lower limit for the hydration shell1 in
> angstroms
> > set max1 10.5;A # higher limit of the hydration shell1 in angstroms
> > set fp1 [open "residence_times_tracker.xvg" w]; # tracks the
> calculation
> > outputs
> > set fp2 [open "residence_times_output.xvg" w]; # tracks the output
> > section
> > set atom_ref ${pol}${unit};A # VMD residue name of the reference
> atom
> > set sel [atomselect top all];A # select all the atomsA
> > set all [$sel num];A # count the number of atoms
> > set fn [molinfo top get numframes]; # count the number of frames
> > set stepsize 40.00;
> > #################################################################
> > # set all the histograms to zero
> > for {set i 0} {$i < $all} {incr i} {
> > set hist1($i) 0; A # counter: If it is within region
> > set hist2($i) 0;A A # note: If it left
> > set hist3($i) 0; A # counter: Sum of those that leftA
> > set hist4($i) 0; # counter
> > }
> > #################################################################
> > # loop through all the frames
> > for {set i 0} {$i < $fn} {incr i} {
> > # do the following for all the frames except the last one
> > if {$i < [expr $fn - 1]} {
> > # go to the next frame
> > animate goto $i;
> > # outputs frame we are on
> > puts $fp1 "Frame $i of $fn: [expr {double(($i * $stepsize)/1000)}] ns
> ";A
> > # make selections in frame and in next frame
> > set sel1 [atomselect top "(name OW and same resid as (resname SOL and
> > (within $max1 of resname $atom_ref))) and not (name OW and same resid
> as
> > (resname SOL and (within $min1 of resname $atom_ref)))" frame $i];A
> > set sel2 [atomselect top "(name OW and same resid as (resname SOL and
> > (within $max1 of resname $atom_ref))) and not (name OW and same resid
> as
> > (resname SOL and (within $min1 of resname $atom_ref)))" frame [expr
> $i +
> > 1]];A
> > A A A $sel1 update; # updates dynamic water selection1 for this
> frame
> > A A A $sel2 update; # updates dynamic water selection2 for this
> frame
> > puts $fp1 "Both selections made";
> > # creates a list of atoms indices from the next frame
> > set list2 [$sel2 get index];
> > puts $fp1 "List of future atom indices made: \n $list2";
> > # for every water oxygen atom that is near the polymer
> > # get the atoms index
> > foreach present [$sel1 get index] {
> > puts $fp1 "ID: $present";
> > # check: Is it near the polymer (Always true in this loop)
> > incr hist1([expr $present - 1]);A # add 1 if present near the area
> of
> > interest
> > puts $fp1 "HIST1: $hist1([expr $present - 1])";
> > # check: Is it coming back?
> > if { $hist2([expr $present - 1]) == 1 } A {
> > incr hist3([expr $present - 1]);A # add 1 if previous hist2 was 1
> > because it left and came back.A
> > puts $fp1 "HIST3: $hist3([expr $present - 1])";
> > }A
> > # check: Will it leave?
> > # compare the atom indices to the list of selection2
> > # and if it is present, it will return a value greater than 1
> > # so add one to our approaches counter
> > # if it has left it returns a -1
> > if { [lsearch -exact $list2 $present] >= 0} {
> > set hist2([expr $present - 1]) 1;
> > incr hist4([expr $present - 1]);
> > puts $fp1 "It will stay in the next frame";A
> > } else {
> > set hist2([expr $present - 1]) 0;
> > puts $fp1 "It will leave in the next frame";
> > }
> > }
> > } else {
> > # go to the next frame
> > animate goto $i;
> > # outputs frame we are on
> > puts $fp1 "Frame $i of $fn: [expr {double(($i * $stepsize)/1000)}] ns
> ";A
> > # make selections in frame and in next frame
> > set sel1 [atomselect top "(name OW and same resid as (resname SOL and
> > (within $max1 of resname $atom_ref))) and not (name OW and same resid
> as
> > (resname SOL and (within $min1 of resname $atom_ref)))" frame $i];A
> > foreach present [$sel1 get index] {
> > puts $fp1 "ID: $present";
> > #Any present? Always in this loop
> > incr hist1([expr $present - 1]); # add 1 because it is present near
> area
> > of interest
> > # check: Is it coming back?
> > if { $hist2([expr $present - 1]) == 1 } A {
> > incr hist3([expr $present - 1]);A # add 1 if previous hist2 was 1
> > because it left and came back.A
> > puts $fp1 "HIST3: $hist3([expr $present - 1])";
> > }A
> > # in order to account for those occasions where the water molecule
> > remained at the last frame
> > # artificially set all atoms present in the area of interest to leave
> in
> > the last frameA
> > set hist2([expr $present - 1]) 1;
> > incr hist4([expr $present - 1]);
> > puts $fp1 "It will stay in the next frame";A
> > }
> > }
> > $sel1 delete; array unset sel1; # unset the selection1A
> > $sel2 delete; array unset sel2; # unset the selection2A
> > }
> > #################################################################
> > # set new constants to zero
> > set hist_final {};A
> > set counter 0;
> > set comebacks 0;
> > # loop through all atoms (non-oxygens as well)
> > for {set i 0} {$i < [llength $hist1([expr $all - 1])]} {incr i} {
> > puts $fp2 "Atom ${i} : [lindex $hist1([expr $all - 1]) $i]";
> > # find the atoms that have residence times (ie if the hist1 counter is
> > greater than zero, a residence time is available)A
> > if {[lindex $hist1([expr $all - 1]) $i] > 0} {
> > # count the number of atoms with residences recorded
> > incr $counter;
> > puts $fp2 "Number of atoms with residences: $counter Counter";
> > # calculate the individual residence times for each oxygen atom that
> was
> > near the molecule
> > puts $fp2 " [lindex $hist1([expr $all - 1]) $i] / [lindex $hist4([expr
> > $all - 1]) $i] = [expr [lindex $hist1([expr $all - 1]) $i] / [lindex
> > $hist4([expr $all - 1]) $i]]";
> > # populate the array hist_final
> > lappend hist_final [expr {double([lindex $hist1([expr $all - 1]) $i])
> /
> > double([lindex $hist4([expr $all - 1]) $i])}];
> > # identify how many times the water molecules returned after leaving
> > set comebacks [expr $comebacks + [lindex $hist3([expr $all - 1]) $i]];
> > } else {
> > puts $fp2 "No residence times to record ${i} - number [lindex
> $hist1([expr
> > $all - 1]) ${i}]";
> > }
> > }
> > puts $fp2 "$counter : Number of oxygen atoms with residences within
> $max1
> > of the polymer over $fn frames";
> > puts $fp2 "$comebacks : Number of times any of the $counter water
> > molecules";
> > puts $fp2 "[expr {[lavg (double($hist_final))] * $stepsize}] ps :
> Average
> > residence time "; A <--------------
> > puts $fp2 "[expr {[stddev (double($hist_final))] * $stepsize}] ps :
> Stdev
> > "; A <-----------
> > close $fp1;
> > close $fp2;
> > #################################################################
> > # unset all the selections for faster execution
> > array unset nth;
> > array unset unit;
> > array unset pol;
> > array unset min1;A
> > array unset max1;A
> > array unset fp1;
> > array unset fp2;
> > array unset atom_ref;A
> > array unset fn;A
> > array unset stepsize;
> > $sel delete; array unset sel;
> > array unset list2;A
> > array unset present;
> > for {set i 0} {$i < $all} {incr i} {
> > array unset hist1($i); A
> > array unset hist2($i);A A
> > array unset hist3($i); A A
> > array unset hist4($i);
> > }
> > array unset all;A
> > array unset hist_final;A
> > array unset counter;
> > array unset comebacks;
> > #################################################################
> > #}
> > END
>
> --
> NIH Center for Macromolecular Modeling and Bioinformatics
> Beckman Institute for Advanced Science and Technology
> University of Illinois, 405 N. Mathews Ave, Urbana, IL 61801
> http://www.ks.uiuc.edu/~johns/ Phone: 217-244-3349
> http://www.ks.uiuc.edu/Research/vmd/
>