proc tracemol {NT seltext} {

#Algorithmn to select NT of tracer molecules
#uniformly distributed throughout the cell volume
#mol load pdb $inputpdb.pdb
#calculate box dimensions (xdim,ydim,zdim)

set box [atomselect top "$seltext"]
set box_dim [measure minmax $box]
set max [lindex $box_dim 1]
set min [lindex $box_dim 0]
set boxsize [vecsub $max $min]
set xdim [lindex $boxsize 0]
set ydim [lindex $boxsize 1]
set zdim [lindex $boxsize 2]

set xmin [lindex $min 0]
set ymin [lindex $min 1]
set zmin [lindex $min 2]
 
#calculate K
set K [expr ceil(pow($NT,0.33))]

#calculate unit cell dimensions
set xu [expr $xdim/$K]
set yu [expr $ydim/$K]
set zu [expr $zdim/$K]

#initiallizing arrays
set X {}
set Y {}
set Z {}

set selection {}
set HHVectors {}
set Hfile [open HHfile.dat w]
set n [molinfo top get numframes]

for {set i 0} {$i < $NT} {incr i} {

#generate 3 random numbers between 1 and K
set nx [expr round(ceil(rand()*$K))]
set ny [expr round(ceil(rand()*$K))]
set nz [expr round(ceil(rand()*$K))]

set dobreak 0
for {set j 0} {$j < $i} {incr j} {
if {[lindex $X $j]==$nx && [lindex $Y $j]==$ny && [lindex $Z $j]==$nz} {
set NT [expr $NT+1]
set dobreak 1
break
}
}

if {$dobreak} {
continue
}

#add values to the arrays
lappend X $nx
lappend Y $ny
lappend Z $nz

#select all atoms inside the selected unit cell
set cell [atomselect top "name OH2 and x>($xmin+($nx-1)*$xu) and x<($xmin+$nx*$xu) and y>($ymin+($ny-1)*$yu) and y<($ymin+$ny*$yu) and z>($zmin+($nz-1)*$zu) and z<($zmin+$nz*$zu)"]
set atomlist [$cell get index]
set atomnum [llength $atomlist]
if {$atomnum==0} {
set NT [expr $NT+1]
puts "no atoms in cell"
continue
}
set atomid [lindex $atomlist 0]
lappend selection $atomid
lappend selection [expr $atomid+1]
lappend selection [expr $atomid+2] 
$cell delete

for {set k 0} {$k<$n} {incr k} {
set H1 [atomselect top "index [expr $atomid+1]" frame $k]
set H2 [atomselect top "index [expr $atomid+2]" frame $k]
set HH1 [lindex [$H1 get {x y z}] 0]
set HH2 [lindex [$H2 get {x y z}] 0]
set HH [vecsub $HH1 $HH2]
lappend HHVectos $HH
open $Hfile a
puts $Hfile "$k $atomid $HH"

$H1 delete
$H2 delete
}
}
#take no of molecules selected and print its value
set total [llength $selection]
puts "$total atoms selected"
puts "selected atoms are [list $selection]"

set atomidlist [lrange $selection 0 end]
set select [atomselect top "index $atomidlist"]
$select writepdb myout.pdb 

animate write pdb traced.pdb beg 0 end [expr $n-1] waitfor all sel $select
$select delete
close $Hfile
$box delete
}