[ESPResSo] Analyzing the Results

[Lorenzo Isella]

Dear All,
Sorry for the many questions I am posting, but I am trying to learn the 
ropes both of TCL and Espresso by running tutorial and example scripts.
I now refer to the example simulation one can find on page 13 of the 
user guide downloadable from the web page:
http://espressowiki.mpip-mainz.mpg.de/wiki/index.php/Documentation
I am coding (actually above all cutting and pasting) the example to 
understand how it works and further modify it.
However, I am having some problems.
The first part of the example works fine (I show the code in the 
following, but I simply added the visualization of particle 
configuration as long as the integration progresses):


set n_part 200; set density 0.7
set box_l [expr pow($n_part/$density,1./3.)]

setmd box_l $box_l $box_l $box_l
setmd periodic 1 1 1

set q 1; set type 0
for {set i 0} { $i < $n_part } {incr i} {
set posx [expr $box_l*[t_random]]
set posy [expr $box_l*[t_random]]
set posz [expr $box_l*[t_random]]
set q [expr -$q]; set type [expr 1-$type]
part $i pos $posx $posy $posz q $q type $type
}

setmd time_step 0.01; setmd skin 0.4
set temp 1; set gamma 1
thermostat langevin $temp $gamma

set sig 1.0; set cut [expr 1.12246*$sig]
set eps 1.0; set shift [expr 0.25*$eps]
inter 0 0 lennard-jones $eps $sig $cut $shift 0
inter 1 0 lennard-jones $eps $sig $cut $shift 0
inter 1 1 lennard-jones $eps $sig $cut $shift 0
inter coulomb 10.0 p3m tunev2 accuracy 1e-3 mesh 32


set integ_steps 200

for {set cap 20} {$cap < 200} {incr cap 20} {
puts "t=[setmd time] E=[analyze energy total]"
inter ljforcecap $cap; integrate $integ_steps
}
inter ljforcecap 0

puts "end of equilibration"


for {set i 0} { $i < 20 } { incr i} {
set temp [expr [analyze energy kinetic]/(1.5*$n_part)]
puts "t=[setmd time] E=[analyze energy total], T=$temp"
integrate $integ_steps

puts "save configuration"

set f [open "config_$i" "w"]
blockfile $f write tclvariable {box_l density}
blockfile $f write variable box_l
blockfile $f write particles {id pos type}
close $f
}

puts "end of integration"



And everything is fine up to now. Then I want to calculate the radial 
density function (rdf). Say I want to have it for configuration 12. Then 
the following snippet works fine:

set f [open "config_12" "r"]
while { [blockfile $f read auto] != "eof" } {}
close $f

puts "ok reading the block file"

set rdf [analyze rdf 0 0 0.9 [expr $box_l/2] 100]
puts "the box size is $box_l"
puts "the radial density is $rdf"
set rlist ""
set rdflist ""
foreach value [lindex $rdf 1] {
lappend rlist [lindex $value 0]
lappend rdflist [lindex $value 1]
}

The documentations shows the following code to add to a second script to 
calculate the rdf for each configuration and take an average:

set cnt 0
for {set i 0} {$i < 100} {incr i} { lappend avg_rdf 0}
foreach filename $argv {
set f [open $filename "r"]
while { [blockfile $f read auto] != "eof" } {}
close $f
set rdf [analyze rdf 0 1 0.9 [expr $box_l/2] 100]
set rlist ""
set rdflist ""
foreach value [lindex $rdf 1] {
lappend rlist [lindex $value 0]
lappend rdflist [lindex $value 1] }
set avg_rdf [vecadd $avg_rdf $rdflist]
incr cnt
}
set avg_rdf [vecscale [expr 1.0/$cnt] $avg_rdf]

But this does not work by itself or added to the previous script...
I suppose it has to do with the set f [open $filename "r"], but I have 
not been able to figure out what to change yet (still learning TCL).
But am I missing the obvious? The code should work without any twisting.
Many thanks

Lorenzo