VMD-L Mailing List

From: John Stone (johns_at_ks.uiuc.edu)
Date: Fri Jan 23 2015 - 23:25:25 CST

Hi,
  In your nwat proc, you're overwriting the contents of selection2
before you delete the previously created selection, you need to delete
it before you overwrite it, so you would do something like this:

foreach satom $slist watom $wlist {
  set selection2 [atomselect top "index $watom"];
  ... rest of the loop ...
  $selection2 delete
}

Cheers,
  John Stone
  vmd_at_ks.uiuc.edu

On Sat, Jan 24, 2015 at 05:16:27AM +0000, MannyEful E wrote:
> Hi John,
>
> The script runs much further than it ever has! This time round it comes to
> a halt during the last frame (ie the 3rd frame). Do you happen to have any
> more suggestions??
>
> I have also tried to unset selections as soon as they are finished being
> used. This speeds up the script but still kills VMD during the third
> frame.
>
> Script (without the comments):
> proc angle {Di Hi Ai} {
>
> set selD [atomselect top "resid $Di and name OW"];
> set selH [atomselect top "index $Hi and name OW"];
> set selA [atomselect top "resid $Ai and name OW"];
> A
> 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];
>
> global M_PI;
>
> set hd [vecsub $D $H];
>
> set ha [vecsub $A $H];A
>
> set cosine [expr [vecdot $hd $ha] / ( [veclength $hd] * [veclength
> $ha])];A
>
> # delete selections we dont need any moreA
> $selD delete; array unset selD;
> $selH delete;array unset selH;
> $selA delete;array unset selA;
> array unset D;
> array unset H;
> array unset A;
> array unset hd;
> array unset ha;
> unset Di;
> unset Hi;
> unset Ai;
>
> return $cosine;
> A
> }
>
> proc nwat {n sel} {
>
> A set l3 {};
> A set selection1 [atomselect top "name OW"];
>
> A set lists [measure contacts 10 $selection1 $sel];
>
> A set wlist [lindex $lists 0];
> A set slist [lindex $lists 1];
>
> A foreach satom $slist watom $wlist {
> A A A set selection2 [atomselect top "index $watom"];
> A A A set wresnum [$selection2 get resid];
> A A A set d [measure bond "$watom $satom"];
> A A A lappend l1 [list $wresnum $d];
> A }
>
> A set l2 [lrange [lsort -index 1 [lsort -unique -index 0 [lsort -index 1
> -decreasing $l1]]] 0 [expr $n - 1]];
>
> A foreach water $l2 {
> A A A lappend l3 "[lindex $water 0]";
> A }
>
> # delete selections we dont need any moreA
> A $selection1 delete;
> A $selection2 delete;
> A array unset n;
> A array unset $sel;
> A array unset lists;
> A array unset wlist;
> A array unset slist;
> A array unset satom;
> A array unset watom;
> A array unset wresnum;
> A array unset d;
> A array unset l1;
> A array unset l2;
>
> A return $l3;
>
> }
>
> A
> ##set q_file [open "q_19Jan2015_sla_260K_0-3_8.txt" w];
> A
> set end 3;
> A
> set total_wat [atomselect top "name OW and (resname SLA) "];A
> A A
> #rotate through all the framesA
> for { set i 0 } { $i < $end } { incr i} {A
> A
> A A A A A A animate goto $i;
>
> A A A A A A $total_wat frame $i;
>
> A A A A A A set num_Totalwat [$total_wat num];
> A A A A A
> A A A A A A set wat_resid [$total_wat get resid];
> A A A A A A set wat_index [$total_wat get index];
> A
> A A A A A A set checkingNums 0;
>
> A A A A A A foreach waterID $wat_index waterRID $wat_resid {
>
> A A A A A A incr checkingNums;
>
> A A A A A A set wat_coord [atomselect top "name OW and index $waterID and
> resid $waterRID"];A
> A A A A A A A A A
> A A A A A A set close_wat [nwat 4 $wat_coord];
>
> A A A A A A set sum 0.0;
>
> A A A A A A set countAngles 0;
>
> A A A A A A for { set k 0 } { $k < 4 } { incr k } {A
> A A A A A A A A A A for { set l [expr $k+1] } { $l < 4 } { incr l } {A
>
> A A A A A A A A A A A A A A incr countAngles;
>
> A A A A A A A A A A A A A A set hi1 "[lindex $close_wat $k]";
>
> A A A A A A A A A A A A A A set atom2 [$wat_coord get {index}];
> A
> A A A A A A A A A A A A A A set hi3 "[lindex $close_wat $l]";A
>
> A A A A A A A A A A A A A A set calcangle [angle $hi1 $atom2 $hi3];
> A
> A A A A A A A A A A A A A A set sum [expr
> {(($calcangle+1.0/3.0)*($calcangle+1.0/3.0))+$sum}];
> A
> A A A A A A A A A A }A
> A A A A A A }
> A
> A
> A A A A A A set q [expr {1.0-((3.0/8.0)*$sum)}];
>
> A A A A A A puts [format "Frame: $i\t ID: $checkingNums of $num_Totalwat
> "];
> A
> ##A A puts $q_file "ID: $checkingNumsA A A AtomIDX: $waterIDA A Q:
> $q";A
>
> A
> A A A A }
> A A A A
> # delete selections we dont need any moreA
> A A A A A $wat_coord delete; array unset wat_coord;
> A A A A A array unset close_wat;
> A A A A A array unset k;
> A A A A A array unset l;
> A A A A A array unset hi1;
> A A A A A array unset atom2;
> A A A A A array unset hi3;
> A A A A A array unset calcangle;
> A A A A A array unset sum;
> A A A A A array unset q;
>
> ##A A puts $q_file "ENDMDL";
>
> # delete selections we dont need any moreA
> A array unset total_wat;
> A array unset num_Totalwat;
> A array unset wat_resid;
> A array unset wat_index;
> A array unset waterID;
> A array unset waterRID;
> A array unset checkingNums;
> A array unset countAngles;
> A
> A A }A
>
> ##A A puts $q_file " ";
>
> ##close $q_file;
>
> $total_wat delete;
> A
> A
> END
>
> On Fri, Jan 23, 2015 at 5:39 PM, John Stone <[1]johns_at_ks.uiuc.edu> wrote:
>
> Hi,
> A Okay, the next problem you have is that in several places in your
> scripts you're creating atom selections which are immediately "leaked"
> because you don't assign them to a variable and delete them:
> A A set lists [measure contacts 10 [atomselect top "name OW"] $sel];
>
> So, for the line above, you need to rewrite it more like this:
> A A set owsel [atomselect top "name OS"]
> A A set lists [measure contacts 10 $owsel $sel];
> A A $owsel delete
>
> I saw several places in your code that have this sort of problem.
>
> Cheers,
> A John Stone
> A [2]vmd_at_ks.uiuc.edu
> On Fri, Jan 23, 2015 at 11:46:57AM +0000, MannyEful E wrote:
> >A A Thank you John,
> >
> >A A I have made the changes you've suggested (the areas labelled with
> >A A arrows).AA And have also tried the same for the nwat process for
> good
> >A A measure.
> >
> >A A Although the script runs faster, vmd still crashes.
> >
> >A A Script:
> >A A
> ##################################################################################
> >A A #: TitleA A A A A AA : q4_calc.tclA A AA A A AA A A AA A A
> AA A A A
> >A A A A AA A A AA A AA A A AA A A AA A A AA A A AA A A AA A A
> AA A AA A A A
> >A A A A AA A A AA A A AA A A AA A A AA A AA A A AA A A AA A A
> AA A A AA A A A
> >A A A A AA A AA A A AA A A AA A A AA AA #
> >A A #: Description : Get the tetrahedral orientation parameter for
> >A A watersA A AA A A AA AA AA A A AA A A AA A A AA A A AA A A
> AA A A A AA A
> >A A A A AA A A AA A A AA A A A AA #
> >A A #: Usage:A A A A A A A A A A AA vmd -dispdev text -gro test.gro
> -e
> >A A q4_calc_bulk.tcl A A AA A A AA A A AA A A AA A A AA A A AA
> A A AA A A A
> >A A A A AA AA AA A A AA A A AA A A AA A AA #
> >A A
> ##################################################################################
> >
> >A A proc angle {Di Hi Ai} {
> >
> >A A # Calculates angles between three points
> >A A # select atoms
> >A A set selD [atomselect top "resid $Di and name OW"];A A A A AA #
> select
> >A A atom1
> >A A set selH [atomselect top "index $Hi and name OW"];A A A AA #
> select atom2
> >A A set selA [atomselect top "resid $Ai and name OW"];A A A A AA #
> select
> >A A atom3
> >A A A
> >A A # obtain corresponding coordinates
> >A A set D [lindex [$selD get {x y z}] 0];
> >A A set H [lindex [$selH get {x y z}] 0];
> >A A set A [lindex [$selA get {x y z}] 0];
> >
> >A A # obtain cosine of the angle between the 3 selected atoms
> >A A A global M_PI;A A A A A A A A A A A A A A A A A A A AA # get pi
> value
> >A A A set hd [vecsub $D $H];A A A AA # calculate hd vector
> >A A A set ha [vecsub $A $H];A A A AA # calculate ha vector
> >A A A set cosine [expr [vecdot $hd $ha] / ( [veclength $hd] *
> [veclength
> >A A $ha])];AA AA # calculate cosine
> >
> >A A # delete variables we do not need any moreAA <----
> >A A A $selD delete; array unset selD;
> >A A A $selH delete;array unset selH;
> >A A A $selA delete;array unset selA;
> >A A A array unset D;
> >A A A array unset H;
> >A A A array unset A;
> >A A A array unset hd;
> >A A A array unset ha;
> >A A A unset Di;
> >A A A unset Hi;
> >A A A unset Ai;
> >
> >A A # return [expr acos($cosine)*(180.0/$M_PI)]; A AA # convert
> cosine to
> >A A angle in degreesAA A
> >A A AA return $cosine
> >A A A
> >A A }
> >
> >A A proc nwat {n sel} {
> >A A AA # Calculates the n nearest neighbours
> >A A AA # reset the l3 variables to become an empty array
> >A A AA set l3 {};
> >
> >A A AA # find the oxygens within 10 angstroms of the selection
> >A A AA set lists [measure contacts 10 [atomselect top "name OW"]
> $sel];
> >
> >A A AA # obtains the list of oxygens indices
> >A A AA set wlist [lindex $lists 0];
> >
> >A A AA # obtain the list of the selection atoms indices
> >A A AA set slist [lindex $lists 1];
> >
> >A A AA # for every listed selection and its corresponding oxygen
> >A A AA # get the residue number for the oxygen
> >A A AA # get the distance between the twoA
> >A A AA # appends this distance and corresponding resids to a list
> called l1
> >A A AA foreach satom $slist watom $wlist {
> >A A A A AA set wresnum [[atomselect top "index $watom"] get resid];
> >A A A A AA set d [measure bond "$watom $satom"];
> >A A A A AA lappend l1 [list $wresnum $d];
> >A A AA }
> >
> >A A AA # make a variable called l1 which sorts out the 4 closest
> distancesA
> >A A AA # appends the 4 atom resids to a list called l3
> >A A AA set l2 [lrange [lsort -index 1 [lsort -unique -index 0 [lsort
> -index 1
> >A A -decreasing $l1]]] 0 [expr $n - 1]];
> >A A AA #puts " F ";A
> >A A AA foreach water $l2 {
> >A A A A AA lappend l3 "[lindex $water 0]";
> >A A AA }
> >
> >A A # delete variables we do not need any moreAA AA <----
> >A A A array unset n;
> >A A A array unset $sel;
> >A A A array unset lists;
> >A A A array unset wlist;
> >A A A array unset slist;
> >A A A array unset satom;
> >A A A array unset watom;
> >A A A array unset wresnum;
> >A A A array unset d;
> >A A A array unset l1;
> >A A A array unset l2;
> >
> >A A AA # returns a list of resids for the 4 nearest
> >A A AA return $l3;
> >
> >A A }
> >
> >A A # assign a filename
> >A A ##AA set q_file [open "q_19Jan2015_sla_260K_0-3_8.txt" w];
> >
> >A A # define the number of frames
> >A A set end 3;
> >A A A
> >A A # make a selection for which q will be calculated
> >A A set total_wat [atomselect top "name OW and (resname SLA) "];A
> >A A A A
> >A A # loop through all the framesA
> >A A for { set i 0 } { $i < $end } { incr i} {A
> >A A A
> >A A A A A A A AA # go to the specified frame
> >A A A A A A A AA animate goto $i;
> >A A A A A A A
> >A A A A A A A AA # selection refreshed for the new frame
> >A A A A A A A AA $total_wat frame $i;
> >
> >A A A A A A A AA # counts the number of atoms selected
> >A A A A A A A AA set num_Totalwat [$total_wat num];
> >A A A A AA A A
> >A A A A A A A AA # get the residue indices and atom indices for
> selection
> >A A A A A A A AA set wat_resid [$total_wat get resid];
> >A A A A A A A AA set wat_index [$total_wat get index];
> >
> >A A A A A A A AA # sets a counter to zero
> >A A A A A A A AA set checkingNums 0;
> >
> >A A A A A A A AA # for each water we do the following
> >A A A A A A A AA foreach waterID $wat_index waterRID $wat_resid {
> >
> >A A A A A A A AA # add 1 to the counter and print out so we know how
> far into
> >A A the loop we have gone
> >A A A A A A A AA incr checkingNums;
> >A A A A A A A AA puts " WATER $checkingNums OF $num_Totalwat ";A
> >A A A A A A A A A A
> >A A A A A A A AA # make a temporary selection using the atom indices
> and
> >A A residue indicesAA A
> >A A A A A A A AA set wat_coord [atomselect top "name OW and index
> $waterID and
> >A A resid $waterRID"];A
> >A A A
> >A A A A A A A AA # find the 4 closest waters to this atom
> >A A A A A A A AA set close_wat [nwat 4 $wat_coord];
> >
> >A A A A A A AA # Find all the possible 3 angle combinations which
> can be made
> >A A with the 4 neighbour whilst wat_coord is at the centre
> >A A AA AA AA # set the sum to zero A A A A A A A A A A A A A A A A
> A A A A A
> >A A A A A A A AA set sum 0.0;
> >
> >A A A A A A A AA # set counter to zero
> >A A A A A A A AA set countAngles 0;
> >A A A A A A A AA for { set k 0 } { $k < 4 } { incr k } {A
> >A A A A A A A A A A A AA for { set l [expr $k+1] } { $l < 4 } { incr
> l } {A
> >
> >A A A A A A A A A A A A A A A AA # add 1 to angles which have been
> calculated.
> >A A The final value should always be 6 for 4 neighbours.
> >A A A A A A A A A A A A A A A AA incr countAngles;
> >
> >A A A A A A A A A A A A A A A AA # select the first atom index
> >A A A A A A A A A A A A A A A AA set hi1 "[lindex $close_wat $k]";
> >
> >A A A A A A A A A A A A A A A AA # select the reference atom index
> >A A A A A A A A A A A A A A A AA set atom2 [[atomselect top "name OW
> and
> >A A (index $waterID and resid $waterRID)"] get {index}];
> >
> >A A A A A A A A A A A A A A A AA # select the third atom indexA
> >A A A A A A A A A A A A A A A AA set hi3 "[lindex $close_wat $l]";A
> >
> >A A A A A A A A A A A A A A A AA # calculate the angle usin proc
> calcangle
> >A A defined earlier
> >A A A A A A A A A A A A A A A AA set calcangle [angle $hi1 $atom2
> $hi3];
> >A A A
> >A A A A A A A A A A A A A A A AA # add up cosine results to those
> from
> >A A previous angles
> >A A A A A A A A A A A A A A A AA set sum [expr
> >A A {(($calcangle+1.0/3.0)*($calcangle+1.0/3.0))+$sum}];
> >A A A
> >A A A A A A A A A A A AA }A
> >A A A A A A A AA }
> >A A A
> >A A A A A A AA # define the q value
> >A A A A A A A AA set q [expr {1.0-((3.0/8.0)*$sum)}];
> >
> >A A A A A A A AA # print theA A A AA FrameA A A A A A AA Loop
> >A A numberA A A A A A AA Atom index A A A A A A AA Q valueA A A AA
> details to
> >A A file
> >A A A A A A A AA puts [format "Frame: $i\t ID: $checkingNums \t
> AtomIDX:
> >A A $waterID \t Q: %.2f"A AA $q];
> >A A A
> >A A ##AA puts $q_file "ID: $checkingNumsA A AA AtomIDX: $waterIDA
> AA Q: $q";A
> >A A # print details to file
> >A A A
> >A A A A A AA }
> >
> >A A A A A A AA # delete variables we do not need any moreAA AA A A
> A A
> >A A A A A A AA $wat_coord delete; array unset wat_coord;
> >A A A A A A AA $close_wat delete; array unset close_wat;
> >A A A A A A AA $k delete; array unset k;
> >A A A A A A AA $l delete; array unset l;
> >A A A A A A AA $hi1 delete; array unset hi1;
> >A A A A A A AA $atom2 delete; array unset atom2;
> >A A A A A A AA $hi3 delete; array unset hi3;
> >A A A A A A AA $calcangle delete; array unset calcangle;
> >A A A A A A AA $sum delete; array unset sum;
> >A A A A A A AA $q delete; array unset q;
> >A A A A A A A A
> >
> >A A ##A AA puts $q_file "END"; # add END line after every frame
> >
> >A A # delete variables we do not need any moreAA A
> >A A A $total_wat delete; array unset total_wat;
> >A A A $num_Totalwat delete; array unset num_Totalwat;
> >A A A $wat_resid $delete; array unset wat_resid;
> >A A A $wat_index delete; array unset wat_index;
> >A A A $waterID delete; array unset waterID;
> >A A A $waterRID delete; array unset waterRID;
> >A A A $checkingNums delete; array unset checkingNums;
> >A A A $countAngles delete; array unset countAngles;
> >A A A
> >A A A AA }A
> >
> >A A ##A AA puts $q_file " "; # add empty newline to file
> >
> >A A ##A AA close $q_file;A A AA # close file
> >
> >A A # delete variables we do not need any moreAA A
> >A A $total_wat delete;
> >A A A
> >A A A
> >
> >A A END
> >A A On Fri, Jan 23, 2015 at 6:20 AM, John Stone
> <[1][3]johns_at_ks.uiuc.edu> wrote:
> >
> >A A A Hi,
> >A A A AA Your return statement in the angle proc occurs before you
> >A A A call the $sel delete commands, so they are not being executed
> at all.
> >A A A Move the return statement so it is the last line of the proc
> and that
> >A A A should cure your problem assuming there aren't other issues
> that I
> >A A A didn't see.
> >
> >A A A Cheers,
> >A A A AA John Stone
> >A A A AA [2][4]vmd_at_ks.uiuc.edu
> >A A A On Fri, Jan 23, 2015 at 05:18:10AM +0000, MannyEful E wrote:
> >A A A >AA AA Hello Everyone,
> >A A A >
> >A A A >AA AA I wrote a script to calculate the tetrahedral
> orientation order
> >A A A parameter
> >A A A >AA AA of liquid water over time. Could anyone advise me on
> what causes
> >A A A VMD to be
> >A A A >AA AA killed despite the use of "$sel delete" / "unset sel"
> / "array
> >A A A unset sel"?
> >A A A >AA AA I imagine that this is a memory leak problem, however
> I can't
> >A A A spot it.
> >A A A >
> >A A A >AA AA Thanks in advanced for your time and help!
> >A A A >
> >A A A >AA AA Script:
> >A A A >AA A
> >A A A
> ##################################################################################
> >A A A >AA AA #: TitleA A A A A AAA : q4_calc.tclA A AAA A A AAA
> A A AAA A A
> >A A A AAA A A A
> >A A A >AA AA A A AAA A A AAA A AAA A A AAA A A AAA A A AAA A
> A AAA A A AAA A A
> >A A A AAA A AAA A A A
> >A A A >AA AA A A AAA A A AAA A A AAA A A AAA A A AAA A AAA A
> A AAA A A AAA A A
> >A A A AAA A A AAA A A A
> >A A A >AA AA A A AAA A AAA A A AAA A A AAA A A AAA AAA #
> >A A A >AA AA #: Description : Get the tetrahedral orientation
> parameter for
> >A A A >AA AA watersA A AAA A A AAA AAA AAA A A AAA A A AAA A
> A AAA A A AAA A A
> >A A A AAA A A A AAA A
> >A A A >AA AA A A AAA A A AAA A A AAA A A A AAA #
> >A A A >AA AA #: Usage:A A A A A A A A A A AAA vmd -dispdev text
> -gro test.gro
> >A A A -e
> >A A A >AA AA q4_calc_bulk.tcl A A AAA A A AAA A A AAA A A AAA
> A A AAA A A AA
> >A A A A A AAA A A A
> >A A A >AA AA A A AAA AAA AAA A A AAA A A AAA A A AAA A AAA
> #
> >A A A >AA A
> >A A A
> ##################################################################################
> >A A A >
> >A A A >AA AA proc angle {Di Hi Ai} {
> >A A A >
> >A A A >AA AA # Calculates angles between three points
> >A A A >AA AA # select atoms
> >A A A >AA AA set selD [atomselect top "resid $Di and name OW"];A A
> A A AAA #
> >A A A select
> >A A A >AA AA atom1
> >A A A >AA AA set selH [atomselect top "index $Hi and name OW"];A A
> A AAA #
> >A A A select atom2
> >A A A >AA AA set selA [atomselect top "resid $Ai and name OW"];A A
> A A AAA #
> >A A A select
> >A A A >AA AA atom3
> >A A A >AA AA A
> >A A A >AA AA # obtain corresponding coordinates
> >A A A >AA AA set D [lindex [$selD get {x y z}] 0];
> >A A A >AA AA set H [lindex [$selH get {x y z}] 0];
> >A A A >AA AA set A [lindex [$selA get {x y z}] 0];
> >A A A >
> >A A A >AA AA # obtain cosine of the angle between the 3 selected
> atoms
> >A A A >AA AA A global M_PI;A A A A A A A A A A A A A A A A A A A
> AAA # get pi
> >A A A value
> >A A A >AA AA A set hd [vecsub $D $H];A A A AAA # calculate hd
> vector
> >A A A >AA AA A set ha [vecsub $A $H];A A A AAA # calculate ha
> vector
> >A A A >AA AA A set cosine [expr [vecdot $hd $ha] / ( [veclength
> $hd] *
> >A A A [veclength
> >A A A >AA AA $ha])];AAA AAA # calculate cosine
> >A A A >
> >A A A >AA AA # return [expr acos($cosine)*(180.0/$M_PI)]; A AAA #
> convert
> >A A A cosine to
> >A A A >AA AA angle in degreesAAA A
> >A A A >AA AA AAA return $cosine
> >A A A >AA AA A
> >A A A >AA AA # delete variables we do not need any moreA
> >A A A >AA AA $Di delete;
> >A A A >AA AA $Hi delete;
> >A A A >AA AA $Ai delete;
> >A A A >AA AA $selD delete;
> >A A A >AA AA $selH delete;
> >A A A >AA AA $selA delete;
> >A A A >AA AA $D delete;
> >A A A >AA AA $H delete;
> >A A A >AA AA $A delete;
> >A A A >AA AA $hd delete;
> >A A A >AA AA $ha delete;
> >A A A >AA AA A
> >A A A >AA AA }
> >A A A >
> >A A A >AA AA proc nwat {n sel} {
> >A A A >AA AA AAA # Calculates the n nearest neighbours
> >A A A >AA AA AAA # reset the l3 variables to become an empty
> array
> >A A A >AA AA AAA set l3 {};
> >A A A >
> >A A A >AA AA AAA # find the oxygens within 10 angstroms of the
> selection
> >A A A >AA AA AAA set lists [measure contacts 10 [atomselect top
> "name OW"]
> >A A A $sel];
> >A A A >
> >A A A >AA AA AAA # obtains the list of oxygens indices
> >A A A >AA AA AAA set wlist [lindex $lists 0];
> >A A A >
> >A A A >AA AA AAA # obtain the list of the selection atoms indices
> >A A A >AA AA AAA set slist [lindex $lists 1];
> >A A A >
> >A A A >AA AA AAA # for every listed selection and its
> corresponding oxygen
> >A A A >AA AA AAA # get the residue number for the oxygen
> >A A A >AA AA AAA # get the distance between the twoA
> >A A A >AA AA AAA # appends this distance and corresponding resids
> to a list
> >A A A called l1
> >A A A >AA AA AAA foreach satom $slist watom $wlist {
> >A A A >AA AA A A AAA set wresnum [[atomselect top "index $watom"]
> get resid];
> >A A A >AA AA A A AAA set d [measure bond "$watom $satom"];
> >A A A >AA AA A A AAA lappend l1 [list $wresnum $d];
> >A A A >AA AA AAA }
> >A A A >
> >A A A >AA AA AAA # make a variable called l1 which sorts out the
> 4 closest
> >A A A distancesA
> >A A A >AA AA AAA # appends the 4 atom resids to a list called l3
> >A A A >AA AA AAA set l2 [lrange [lsort -index 1 [lsort -unique
> -index 0 [lsort
> >A A A -index 1
> >A A A >AA AA -decreasing $l1]]] 0 [expr $n - 1]];
> >A A A >AA AA AAA #puts " F ";A
> >A A A >AA AA AAA foreach water $l2 {
> >A A A >AA AA A A AAA lappend l3 "[lindex $water 0]";
> >A A A >AA AA AAA }
> >A A A >
> >A A A >AA AA AAA # returns a list of resids for the 4 nearest
> >A A A >AA AA AAA return $l3;
> >A A A >
> >A A A >AA AA # delete variables we do not need any moreAAA A
> >A A A >AA AA array unset $selD;
> >A A A >AA AA array unset $selA;
> >A A A >AA AA array unset $selH;
> >A A A >AA AA array unset $D;
> >A A A >AA AA array unset $A;
> >A A A >AA AA array unset $H;
> >A A A >AA AA $Di delete; unset Di;
> >A A A >AA AA $Hi delete; unset Hi;
> >A A A >AA AA $Ai delete; unset Ai;
> >A A A >AA AA $selD delete; array unset selD;
> >A A A >AA AA $selH delete;array unset selH;
> >A A A >AA AA $selA delete;array unset selA;
> >A A A >AA AA $D delete;array unset D;
> >A A A >AA AA $H delete;array unset A;
> >A A A >AA AA $A delete; array unset A;
> >A A A >AA AA $hd delete;array unset hd;
> >A A A >AA AA $ha delete;array unset ha;
> >A A A >AA AA A
> >A A A >AA AA }
> >A A A >
> >A A A >AA AA # assign a filename
> >A A A >AA AA ##AAA set q_file [open
> "q_19Jan2015_sla_260K_0-3_8.txt" w];
> >A A A >
> >A A A >AA AA # define the number of frames
> >A A A >AA AA set end 3;
> >A A A >AA AA A
> >A A A >AA AA # make a selection for which q will be calculated
> >A A A >AA AA set total_wat [atomselect top "name OW and (resname
> SLA) "];A
> >A A A >AA AA A A
> >A A A >AA AA # loop through all the framesA
> >A A A >AA AA for { set i 0 } { $i < $end } { incr i} {A
> >A A A >AA AA A
> >A A A >AA AA A A A A A AAA # go to the specified frame
> >A A A >AA AA A A A A A AAA animate goto $i;
> >A A A >AA AA A A A A A
> >A A A >AA AA A A A A A AAA # selection refreshed for the new
> frame
> >A A A >AA AA A A A A A AAA $total_wat frame $i;
> >A A A >
> >A A A >AA AA A A A A A AAA # counts the number of atoms selected
> >A A A >AA AA A A A A A AAA set num_Totalwat [$total_wat num];
> >A A A >AA AA A A AAA A A
> >A A A >AA AA A A A A A AAA # get the residue indices and atom
> indices for
> >A A A selection
> >A A A >AA AA A A A A A AAA set wat_resid [$total_wat get resid];
> >A A A >AA AA A A A A A AAA set wat_index [$total_wat get index];
> >A A A >
> >A A A >AA AA A A A A A AAA # sets a counter to zero
> >A A A >AA AA A A A A A AAA set checkingNums 0;
> >A A A >
> >A A A >AA AA A A A A A AAA # for each water we do the following
> >A A A >AA AA A A A A A AAA foreach waterID $wat_index waterRID
> $wat_resid {
> >A A A >
> >A A A >AA AA A A A A A AAA # add 1 to the counter and print out
> so we know how
> >A A A far into
> >A A A >AA AA the loop we have gone
> >A A A >AA AA A A A A A AAA incr checkingNums;
> >A A A >AA AA A A A A A AAA puts " WATER $checkingNums OF
> $num_Totalwat ";A
> >A A A >AA AA A A A A A A A A
> >A A A >AA AA A A A A A AAA # make a temporary selection using the
> atom indices
> >A A A and
> >A A A >AA AA residue indicesAAA A
> >A A A >AA AA A A A A A AAA set wat_coord [atomselect top "name OW
> and index
> >A A A $waterID and
> >A A A >AA AA resid $waterRID"];A
> >A A A >AA AA A
> >A A A >AA AA A A A A A AAA # find the 4 closest waters to this
> atom
> >A A A >AA AA A A A A A AAA set close_wat [nwat 4 $wat_coord];
> >A A A >
> >A A A >AA AA A A A A AAA # Find all the possible 3 angle
> combinations which
> >A A A can be made
> >A A A >AA AA with the 4 neighbour whilst wat_coord is at the
> centre
> >A A A >AA AA AAA AAA AAA # set the sum to zero A A A A A A A A
> A A A A A A A A
> >A A A A A A A A
> >A A A >AA AA A A A A A AAA set sum 0.0;
> >A A A >
> >A A A >AA AA A A A A A AAA # set counter to zero
> >A A A >AA AA A A A A A AAA set countAngles 0;
> >A A A >AA AA A A A A A AAA for { set k 0 } { $k < 4 } { incr k }
> {A
> >A A A >AA AA A A A A A A A A A AAA for { set l [expr $k+1] } { $l
> < 4 } { incr
> >A A A l } {A
> >A A A >
> >A A A >AA AA A A A A A A A A A A A A A AAA # add 1 to angles
> which have been
> >A A A calculated.
> >A A A >AA AA The final value should always be 6 for 4 neighbours.
> >A A A >AA AA A A A A A A A A A A A A A AAA incr countAngles;
> >A A A >
> >A A A >AA AA A A A A A A A A A A A A A AAA # select the first
> atom index
> >A A A >AA AA A A A A A A A A A A A A A AAA set hi1 "[lindex
> $close_wat $k]";
> >A A A >
> >A A A >AA AA A A A A A A A A A A A A A AAA # select the reference
> atom index
> >A A A >AA AA A A A A A A A A A A A A A AAA set atom2 [[atomselect
> top "name OW
> >A A A and
> >A A A >AA AA (index $waterID and resid $waterRID)"] get {index}];
> >A A A >
> >A A A >AA AA A A A A A A A A A A A A A AAA # select the third
> atom indexA
> >A A A >AA AA A A A A A A A A A A A A A AAA set hi3 "[lindex
> $close_wat $l]";A
> >A A A >
> >A A A >AA AA A A A A A A A A A A A A A AAA # calculate the angle
> usin proc
> >A A A calcangle
> >A A A >AA AA defined earlier
> >A A A >AA AA A A A A A A A A A A A A A AAA set calcangle [angle
> $hi1 $atom2
> >A A A $hi3];
> >A A A >AA AA A
> >A A A >AA AA A A A A A A A A A A A A A AAA # add up cosine
> results to those
> >A A A from
> >A A A >AA AA previous angles
> >A A A >AA AA A A A A A A A A A A A A A AAA set sum [expr
> >A A A >AA AA {(($calcangle+1.0/3.0)*($calcangle+1.0/3.0))+$sum}];
> >A A A >AA AA A
> >A A A >AA AA A A A A A A A A A AAA }A
> >A A A >AA AA A A A A A AAA }
> >A A A >AA AA A
> >A A A >AA AA A A A A AAA # define the q value
> >A A A >AA AA A A A A A AAA set q [expr {1.0-((3.0/8.0)*$sum)}];
> >A A A >
> >A A A >AA AA A A A A A AAA # print theA A A AAA FrameA A A A A A
> AAA Loop
> >A A A >AA AA numberA A A A A A AAA Atom index A A A A A A AAA Q
> valueA A A AA
> >A A A details to
> >A A A >AA AA file
> >A A A >AA AA A A A A A AAA puts [format "Frame: $i\t ID:
> $checkingNums \t
> >A A A AtomIDX:
> >A A A >AA AA $waterID \t Q: %.2f"A AAA $q];
> >A A A >AA AA A
> >A A A >AA AA ##AAA puts $q_file "ID: $checkingNumsA A AAA
> AtomIDX: $waterIDA
> >A A A AAA Q: $q";A
> >A A A >AA AA # print details to file
> >A A A >AA AA A
> >A A A >AA AA A A A AAA }
> >A A A >
> >A A A >AA AA A A A A AAA # delete variables we do not need any
> moreAAA AAA A A
> >A A A A A
> >A A A >AA AA A A A A AAA $wat_coord delete; array unset
> wat_coord;
> >A A A >AA AA A A A A AAA $close_wat delete; array unset
> close_wat;
> >A A A >AA AA A A A A AAA $k delete; array unset k;
> >A A A >AA AA A A A A AAA $l delete; array unset l;
> >A A A >AA AA A A A A AAA $hi1 delete; array unset hi1;
> >A A A >AA AA A A A A AAA $atom2 delete; array unset atom2;
> >A A A >AA AA A A A A AAA $hi3 delete; array unset hi3;
> >A A A >AA AA A A A A AAA $calcangle delete; array unset
> calcangle;
> >A A A >AA AA A A A A AAA $sum delete; array unset sum;
> >A A A >AA AA A A A A AAA $q delete; array unset q;
> >A A A >AA AA A A A A A A
> >A A A >
> >A A A >AA AA ##A AAA puts $q_file "END"; # add END line after
> every frame
> >A A A >
> >A A A >AA AA # delete variables we do not need any moreAAA A
> >A A A >AA AA A $total_wat delete; array unset total_wat;
> >A A A >AA AA A $num_Totalwat delete; array unset num_Totalwat;
> >A A A >AA AA A $wat_resid $delete; array unset wat_resid;
> >A A A >AA AA A $wat_index delete; array unset wat_index;
> >A A A >AA AA A $waterID delete; array unset waterID;
> >A A A >AA AA A $waterRID delete; array unset waterRID;
> >A A A >AA AA A $checkingNums delete; array unset checkingNums;
> >A A A >AA AA A $countAngles delete; array unset countAngles;
> >A A A >AA AA A
> >A A A >AA AA A AAA }A
> >A A A >
> >A A A >AA AA ##A AAA puts $q_file " "; # add empty newline to
> file
> >A A A >
> >A A A >AA AA ##A AAA close $q_file;A A AAA # close file
> >A A A >
> >A A A >AA AA # delete variables we do not need any moreAAA A
> >A A A >AA AA $total_wat delete;
> >A A A >AA AA A
> >A A A >AA AA A
> >A A A >
> >A A A >AA AA END
> >A A A --
> >A A A NIH Center for Macromolecular Modeling and Bioinformatics
> >A A A Beckman Institute for Advanced Science and Technology
> >A A A University of Illinois, 405 N. Mathews Ave, Urbana, IL 61801
> >A A A [3][5]http://www.ks.uiuc.edu/~johns/AA AA AA AA AA A
> Phone: [4]217-244-3349
> >A A A [5][6]http://www.ks.uiuc.edu/Research/vmd/
> >
> > References
> >
> >A A Visible links
> >A A 1. mailto:[7]johns_at_ks.uiuc.edu
> >A A 2. mailto:[8]vmd_at_ks.uiuc.edu
> >A A 3. [9]http://www.ks.uiuc.edu/~johns/
> >A A 4. tel:[10]217-244-3349
> >A A 5. [11]http://www.ks.uiuc.edu/Research/vmd/
> --
> NIH Center for Macromolecular Modeling and Bioinformatics
> Beckman Institute for Advanced Science and Technology
> University of Illinois, 405 N. Mathews Ave, Urbana, IL 61801
> [12]http://www.ks.uiuc.edu/~johns/A A A A A A Phone:
> [13]217-244-3349
> [14]http://www.ks.uiuc.edu/Research/vmd/
>
> References
>
> Visible links
> 1. mailto:johns_at_ks.uiuc.edu
> 2. mailto:vmd_at_ks.uiuc.edu
> 3. mailto:johns_at_ks.uiuc.edu
> 4. mailto:vmd_at_ks.uiuc.edu
> 5. http://www.ks.uiuc.edu/~johns/A
> 6. http://www.ks.uiuc.edu/Research/vmd/
> 7. mailto:johns_at_ks.uiuc.edu
> 8. mailto:vmd_at_ks.uiuc.edu
> 9. http://www.ks.uiuc.edu/~johns/
> 10. tel:217-244-3349
> 11. http://www.ks.uiuc.edu/Research/vmd/
> 12. http://www.ks.uiuc.edu/~johns/
> 13. tel:217-244-3349
> 14. http://www.ks.uiuc.edu/Research/vmd/ 

-- 
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/