From: MannyEful E (mannyeful_at_gmail.com)
Date: Thu Jan 22 2015 - 23:18:10 CST

Hello Everyone,

I wrote a script to calculate the tetrahedral orientation order parameter
of liquid water over time. Could anyone advise me on what causes VMD to be
killed despite the use of "$sel delete" / "unset sel" / "array unset sel"?
I imagine that this is a memory leak problem, however I can't spot it.

Thanks in advanced for your time and help!

Script:
##################################################################################
#: Title : q4_calc.tcl

      #
#: Description : Get the tetrahedral orientation parameter for waters
                                                    #
#: Usage: vmd -dispdev text -gro test.gro -e q4_calc_bulk.tcl
                                                       #
##################################################################################

proc angle {Di Hi Ai} {

# Calculates angles between three points
# select atoms
set selD [atomselect top "resid $Di and name OW"]; # select atom1
set selH [atomselect top "index $Hi and name OW"]; # select atom2
set selA [atomselect top "resid $Ai and name OW"]; # select atom3

# obtain corresponding coordinates
set D [lindex [$selD get {x y z}] 0];
set H [lindex [$selH get {x y z}] 0];
set A [lindex [$selA get {x y z}] 0];

# obtain cosine of the angle between the 3 selected atoms
 global M_PI; # get pi value
 set hd [vecsub $D $H]; # calculate hd vector
 set ha [vecsub $A $H]; # calculate ha vector
 set cosine [expr [vecdot $hd $ha] / ( [veclength $hd] * [veclength
$ha])]; # calculate cosine

# return [expr acos($cosine)*(180.0/$M_PI)]; # convert cosine to angle
in degrees
  return $cosine

# delete variables we do not need any more
$Di delete;
$Hi delete;
$Ai delete;
$selD delete;
$selH delete;
$selA delete;
$D delete;
$H delete;
$A delete;
$hd delete;
$ha delete;

}

proc nwat {n sel} {
  # Calculates the n nearest neighbours
  # reset the l3 variables to become an empty array
  set l3 {};

  # find the oxygens within 10 angstroms of the selection
  set lists [measure contacts 10 [atomselect top "name OW"] $sel];

  # obtains the list of oxygens indices
  set wlist [lindex $lists 0];

  # obtain the list of the selection atoms indices
  set slist [lindex $lists 1];

  # for every listed selection and its corresponding oxygen
  # get the residue number for the oxygen
  # get the distance between the two
  # appends this distance and corresponding resids to a list called l1
  foreach satom $slist watom $wlist {
    set wresnum [[atomselect top "index $watom"] get resid];
    set d [measure bond "$watom $satom"];
    lappend l1 [list $wresnum $d];
  }

  # make a variable called l1 which sorts out the 4 closest distances
  # appends the 4 atom resids to a list called l3
  set l2 [lrange [lsort -index 1 [lsort -unique -index 0 [lsort -index 1
-decreasing $l1]]] 0 [expr $n - 1]];
  #puts " F ";
  foreach water $l2 {
    lappend l3 "[lindex $water 0]";
  }

  # returns a list of resids for the 4 nearest
  return $l3;

# delete variables we do not need any more
array unset $selD;
array unset $selA;
array unset $selH;
array unset $D;
array unset $A;
array unset $H;
$Di delete; unset Di;
$Hi delete; unset Hi;
$Ai delete; unset Ai;
$selD delete; array unset selD;
$selH delete;array unset selH;
$selA delete;array unset selA;
$D delete;array unset D;
$H delete;array unset A;
$A delete; array unset A;
$hd delete;array unset hd;
$ha delete;array unset ha;

}

# assign a filename
## set q_file [open "q_19Jan2015_sla_260K_0-3_8.txt" w];

# define the number of frames
set end 3;

# make a selection for which q will be calculated
set total_wat [atomselect top "name OW and (resname SLA) "];

# loop through all the frames
for { set i 0 } { $i < $end } { incr i} {

       # go to the specified frame
       animate goto $i;

       # selection refreshed for the new frame
       $total_wat frame $i;

       # counts the number of atoms selected
       set num_Totalwat [$total_wat num];

       # get the residue indices and atom indices for selection
       set wat_resid [$total_wat get resid];
       set wat_index [$total_wat get index];

       # sets a counter to zero
       set checkingNums 0;

       # for each water we do the following
       foreach waterID $wat_index waterRID $wat_resid {

       # add 1 to the counter and print out so we know how far into the
loop we have gone
       incr checkingNums;
       puts " WATER $checkingNums OF $num_Totalwat ";

       # make a temporary selection using the atom indices and residue
indices
       set wat_coord [atomselect top "name OW and index $waterID and resid
$waterRID"];

       # find the 4 closest waters to this atom
       set close_wat [nwat 4 $wat_coord];

      # Find all the possible 3 angle combinations which can be made with
the 4 neighbour whilst wat_coord is at the centre
      # set the sum to zero
       set sum 0.0;

       # set counter to zero
       set countAngles 0;
       for { set k 0 } { $k < 4 } { incr k } {
           for { set l [expr $k+1] } { $l < 4 } { incr l } {

               # add 1 to angles which have been calculated. The final
value should always be 6 for 4 neighbours.
               incr countAngles;

               # select the first atom index
               set hi1 "[lindex $close_wat $k]";

               # select the reference atom index
               set atom2 [[atomselect top "name OW and (index $waterID and
resid $waterRID)"] get {index}];

               # select the third atom index
               set hi3 "[lindex $close_wat $l]";

               # calculate the angle usin proc calcangle defined earlier
               set calcangle [angle $hi1 $atom2 $hi3];

               # add up cosine results to those from previous angles
               set sum [expr
{(($calcangle+1.0/3.0)*($calcangle+1.0/3.0))+$sum}];

           }
       }

      # define the q value
       set q [expr {1.0-((3.0/8.0)*$sum)}];

       # print the Frame Loop number Atom index Q
value details to file
       puts [format "Frame: $i\t ID: $checkingNums \t AtomIDX: $waterID \t
Q: %.2f" $q];

## puts $q_file "ID: $checkingNums AtomIDX: $waterID Q: $q"; # print
details to file

     }

      # delete variables we do not need any more
      $wat_coord delete; array unset wat_coord;
      $close_wat delete; array unset close_wat;
      $k delete; array unset k;
      $l delete; array unset l;
      $hi1 delete; array unset hi1;
      $atom2 delete; array unset atom2;
      $hi3 delete; array unset hi3;
      $calcangle delete; array unset calcangle;
      $sum delete; array unset sum;
      $q delete; array unset q;

## puts $q_file "END"; # add END line after every frame

# delete variables we do not need any more
 $total_wat delete; array unset total_wat;
 $num_Totalwat delete; array unset num_Totalwat;
 $wat_resid $delete; array unset wat_resid;
 $wat_index delete; array unset wat_index;
 $waterID delete; array unset waterID;
 $waterRID delete; array unset waterRID;
 $checkingNums delete; array unset checkingNums;
 $countAngles delete; array unset countAngles;

   }

## puts $q_file " "; # add empty newline to file

## close $q_file; # close file

# delete variables we do not need any more
$total_wat delete;

END