Adding a new per atom variable in the source code

    Dear Axel,

    I'm trying out both the fix:array and the find_custom options you
    proposed.

    I'm now able to define new variables through the lammps_command
    function and find them again through the code snippet below.

    int tmp = -1;
    int n = lmp->atom->find_custom("d_xDuDt",tmp);
    double *dvector = lmp->atom->dvector[n];

    How can I now write to individual atoms?

    Given that I have nlocal atoms on my MPI rank and I know the index
    of the atom of interest "atomI", I'd like to write a new double to
    dvector[atomI] or to lmp->atom->dvector[n][atomI]. If tried it in
    this manner but I get an error and gdb isn't very informative about
    it (only gives me the line).

Long story short, problem in chair, not in computer. Thanks for the
elaborate reply though. I had to search for "xDuDt" rather than
"d_xDuDt", otherwise you get -1 as a dvector index and the code crashed
due to out of bound array indexing. Could be prevented by a simple "if(n
< 0)" check.

​it is not clear what you are referring to as "index". atoms in LAMMPS
have an "atom id", which is stored in atom->tag. this "tag" value is
the atom id used in a data file or the "id" field in a dump file. these
allow to address atoms globally.

With index I mean the index of the atom in a for loop over local atoms
on the rank/ local domain.

I keep track of them by using a lmpCpuId variable which is just a place
holder for the MPI_Comm_rank() output. This way I only send stuff to the
local domain atoms. The reason this is needed is that OpenFoam has a
mesh based decomposition and a particle may reside on a different processor.

In the coupling lib function I match the "atom id" between the incoming
info and the local atoms by sorting the tag arrays and looping over the
indexes (I hope the way I put it here make sense to you). That way I
don't have to care about the order are in, I jsut make sure the atom
od's are unique (they should be of course).

This was all pre-implemented in the sediFoam framework by Heng Xiao et
al. ( GitHub - xiaoh/sediFoam: CFD-DEM Solver with emphasis on sediment transport ). I just ported the whole thing
to the latest lammps version and modified some things on the openfoam side.

however, the per-atom storage, is associated with atoms and they are
distributed and in a different order. you have the "local" atoms​ (i.e.
atoms owned by a subdomain == MPI rank) between index 0 and nlocal and
"ghost" atoms (copies of atoms from neighboring subdomains up to the
communication cutoff) between nlocal and nlocal+nghost. to get the
global atom id from a local or ghost atom, you look at atom->tag[i] you
will get the global atom id of an atom, and the function atom->map(id)
will return the local atom index or -1 from a global atom id. index ==
-1 means, this atom is not present on the local MPI rank.

Question though, will the following command give me back ghost atom tags
or only local atom tags ?

LAMMPS *lmp = (LAMMPS * ) ptr;
int *tag = lmp->atom->tag

if you are looping over atom indices from neighbor lists, be aware that
those are using local indices and in the jlist, atoms indices may be
"decorated" with extra bits, when an atom is connected via a bond or
angle or dihedral interaction and "special_bond" scaling factors are
applied.

I just do granular stuff without bonds etc, but I'll look into this as a
precaution.

if you want gdb to be more helpful, please compile without optimization
and with -g. then gcc won't optimize away loop indices or other
temporary variables and you can look up their values as expected.

I still had some optimizing flags in the makefile, lazy of me sorry.
Removed them and now it's more helpful.

​[...]​

Question though, will the following command give me back ghost atom tags
or only local atom tags ?

LAMMPS *lmp = (LAMMPS * ) ptr;
int *tag = lmp->atom->tag

​this will work for local and ghost atoms​

> if you are looping over atom indices from neighbor lists, be aware that
> those are using local indices and in the jlist, atoms indices may be
> "decorated" with extra bits, when an atom is connected via a bond or
> angle or dihedral interaction and "special_bond" scaling factors are
> applied.
>
I just do granular stuff without bonds etc, but I'll look into this as a
precaution.

​ignoring the special_bonds bits is done by simply applying & NEIGHMASK to
the innerloop indices (usually indexed by j)​.

​axel.​