VMD-L Mailing List

From: MannyEful E (mannyeful_at_gmail.com)
Date: Wed Mar 30 2016 - 19:51:48 CDT

Hello Everyone,

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:
 invalid character "." in expression "...st +])/[llength $list]."

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).
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. (The
errors crop up here when I try to output my results in double precision,
whilst using a proc function.)

#################### The Residence Script: #####################

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].}
}

 # process used to calculate the stdev of values in an array
proc stddev list {
 set m [lavg $list]
 expr {sqrt([mean2 $list]-$m*$m)}
 }

#################################################################

#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; # lower limit for the hydration shell1 in angstroms
set max1 10.5; # 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}; # VMD residue name of the reference atom

set sel [atomselect top all]; # select all the atoms
set all [$sel num]; # 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; # counter: If it is within region
set hist2($i) 0; # note: If it left
set hist3($i) 0; # counter: Sum of those that left
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 ";

# 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];
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]];

     $sel1 update; # updates dynamic water selection1 for this frame
     $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]); # 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 } {

incr hist3([expr $present - 1]); # add 1 if previous hist2 was 1 because
it left and came back.
puts $fp1 "HIST3: $hist3([expr $present - 1])";

}

# 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";
} 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 ";

# 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];

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 } {

incr hist3([expr $present - 1]); # add 1 if previous hist2 was 1 because
it left and came back.
puts $fp1 "HIST3: $hist3([expr $present - 1])";

}
# 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 frame
set hist2([expr $present - 1]) 1;
incr hist4([expr $present - 1]);
puts $fp1 "It will stay in the next frame";

}

}

$sel1 delete; array unset sel1; # unset the selection1
$sel2 delete; array unset sel2; # unset the selection2

}

#################################################################

# set new constants to zero
set hist_final {};
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)
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 "; <--------------
puts $fp2 "[expr {[stddev (double($hist_final))] * $stepsize}] ps : Stdev
"; <-----------

close $fp1;
close $fp2;

#################################################################

# unset all the selections for faster execution
array unset nth;
array unset unit;
array unset pol;

array unset min1;
array unset max1;
array unset fp1;
array unset fp2;
array unset atom_ref;
array unset fn;
array unset stepsize;

$sel delete; array unset sel;

array unset list2;
array unset present;

for {set i 0} {$i < $all} {incr i} {
array unset hist1($i);
array unset hist2($i);
array unset hist3($i);
array unset hist4($i);
}

array unset all;

array unset hist_final;
array unset counter;
array unset comebacks;

#################################################################

#}

END