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