variance of the position of each atom during a run

Dear all
I want to calculate the variance of the position of each atom during a
run. First I calculate the mean position and store it in a variable
called "mean_pos"

...
fix mean_posf all ave/atom 1 1000 1000 x y
run 1000
variable mean_pos0 atom f_mean_posf
variable mean_pos atom ${mean_pos0} # just to store this per
atom vector for later use
run 0

then I clear everything and restart from initial configuration.

...
atom_modify map array
....
variable sigma2 atom (x-mean_pos[1])^2+(y-mean_pos[2])^2
fix VARf all ave/atom 1 1000 1000 v_sigma2
run 1000
dump VARd all custom 1 test.dump id type x y f_VARf
dump_modify VARd format "%d %d %20.15g %20.15g %20.15g" flush yes
first yes sort id
run 0

But I get the the error "ERROR on proc 0: Substitution for illegal
variable" while lammps try to execute "variable mean_pos atom
\{mean\_pos0\}" and if I try to remove in ${mean_pos0} to have "variable
mean_pos atom mean_pos0" I get an error while while lammps tries to
substitute a value for variable sigma2.
(ERROR on proc 0: Invalid atom vector in variable formula)

Do you know where I'm doing a mistake?

Best,
Nima

variable mean_pos atom ${mean_pos0} # just to store this per atom vector for later use

This will not work. Atom and equal style varaibles don't store
anything. They just store formulas. Nor will per-atom quantities
in the fix mean_posf be retained after you do a "clear". That will
destroy the fix and what it stores.

What is it you want to do with old positions?

Steve

I want to calculate the variance of the particle positions during a short run. At the end I’ll have a per-atom quantity which is the variance of particle positions.

First I do a run in order to find the mean position of each particle. I should store this per-atom quantity in a variable “mean_pos” to be later used. Then I do the same run again, this time I have mean_pos stored, so I can calculate the variance with fix ave/atom. in other words in 2D I want the following time averaged per-atom quantity:

< (x-mean_pos[1])^2+(y-mean_pos[2])^2 >

Can this algorithm be implemented in lammps? Or you might have a better suggestion…

Best,
Nima

I want to calculate the variance of the particle positions during a short
run. At the end I'll have a per-atom quantity which is the variance of
particle positions.

First I do a run in order to find the mean position of each particle. I
should store this per-atom quantity in a variable "mean_pos" to be later
used. Then I do the same run again, this time I have mean_pos stored, so I
can calculate the variance with fix ave/atom. in other words in 2D I want
the following time averaged per-atom quantity:

< (x-mean_pos[1])^2+(y-mean_pos[2])^2 >

Can this algorithm be implemented in lammps? Or you might have a better
suggestion...

just do it in post-processing after running a regular trajectory.
VMD for example has a "measure rmsf" command that can be
used for that purpose.

axel.

Thanks Axel,
I already thought to do post processing. However it is very much preferable for my application to do everything inside the input script. I have a large system and there will be even a second averaging over an ensemble of ~1000 samples involved. It saves a lot of disk space and effort if it can be simply implemented in Lammps.

If it would be possible to save a per-atom quantity in a variable then this task could have been done in Lammps, right?

the other possibility is to have a loop over particles and just use “equal” style instead of “atom” style in variable command. Would you do that like this?

Best,
Nima

Again, variables don't store values, they store formulas. The
fix store/state command will store a per-atom vector(s) of
values, like the position at some point in time. Then you
can write a variable with a per-atom formula that computes
the variance you list below, using the values stored by
the fix as the reference state. Depending on when/how you
evaluate that formula it will then calculate the variance
at later time relative to the stored values. You could use
that as input to a fix ave/atom to average and output it.

Steve

Thanks, fix store/sate is the solution. Now I can store the mean position of particles.

In the second step I want to restore the initial condition in order to perform an identical run again (having the mean positions, the variance will be calculated during this second run). Unfortunatly as Steve mentioned, this cannot be done via ( clear & read_restart ) commands, since the stored state in the fix store/state will be lost.

Is there a way to reset all velocities and positions to initial values, without clearing everything?

