Extending `stress/cartesian` to include bond forces: how to grow `num_bond`/`bond_atom` arrays to include ghost atoms?

LAMMPS LAMMPS Development
1

As discussed at the end of this bug report, I am trying to extend compute stress/cartesian to include bond forces.

Initial implementation

Initially, I came up with this code:

--- a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp
+++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp
@@ -361,6 +362,15 @@ void ComputeStressCartesian::compute_array()
       // Check if inside cut-off
       if (rsq >= cutsq[itype][jtype]) continue;
       pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
+      // Loop over all bonds, checking if bond is between current pair
+      for (int bondi = 0; bondi < neighbor->nbondlist; bondi++) {
+        int *bond = neighbor->bondlist[bondi];
+        if (bond[0] == i && bond[1] == j) {
+          // Add bond force to pair force
+          force->bond->single(bond[2], rsq, i, j, fpair);
+          break;
+        }
+      }
       compute_pressure(fpair, xi1, xi2, delx, dely, delz);
     }
   }

This works, but unfortunately is excruciatingly slow. (not unexpectedly, obviously, since what I’m doing here is looping over all bonds in the system for each pair of atoms in the neighbor list.

Note that this code requires you to set special_bonds to > 0, since special neighbors will be excluded from the neighbor lists when set to exactly 0. Hence, instead of special_bonds lj/coul 0.0 1.0 1.0, for this to work it is required to set special_bonds lj/coul 1.0e-100 1.0 1.0 instead.

Optimization

In order to improve the performance, I tried another approach where instead of looping through all bonds in the system, I specifically loop over the bonds of the current atom:

--- a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp
+++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp
@@ -361,6 +362,15 @@ void ComputeStressCartesian::compute_array()
       // Check if inside cut-off
       if (rsq >= cutsq[itype][jtype]) continue;
       pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
+      if (atom->molecular == Atom::MOLECULAR) compute_bond_force(i, j, rsq, fpair);
       compute_pressure(fpair, xi1, xi2, delx, dely, delz);
     }
   }
@@ -409,6 +407,34 @@ void ComputeStressCartesian::compute_array()
   }
 }
 
+void ComputeStressCartesian::compute_bond_force(int atom1, int atom2, double rsq, double& f) {
+  // Loop over bonds of atom1
+  for (int i_bond = 0; i_bond < atom->num_bond[atom1]; i_bond++) {
+    // Get local index of bonded atom
+    int bonded_atom = atom->map(atom->bond_atom[atom1][i_bond]);
+    if (atom2 == bonded_atom) {
+      // Add bond force to pair force
+      force->bond->single(atom->bond_type[atom1][i_bond], rsq, atom1, atom2, f);
+      return;
+    }
+  }

This works if newton_bond=false, since in that case the num_bond/bond_atom contains bonds ‘symmetrically’ (for both atoms) (as I found out in a different thread, so I only need to loop over bonds of the ‘left’ atom.

However, for newton_bond=true, I need to loop over bonds of both the ‘left’ and ‘right’ atom in the pair. This is more complicated since now the right atom can be a ghost atom, in which case num_bond/bond_atom does not contain it and I need to manually communicate this information.

As a first try, I came up with this:

--- a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp
+++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp
@@ -176,6 +176,8 @@ ComputeStressCartesian::ComputeStressCartesian(LAMMPS *lmp, int narg, char **arg
   memory->create(tpcyy, nbins1 * nbins2, "tpcyy");
   memory->create(tpczz, nbins1 * nbins2, "tpczz");
   memory->create(array, size_array_rows, size_array_cols, "stress:cartesian:output");
+
+  comm_forward = atom->bond_per_atom + 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -252,6 +254,9 @@ void ComputeStressCartesian::compute_array()
   // invoke half neighbor list (will copy or build if necessary)
   neighbor->build_one(list);
 
+  // We have to communicate bonds of ghost atoms if newton is on for bonds
+  if (force->newton_bond) comm->forward_comm(this);
+
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
@@ -361,16 +366,9 @@ void ComputeStressCartesian::compute_array()
 
       // Check if inside cut-off
       if (rsq >= cutsq[itype][jtype]) continue;
+
       pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
-      // Loop over all bonds, checking if bond is between current pair
-      for (int bondi = 0; bondi < neighbor->nbondlist; bondi++) {
-        int *bond = neighbor->bondlist[bondi];
-        if (bond[0] == i && bond[1] == j) {
-          // Add bond force to pair force
-          force->bond->single(bond[2], rsq, i, j, fpair);
-          break;
-        }
-      }
+      if (atom->molecular == Atom::MOLECULAR) compute_bond_force(i, j, rsq, fpair);
       compute_pressure(fpair, xi1, xi2, delx, dely, delz);
     }
   }
@@ -409,6 +407,34 @@ void ComputeStressCartesian::compute_array()
   }
 }
 
