Dear Loronzo,
I think you problem is the line
"foreach filename $argv {"
it take the argument of the command line (stored in $argv) as files.
I think you should as the names of the file (config_12 etc.) behind the
script.
For later use I would prefer <rdf> (internal averaging) anyway.
Bye,
Christoph
PS: Never mind to ask question, that help us to improve Espresso (and
the comments in the scripts). ;-)
Lorenzo Isella wrote:
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
_______________________________________________
ESPResSo mailing list
address@hidden
https://fias.uni-frankfurt.de/mailman/listinfo/espresso
This email was Anti Virus checked by Astaro Security Gateway.
http://www.astaro.com
--
Christoph Junghans
Max Planck Institute for Polymer Research
Theory Group
POBox 3148
D 55021 Mainz, Germany
Phone: +49 6131 379 335
Web: http://www.mpip-mainz.mpg.de/~junghans