Best,
Nima

Thanks, fix store/sate is the solution. Now I can store the mean position of
particles.

In the second step I want to restore the initial condition in order to
perform an identical run again (having the mean positions, the variance will
be calculated during this second run). Unfortunatly as Steve mentioned, this
cannot be done via ( clear & read_restart ) commands, since the stored state
in the fix store/state will be lost.

Is there a way to reset all velocities and positions to initial values,
without clearing everything?

why do you need to reset?

your calculation only makes sense, if your system is in equilibrium.
thus for converged results you first have to run long enough to get
your average positions (i assume there is no diffusion) and then
you can just continue the trajectory. after all, even if you reset, to
the initial conditions, you will not get an identical trajectory!
trajectories are only bitwise reproducible if you use fixed point
math and nve integration.

cheers,
    axel.

Actually I want to calculate the variance on short time scales (much smaller than the alpha relaxation time of the system), so continuing without restarting from initial positions will probably add a bias to the calculated value.

I think the ability to restart from given configuration without clearing everything would be in general useful. for example if you want to calculate the iso-configuration ensemble average of a quantity (like propensity in the context of supercooled liquids which was introduced by P. Harrowell) you would need to restart always from same configuration but with a randomly chosen set of velocities from MB distribution.

If there is no way to do this at the moment, it may be possible: with a generalization of “set” command to make it capable of setting all velocities and/or positions to corresponding per-atom variables. What do you think?

Best,
Nima

You can always write a restart file at any point during your
input script (via the write_restart command), which effectively
allows you to restart from that state point. Either in another
acript, or via clear, read_restart. If you are simply wanting
to archive the coords stored in fix store/state, then those
are stored to the restart file, so if you restarted, you would
also retrieve the earlier stored state.

Steve

The problem is that, I want to restart from the initial point. The point where fix store/state (or fix ave/atom) has not yet been called.

Nima

Can you save a restart file at that point, then read it back in?

Steve

I can always write restart files and read them back but,
If I write the restart file at the beginning where I want my system to be restated from again. by reading back that restart, I don’t have the desired information calculated by fix ave/atom at later times.
if I write the restart at the end of the first run, then I have the information of the fix but the positions are not those I want to restart from.

The solution could be to reset just positions & velocities to initial values without killing all information stored in fixes.

Nima

I'm not clear why it would be useful to have old positions and the
new fix ave/atom averages. What could you do with both
those quantities?

To restore old positions, some new code would have to be written.
It's not straightforward b/c of issues like: what if the box size
changed? what if the count of atoms has changed? What if
atoms are a long way from where they started?

Steve

One application is to do iso-configurational averages. members of isoconfigurational ensemble have the same configuration but with different set of velocities chosen from MB distribution. This idea helped for example to understand better the dynamical heterogeneity in supercooled liquids and whether it has any structural origin. If you want to do an iso-configurational average, you’ll have to keep information stored in fix ave/atom and to be able to restore to original configuration (new member of ensemble).

regarding the difficulties you mentioned, would it make it easier to restore to original configuration via a restart file, and then feeding in the data which was stored in fixes? in other words adding the option to fix ave/atom to be able to start averaging instead of scratch from previous point, that could be read in from disk in the worst case.

The other option -which I’ve tried to avoid so far- is just do post processing.

Best,
Nima

regarding the difficulties you mentioned, would it make it easier to restore
to original configuration via a restart file, and then feeding in the data
which was stored in fixes? in other words adding the option to fix ave/atom
to be able to start averaging instead of scratch from previous point, that
could be read in from disk in the worst case.

I don't think of fix ave/atom as storing old info. Rather it accumulates
info for a short amount of time, then writes it to a file. Again and again.
So I don't see any great advantage to trying to make fix ave/atom continue
on. Can you just write its averages to a file (it does this for you).
Then when
you restart the old config, give it new velocities, you also open a new
fix ave/atom file and write the ave for that run to a new file. Rinse
and repeat.

The other option -which I've tried to avoid so far- is just do post
processing.

If the suggestions I've made don't do what you want, then post-processing
is probably your best option.

Steve