how to calculate the ensemble averge of a variable

Dear all:

Using lammps, it is easy to output variables at each timestep, for example, pxy, volume and so on. Then how to calculate the ensemble average for each variable using lammps? Using fortran code, it has something like: sum_i=0; sum_i=sum_i+xi for x output at different step i; then avg_i=sum_i/real(i). I tried to use fix ave/time to do that by changing Nevery, Nrepeat and Nfreq values but could not get satisfying result. Can I get some hint here? thanks a lot.

Laura

2011/10/12 zhanglaura <[email protected]...>:

Dear all:

Using lammps, it is easy to output variables at each timestep, for example,
pxy, volume and so on. Then how to calculate the ensemble average for each
variable using lammps? Using fortran code, it has something like: sum_i=0;
sum_i=sum_i+xi for x output at different step i; then avg_i=sum_i/real(i). I
tried to use fix ave/time to do that by changing Nevery, Nrepeat and Nfreq
values but could not get satisfying result. Can I get some hint here? thanks

what is a "non-satisfying" result?

the explanation of the three numbers is pretty detailed.
there is very little that can be added there without just
restating what is there. if you look at this very carefully
and perhaps write down a few examples of what would
have to happen when, you should be able to figure it out.

other than that, perhaps you can give us an example
of what you tried and what you expected and what
you got instead.

the more specifics you provide, the easier it is to
give specific advice. the general advice is simply:
re-read the documentation. it _is_ explained there,
you probably just haven't see it from the right angle yet.

cheers,
     axel.

Hi, Axel:

I did one test with important part of the codes like following:


variable pxy equal pxy
variable pxz equal pxz
variable pyz equal pyz

fix 3 all ave/spatial 100 10 1000 y 0.0 0.05 vx vz norm sample file vel.profile units reduced

fix 4 all ave/time 1000 1 1000 v_pxy file avg_pxy.file

thermo 1000

thermo_style custom step temp pe etotal press v_pxy v_pxz v_pyz vol

Then I got the data in avg_pxy.file being exactly the same as in log.lammps file which output the thermo_style contents. It did not do average based on data output every 1000 at all. That is why I said it was not satisfying.

Another test I did is following with all other parts the same as above:

fix 4 all ave/time 1 1 1 v_pxy file avg_pxy.file

The avg_pxy.file showed data fluctuating even larger than the original pxy output from log.lammps file. So I am wondering what I should do. Thanks a lot.

Laura

2011/10/12 zhanglaura <[email protected]...>:

Hi, Axel:

I did one test with important part of the codes like following:

.....
variable pxy equal pxy
variable pxz equal pxz
variable pyz equal pyz

fix 3 all ave/spatial 100 10 1000 y 0.0 0.05 vx vz norm sample file
vel.profile units reduced

fix 4 all ave/time 1000 1 1000 v_pxy file avg_pxy.file

thermo 1000

thermo_style custom step temp pe etotal press v_pxy v_pxz v_pyz vol

......

Then I got the data in avg_pxy.file being exactly the same as in log.lammps
file which output the thermo_style contents. It did not do average based on
data output every 1000 at all. That is why I said it was not satisfying.

Another test I did is following with all other parts the same as above:

fix 4 all ave/time 1 1 1 v_pxy file avg_pxy.file

The avg_pxy.file showed data fluctuating even larger than the original pxy
output from log.lammps file. So I am wondering what I should do. Thanks a
lot.

neither of those choices make sense to me.
have a look at the example that is given in the
documentation:

For example, if Nevery=2, Nrepeat=6, and Nfreq=100, then values on
timesteps 90,92,94,96,98,100 will be used to compute the final average
on timestep 100. Similarly for timesteps 190,192,194,196,198,200 on
timestep 200, etc. If Nrepeat=1 and Nfreq = 100, then no time
averaging is done; values are simply generated on timesteps
100,200,etc

if you want o average over all frames, you have
to set Nevery to 1, if you want output every 1000 frames,
you have to set Nfreq to 1000, the maximum you can
average over is Nfreq/Nevery, i.e. 1000.
so i would try using

fix 4 all ave/time 1 1000 1000 v_pxy file avg_pxy.file

since immediately succeeding values in an MD are highly
correlated, you don't need to have all frames included, so
doing:

fix 4 all ave/time 10 100 1000 v_pxy file avg_pxy.file

should give you the statistically equivalent output.
if you want to average over a longer interval, you have
to scale the numbers accordingly.

is that more helpful?
i've been really just paraphrasing
the example in the documentation,
since there is not much else to say.

axel.

You can have fix ave/time just print un-averaged
values one by one directly to a file. Then if you
average them you should be able to tell exactly
what it averaged.

Steve