+void ComputeStressCartesian::compute_bond_force(int atom1, int atom2, double rsq, double& f) {
+  // Loop over bonds of atom1
+  for (int i_bond = 0; i_bond < atom->num_bond[atom1]; i_bond++) {
+    // Get local index of bonded atom
+    int bonded_atom = atom->map(atom->bond_atom[atom1][i_bond]);
+    if (atom2 == bonded_atom) {
+      // Add bond force to pair force
+      force->bond->single(atom->bond_type[atom1][i_bond], rsq, atom1, atom2, f);
+      return;
+    }
+  }
+
+  // Only have to loop over bonds of atom2 if newton is on (num_bond and bond_atom only contain 'owned' bonds if
+  // newton is on)
+  if (!force->newton_bond) return;
+
+  for (int i_bond = 0; i_bond < atom->num_bond[atom2]; i_bond++) {
+    // Get local index of bonded atom
+    int bonded_atom = atom->map(atom->bond_atom[atom2][i_bond]);
+    if (atom1 == bonded_atom) {
+      // Add bond force to pair force
+      force->bond->single(atom->bond_type[atom2][i_bond], rsq, atom1, atom2, f);
+      return;
+    }
+  }
+
+}
+
 void ComputeStressCartesian::compute_pressure(double fpair, double xi, double yi, double delx,
                                               double dely, double delz)
 {
@@ -488,3 +514,48 @@ double ComputeStressCartesian::memory_usage()
 {
   return (14.0 + dims + 7) * (double) (nbins1 * nbins2) * sizeof(double);
 }
+
+int ComputeStressCartesian::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) {
+  // n: number of ghost atoms
+  // list: atom ids
+
+  // Pack num_bond and bond_atom arrays for each ghost atom
+
+  int m = 0;
+
+  // Loop over ghost atoms
+  for (int i = 0; i < n; i++) {
+    int atom1 = list[i];
+    int num_bond = atom->num_bond[atom1];
+    buf[m++] = ubuf(num_bond).d;
+
+    // Loop over bonds for ghost atom
+    for (int j = 0; j < num_bond; j++) {
+      // Get global index of bonded atom
+      buf[m++] = ubuf(atom->bond_atom[atom1][j]).d;
+    }
+  }
+
+  return m;
+}
+
+void ComputeStressCartesian::unpack_forward_comm(int n, int first, double *buf) {
+  int m = 0;
+
+  int last = first + n;
+  for (int i = first; i < last; i++) {
+    int num_bond = (int) ubuf(buf[m++]).i;
+    atom->num_bond[i] = num_bond;
+    for (int j = 0; j < num_bond; i++) {
+      // Get global index of bonded atom from buffer
+      atom->bond_atom[i][j] = (int) ubuf(buf[m++]).i;
+    }
+  }
+}

However, this segfaults in unpack_forward_comm(), where it is assigning atom->num_bond[i] = num_bond; I presume this array is only nlocal large when newton_bond=true, so the ghost atoms don’t fit.

This is where I’m stuck. I presume I should grow the num_bond/bond_atom arrays, but haven’t yet found out how. That’s my question: how should I grow num_bond/bond_atom when newton is on, so that I can communicate the bond info of ghost atoms? Or this there a better way to do this?

2

Why not do what LAMMPS does, and run through bonds after running through pairs, instead of during? That way you don’t need to tweak special_bonds, and you can take whatever logic the pair calculation uses to deal with newton flags and reuse it with different newton_bonds settings. The compute_pressure function (which really should be modernized with tallies) seems to take all the necessary quantities without worrying about whether fpair has been calculated from a pair style or bond style.

