variable atom command computed for a group instead of all?

Is it possible to calculate a “variable” over a “group” of atoms, rather than all atoms?

I’m trying to have LAMMPS compute the normalized velocity autocorrelation function for a subset of atoms in my model (atoms with type=4):

group diff type == 4

fix v0 diff store/state 0 vx vy vz

variable num atom “f_v0[1]*vx + f_v0[2]*vy + f_v0[3]*vz”

variable den atom “f_v0[1]*f_v0[1] + f_v0[2]*f_v0[2] + f_v0[3]*f_v0[3]”

variable cv atom “v_num/v_den”

compute 5 diff reduce ave v_cv

The fix store/state stores the initial velocities of the atoms in group “diff” but also adds a ‘0’ in the array for all atoms not in “diff”. Because the “variable atom” command is implicit over all atoms in the array f_v0, my simulation halts with a “divide by zero” error when I get down to the line where v_cv is computed. If v_num and v_den were computed only over the group “diff” I could avoid this error. Is this possible?

Also, if anybody has a better way to compute velocity autocorrelation without post-processing, please let me know.

Thanks!

-Doug Spearot

There's no way to do this currently that I can think of.

But it's a good
idea. We could probably add an optional arg to
the variable atom command to allow specification
of a group, then
have the variable loop skip atoms not in the group
and simply put a zero in the result, like the computes
do. I'll add it to the to-do list - shouldn't be too hard.

Steve

Thanks Steve. I have a fortran script that will compute VACF for me, but it requires me to output velocities to a dump file and then perform a post-processing step. I was hoping I could save time and file space by computing VACF automatically within LAMMPS. I think what you propose will work.

-Doug

Douglas,

This is really a good way to compute VACF. I have the same pain of dumping huge data and do post-processing for getting phonon density of state. On-fly computing is really cool. Steve may think about it.

Ajing

Douglas,

This is really a good way to compute VACF. I have the same pain of dumping
huge data and do post-processing for getting phonon density of state.
On-fly computing is really cool. Steve may think about it.

i have a suggestion on the topic.

the main objection against implementing this has been the excessive
storage requirements needed for VACFs and in some sense, those
are somewhat tricky to analyse, as it can be fairly difficult to decide
what is converged or not.

there is however, a way to do this with less memory consumption
and next to no artifacts when doing the transform to frequency space
right away and then assemble convolutions. the only downside
of this approach is, that one has to define up front over how many
time steps data is to be accumulated, as that determines the
resolution of the spectral densities (one can still run longer trajectories,
but the analysis would have to be done in chunks). i have a code that
does this in post processing (as a standalone multi-threaded fortran
code and as a plugin for VMD written in multi-threaded C) and that
could be adapted to be added to LAMMPS as a fix.

however, it would make the most sense, if that would be
implemented in a way so that it would parallelize well and
can be particularly used efficiently on arbitrarily large systems
on a large parallel machine. for people running on desktops
or small clusters, i think it would currently still be better/easier/faster
to do this in postprocessing. and just tell LAMMPS to dump out
the desired data in parallel to binary files and then teach the
analysis code to read those, since disk space is so cheap.

the biggest drawback of the method that i would suggest is
that you lose the original VACFs and directly assemble the
spectral density (i.e. no phase info). at least in the code that
i've been maintaining and porting now for many years.

i'd be willing to put some effort into it (both approaches), but
can only do so with significant help of people around that
would do tests and validation of the code and provide a representative
set of usage examples that help to narrow down, how much
flexibility would be needed and how much performance would
have to be compromised for that.

cheers,
    axel.

Axel,

I see your points. There are two things. One is energy spectral. The other is VACFs. The former one contains space info but the later one does not. The downside of implementing this directly in LAMMPS with a " fix " would be one has to assemble all atoms or a group atoms in a parallel way, increasing communication time between each processor. This is true for energy spectral but may not be a big deal for VACFs since one can just integrate over time of each atoms (no needs to do communication). In terms of convergence, I think that is up to user to decide. LAMMPS can give this option for users to choose.

cheers,

Ajing

Axel,

I see your points. There are two things. One is energy spectral. The other
is VACFs. The former one contains space info but the later one does not. The
downside of implementing this directly in LAMMPS with a " fix " would be one
has to assemble all atoms or a group atoms in a parallel way, increasing
communication time between each processor. This is true for energy spectral
but may not be a big deal for VACFs since one can just integrate over time
of each atoms (no needs to do communication). In terms of convergence, I

that is not correct. you are discounting the fact that atoms move between
domains, so you would have to move the collected VACF data, i.e. the
entire velocity history with it. that would significantly increase the
communication
cost of every re-neighboring.

think that is up to user to decide. LAMMPS can give this option for users to
choose.

of course it has to be up to the user, but with VACFs the problem is that you
cannot easily tell this from the VACF itself and it also is an O(N^2) scaling
procedure. if you assemble the spectrum directly, then you exclude this up
front. it should in principle be possible to reconstruct the VACFs from the
spectrum, if one would assemble it with the phase information intact. i think
this would be worth a try, but since i didn't write the original code (i only
translated it to c and parallelized the fortran and c versions), i never looked
into the exact details of the algorithm.

as i wrote before, i am willing to look into doing the "dirty work",
but i need input from several people so that it is written in a very
generic and reusable way.

axel.

Looked into this and (re)discovered that LAMMPS is
already setup to do this correctly. Whenever an atom-style
variable is invoked, the caller has a group defined (e.g. a dump,
compute, or fix), so the variable should only be invoked on
atoms in the group, and return 0.0 otherwise.

There were a couple logic bugs in how this was working, so
it was broken. But the 23Jun11 patch fixes it, and your
commands below now appear to work.

Steve

Steve-

I downloaded the Jun23 version of LAMMPS and I'm having trouble confirming that this is working. With the following set of commands (note the thermo style output):

group diff type == 4
fix v0 diff store/state 0 vx vy vz
variable num atom "f_v0[1]*vx + f_v0[2]*vy + f_v0[3]*vz"
variable den atom "f_v0[1]*f_v0[1] + f_v0[2]*f_v0[2] + f_v0[3]*f_v0[3]"
variable cv atom "v_num/v_den"
compute 5 diff reduce ave v_cv
thermo_style custom step temp pe etotal press vol c_5

I still get a divide by zero error at the 0 timestep (first time compute 5 is evaluated):

ERROR on proc 12: Divide by 0 in variable formula
ERROR on proc 15: Divide by 0 in variable formula

This implies to me that when the compute and variable commands are invoked at the initial thermo output, it is still evaluating the variable cv over all atoms (some of which have zero initial velocity) rather than over the group diff (where all of the atoms have a non-zero velocity).

Thanks,

-Doug

Doug told me this was a problem on his production machine using
the wrong version. So the patch for this does indeed work.

Steve