Reverse Communication for forces

LAMMPS LAMMPS Development
1

Dear LAMMPS Users,

I have developed a new fix that biases the number of water molecules in the hydration layer surrounding specific atom types or IDs. This fix applies a force to water molecules within a certain distance of the selected target atoms, while simultaneously applying an equal but opposite force to the predefined atoms themselves. Since the target atom types could include ghost atoms, I added pack/unpack_reverse_comm after the post_force calculations.

When I run the code, the simulation box is divided into rectangular regions corresponding to the number of MPI processes used (snapshot attached!). However, if I remove the reverse communication in my code, these regions do not form, but I’m certain the calculation is incorrect without it.

Am I misunderstanding how force communication should be handled? I’ve included the relevant code snippet below.

Best,
Hadi

int Fixapplyforce::pack_reverse_comm(int n, int first, double *buf)
{
int i;
int m = 0;
int last = first + n;

for (i = first; i < last; i++) {
buf[m++] = atom->f[i][0];
buf[m++] = atom->f[i][1];
buf[m++] = atom->f[i][2];
}

return m;
}

void Fixapplyforce::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j;
int m = 0;

for (i = 0; i < n; i++) {
j = list[i];
atom->f[j][0] += buf[m++];
atom->f[j][1] += buf[m++];
atom->f[j][2] += buf[m++];
}
}

test1
test14000×3000 5.27 MB

2

Actually, that piece is only showing the pack/unpack. What is more important is what you do to compute the forces.

Here are two general thoughts about this issue:

  • since you are using the original force array, you probably have to set the forces on the ghost atoms to zero before you compute and apply forces. The reverse force operation on the force array should have already happened before the post force step, but that sums the forces from the ghost, yet does not set the ghost atom forces to zero (for efficiency reasons). Thus you have to wipe them to avoid double counting.
  • when adding forces to pairs of atoms (do you use a neighbor list for that?) where one of the two atoms is a ghost atom, you have to make certain to not add those forces twice. If you use a neighbor list, it depends on the type of neighbor list and whether you have newton on for pairs or off. same as for pair styles. If you use a full neighbor list and thus process each pair twice, you need to scale the forces and energies by 0.5 and don’t need the reverse communication. For a half neighbor list with newton pair off all pairs straddling the domain boundaries are listed twice (once on each sub-domain) and thus you only need to use the forces on the local atoms (but scale the energy contribution by 0.5) and use no reverse communication. Only with a half neighbor list and newton pair on, you need the reverse communication.

For reference please see the recent LAMMPS paper and 4.4. Parallel algorithms — LAMMPS documentation and 4.6. Communication patterns — LAMMPS documentation.

3

Thank you Axel.

Your points were crucial in helping me resolve the issue, as I implemented them into my code. I created two versions of the code—one with a neighbor list and one without—but both had the same problem initially. Setting the initial force value for the ghost atom was the key to fixing it.

Thanks again,
Hadi