3

Well, frankly, I’m not sure :slight_smile: It might be that there is a simpler way to accomplish this.

Because the computation inherently deals with pairs of atoms, I thought the easiest way to add bond forces is to check if a pair of atoms is bonded, and if it is, to simply add the bond force to the pair force fpair. Also, the code in compute_array() that deals with binning and ghost atoms / newton is a bit arcane to me, and this way I don’t have to modify/re-implement it.

4

That is a very, VERY bad idea. If you change that information, you change the system topology and thus the force field. More importantly, these kind of properties are locked in when you define the box and cannot be changed. This is why commands like create_box or read_data have all those extra flags to reserve extra space in case it would be required later.

5

I see. What would be the proper course of action to communicate bond information of ghost atoms in the case that newton is on? Just store it in a separate array private to the ComputeStressCartesian class?

6

Hmm, I am now trying to store the bond information in separate arrays.

The pack_forward_comm() and unpack_forward_comm() methods are still relatively the same:

@@ -488,3 +535,52 @@ double ComputeStressCartesian::memory_usage()
 {
   return (14.0 + dims + 7) * (double) (nbins1 * nbins2) * sizeof(double);
 }
+
+int ComputeStressCartesian::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) {
+  int m = 0;
+
+  for (int i = 0; i < n; i++) {
+    int atom1 = list[i];
+
+    int num_bond = atom->num_bond[atom1];
+    buf[m++] = ubuf(num_bond).d;
+
+    // Loop over bonds for ghost atom
+    for (int j = 0; j < num_bond; j++) {
+      // Get global index of bonded atom
+      buf[m++] = ubuf(atom->bond_atom[atom1][j]).d;
+      // Get bond type
+      buf[m++] = ubuf(atom->bond_type[atom1][j]).d;
+    }
+  }
+
+  return m;
+}
+
+void ComputeStressCartesian::unpack_forward_comm(int n, int first, double *buf) {
+  int m = 0;
+
+  for (int i = 0; i < n; i++) {
+    int num_bond = (int) ubuf(buf[m++]).i;
+    ghost_num_bond[i] = num_bond;
+    for (int j = 0; j < num_bond; j++) {
+      // Get global index of bonded atom from buffer
+      ghost_bond_atom[i][j] = (int) ubuf(buf[m++]).i;
+      // Get bond type from buffer
+      ghost_bond_type[i][j] = (int) ubuf(buf[m++]).i;
+    }
+  }
+}

However, I now encounter a different issue: in pack_forward_comm(), I access atom->num_bond[] out of bounds. In other words, the atom1 get get from list, does not exist in atom->num_bond. I don’t really get why. As far as I understood, pack_forward_comm() is called with a number of ‘owned’ atoms (n) needing to be communicated to other processors (where they are ghost atoms), which you get from the list array. I have indeed checked that at some point atom1 > atom->nlocal. I don’t get how this can be, if they are supposed to be ‘owned’ atoms. Am I missing something?

7

For one thing, your unpack is not correct, you need to use the first variable, not 0. Since LAMMPS does multiple passes of pack and unpack, it could show up in the next pack.

Here is a good simple example:

8

The code currently doesn’t use comm->forward. Either it’s already wrong (in which case you need to fix that first), or it’s already right – in which case you should just re-use the pair logic for bond forces, which I imagine aren’t any different from pair forces in the stress calculation. (And if you’re worrying about bond forces, what about angle and dihedral forces?)

It looks like the tpcDD arrays (D = x, y, or z) cover the entire simulation box, which is how the MPI_Allreduce calls in compute_array work (to-do: check if MPI_IN_PLACE maintains performance, because that could halve the memory requirements, and also see if the tXXXX arrays can be collapsed into one to save on the overhead of multiple MPI_Allreduce calls).

This makes sense because the number of bins should never be large enough to make the memory requirements prohibitive, and the logic of “localising bins” to a proc’s spatial subdomain would take some work (especially without guarantees that the bin boundaries will align neatly with a proc’s domain boundaries) (although it’s not too difficult to adapt from other parts of the LAMMPS code, notably the PPPM, but it was certainly judged work not worth doing by the original devs) (since stress is accumulated into bins across the line between particles, for a large enough cutoff, procs will always have to deposit entries into each others’ boxes, which would be prohibitive without doing things this way).

