• ## Outreach

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.

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