Variance of density/number using chunks

LAMMPS LAMMPS General Discussion
1

I would like to compute the variance in the number of atoms, binned a 3d grid using chunk/atom.

I know that I can easily compute the density in every bin using the chunks, e.g.

compute cc1 all chunk/atom  bin/3d x 0 1 y 0 1 z 0 1 units box
fix 1 all ave/chunk 100 1 100 cc1 density/number ave running file output.txt overwrite

This will return a time average of the 3d density field.

However, what is the best way to calculate the time average of the variance of the density (or number of particles per bin) as a 3d field using the chunks (e.g. cc1 above)?

2

I cannot think of a “best” way. You have the following options:

  • use file ave/chunk with 1 1 1 settings and post-process the resulting output. But that file can become large for a long simulation. Since data from adjacent time steps are highly correlated, it may not be very meaningful to output every step anyway.
  • create a new fix based on fix ave/chunk that computes the variance instead of the average
  • modify the existing fix ave/chunk to honor a new keyword that will compute the variance in addition to the average.

In general, the value of the variance is going to be tainted by the limited “statistical efficiency” in MD simulations since the collected samples are correlated and not statistically independent, while the formula to compute the variance assumes statistically independent samples. A simple test for this would be running the same simulation and analyzing variance from the output for different values of nfreq. It should reach a plateau for a sufficiently large number (which may be 1000s of steps or more).

1 Like
3

I don’t know if this code is cursed, but:

compute cz all chunk/atom bin/1d z lower 0.2 units reduced
variable counter atom 1
compute nperchunk      all reduce/chunk cz sum v_counter # or ave?
compute nperchunk_atom all chunk/spread/atom cz c_nperchunk

fix ave_nperchunk   all ave/chunk 100 10 1000 cz c_nperchunk_atom # calculates average number of atoms in chunk
fix ave_nsqperchunk all ave/chunk 100 10 1000 cz c_nperchunk_atom norm none # calculates average (number of atoms)^2 in chunk

Then \langle n^2\rangle - \langle n\rangle^2 gives you the variance.

4

Better yet:

compute cz      all chunk/atom bin/1d z lower 0.2 units reduced
variable counter atom 1
compute n1      all  reduce/chunk cz sum v_counter
compute n1a     all  chunk/spread/atom cz c_n1
compute n2      all  reduce/chunk cz sum c_n1a
compute n2a     all  chunk/spread/atom cz c_n2
# I can do this all day ...

fix ave all ave/chunk 100 10 1000 cz c_n1a c_n2a