Adding custom per-atom property to ghost atom communication

Dear LAMMPS users and developers,

I am currently writing a 'compute_pressure/atom' script, which calculates the average local pressure of every atom in a group.
To do this, one needs the average local virial. Therefore I have built a 'compute_virial/atom' script (derived from the compute stress/atom script).

In my pressure script, I create a full neighbour list and gather the average virial per atom by including the virial contribution of every atom in range of the selected atom.
See the following excerpt:

   // compute local sum of virials for each atom in group
   // add virial of atom j where distance is smaller than rcut

   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     if (mask[i] & groupbit) {
       jlist = firstneigh[i];
       jnum = numneigh[i];

       // start of loop

       for (jj = 0; jj < jnum; jj++) {
         j = jlist[jj];
         j &= NEIGHMASK;

         // dsq(i,j) = squared distance of atom i to atom j, rsq = fixed distance

         if (dsq(i,j) < rsq) {
           count[i]++;
           pressure[i][6] += virial[j][0]+virial[j][1]+virial[j][2];
         }
       }
     }
   }

   // should include this:
   // if (force->newton_pair) comm->reverse_comm_compute(this);

   for (i = 0; i < nall; i++)
     pressure[i][6] /= count[i];

where pressure[i][6] is the atomic pressure scalar.

This approach yields a rather unsatisfying result:
Ghost atoms do not hold information on the atomic virial and therefore do not contribute to the sum.
Therefore, atomic pressures near cell borders are generally lower compared to the rest.

To fix this, I would need to call 'comm->reverse_comm_compute(this)'.
As stated, ghost atoms do not hold virial information and therefore this doesn't alter the result.

How can I implement the communication of the custom 'virial/atom' tensor?
I would gladly share the scripts, once they are finished.

Thank you in advance.
C. Mink

Dear LAMMPS users and developers,

I am currently writing a 'compute_pressure/atom' script,
which calculates the average local pressure of every atom
in a group.
To do this, one needs the average local virial. Therefore
I have built a 'compute_virial/atom' script (derived from
the compute stress/atom script).

In my pressure script, I create a full neighbour list and
gather the average virial per atom by including the virial
contribution of every atom in range of the selected atom.
See the following excerpt:

​[...]

This approach yields a rather unsatisfying result:
Ghost atoms do not hold information on the atomic virial
and therefore do not contribute to the sum.
Therefore, atomic pressures near cell borders are
generally lower compared to the rest.

To fix this, I would need to call
'comm->reverse_comm_compute(this)'.
As stated, ghost atoms do not hold virial information and
therefore this doesn't alter the result.

How can I implement the communication of the custom
'virial/atom' tensor?

​you are confusing the two communication options. a reverse communication
collects data from ghost atoms (that is stored on them because of having
half neighbor lists with newton's 3rd law enabled). what you seem to be
looking for is a _forward_ communication, which distributes information
from local atoms to the corresponding ghost atoms. in order to make that
work, you need to add the corresponding pack/unpack methods to your compute
style implementation. a quick check with grep reveals, that you may look at
compute cluster/atom and voronoi/atom for examples on how to use that.

axel.​

Not clear on how what you are doing is different than

what compute stress/atom already does.

You have this line:
pressure[i][6] += virial[j][0]+virial[j][1]+virial[j][2];

What is virial[i][j] and where does it come from? Does it include

ghost atom values?

In compute stress/atom, it accesses vatom[i][j] which is

the per-atom virial from pair, bond, angle, etc interactions.

Those values may or may not have ghost atom contributions

and the code is careful to include/exlcude them, and then

do reverse communication to sum up the ghost contributions.

Unless you do something similar in your code, you likely

will not get the right answer.

Steve

(I hope this will be displayed in the thread I started yesterday).

@Axel Kohlmeyer:

Yes, I was mistaking forward and backward communication.
Thanks to your explanation I was able to point out the problem and fix it.
The local pressure detector is now working as it should be.

@Steve Plimpton:

Out of confusion I did not state my problem in an appropriate clarity.

virial[i][j] references the 'compute virial/atom' script I mentioned earlier.
I copy&pasted the virial calculation from 'stress/atom' and adapted it to a 'per atom' calculation for two reasons:
i) having the atomic virial precomputed
ii) general access for future purposes

My first draft of the 'pressure/atom' script was heavily based on 'stress/atom', but I have made major changes:
-a pressure scalar based on the average local virial and not the atomic virial is added
-the kinetic energy was replaced by the local average temperature (based on another script of mine)
-as suggested in the manual, the result is divided by the volume gathered via 'voronoi/atom', thus yielding a true pressure unit

Thank you, too!