So the code’s logic is that if a proc holds a local<->ghost pair, it can simply deposit the accumulated (F.displacement) into the correct bin, because all procs are tallying into all bins across the simulation box. You can see this at work in the tpcDD[ADDRESS] += (fpair * delD * delD * (lb - la)); lines in compute_pressure. Based on that logic, you can loop over bonds and call the same compute_pressure function on whatever the equivalent of pair->single is for bonds. Then you avoid having to communicate bond information because with newton on the proc which “knows” which type a bond is can also be the proc in charge of calculating the stress contribution, regardless of which atoms in the bond are local or ghost.

The main thing to check is then that with newton off, local<->ghost bonds aren’t double counted. For that you can copy the if (newton_pair == 0 && j >= nlocal) logic in the pair calculation (replacing newton_pair with newton_bond, and omitting the spatial checks which are required when itag == jtag – since atoms can’t be bonded to themselves).

9

I was using first in my previous example, where I tried to unpack it into atom->num_bond and atom->bond_atom. I figured the purpose of first is to map i to local atom ids. Since in this last example I was unpacking it into my ‘own’ arrays, I thought it was not required.

10

I see, that makes sense, thank you for the explanation.

11

Thank you, that makes sense! I just gave it a try and it seems to work, and it’s more or less as performant as my previous ‘optimization’ using bond_atom:

--- a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp
+++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp
@@ -14,6 +14,7 @@
 #include "compute_stress_cartesian.h"
 
 #include "atom.h"
+#include "bond.h"
 #include "citeme.h"
 #include "comm.h"
 #include "domain.h"
@@ -365,6 +366,38 @@ void ComputeStressCartesian::compute_array()
     }
   }
 
+  // Loop over all bonds
+  for (int i_bond = 0; i_bond < neighbor->nbondlist; i_bond++) {
+    // i == atom1, j == atom2
+    int i = neighbor->bondlist[i_bond][0];
+    int j = neighbor->bondlist[i_bond][1];
+    int type  = neighbor->bondlist[i_bond][2];
+
+    // Skip if one of both atoms is not in group
+    if (!(atom->mask[i] & groupbit)) continue;
+    if (!(atom->mask[j] & groupbit)) continue;
+
+    // if newton_bond is off and atom2 is a ghost atom, only compute this on one processor
+    if (!force->newton_bond && j >= atom->nlocal) {
+      if (tag[i] > tag[j]) {
+        if ((tag[i] + tag[j]) % 2 == 0) continue;
+      } else if (tag[i] < tag[j]) {
+        if ((tag[i] < tag[j]) % 2 == 1) continue;
+      }
+    }
+
+    double dx = x[j][0] - x[i][0];
+    double dy = x[j][1] - x[i][1];
+    double dz = x[j][2] - x[i][2];
+    double rsq = dx*dx + dy*dy + dz*dz;
+    double xi = x[i][dir1] - boxlo[dir1];
+    double yi = x[i][dir2] - boxlo[dir2];
+
+    double fbond;
+    force->bond->single(type, rsq, i, j, fbond);
+    compute_pressure(fbond, xi, yi, dx, dy, dz);
+  }
+
   // normalize pressure
   for (bin = 0; bin < nbins1 * nbins2; bin++) {
     tdens[bin] *= invV;

This has the advantage that it’s simpler, works without the special_bonds hack, and doesn’t require communication (even for newton = on).

I have to test it a bit more but if it all looks good I’ll cook up a PR.

As for the angle and dihedral forces: for my use case, I don’t have those. Granted, while I’m at it, it would be nice for completeness to also implement angle and dihedral forces in stress/cartesian, but I’m not sure if that plays nicely with this algorithm. The algorithm deals with localizing forces between pairs of atoms in certain bins, but once you start looking at angle and dihedral forces, they are not between pairs anymore but between 3 and 4 atoms. I’m not sure how to deal with that…

Thanks once again for all the help!