Re: [ESPResSo] Analyzing the Results

[Lorenzo Isella]

Hello,
Thanks for your suggestion, but I am not quite there yet.
Basically, once I know how to read and take statistics on a
configuration, I would like to do that in a loop.
As trivial as it can sound, I have not been able to code that yet.
Any suggestions? It must be a one-liner, but it is taking me more than
I thought to get the script working.
Kind Regards

Lorenzo

On 08/10/2007, Christoph Junghans <address@hidden> wrote:
> 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
>
>