maint(kokkosNL): nicer failure modes for large cutoffs in dense systems

examples/PACKAGES/metatomic/in.kokkos.pair_metatomic

@@ -20,6 +20,10 @@ velocity all create 123 42
 pair_style metatomic/kk nickel-lj.pt
 pair_coeff * * 28
 
+# For dense systems with large cutoffs (e.g. PET-MAD, ~5.5 A interaction range),
+# pair_style metatomic will auto-adjust binsize, one, and page. To override:
+# neigh_modify one 100000 page 1000000 binsize 5.5
+
 # pair_style hybrid/overlay/kk metatomic_1/kk nickel-lj.pt metatomic_2/kk nickel-lj.pt
 # pair_coeff * * metatomic_1/kk 28
 # pair_coeff * * metatomic_2/kk 28

src/KOKKOS/nbin_kokkos.cpp

@@ -17,6 +17,7 @@
 #include "atom_kokkos.h"
 #include "atom_masks.h"
 #include "comm.h"
+#include "error.h"
 #include "kokkos.h"
 #include "memory_kokkos.h"
 #include "update.h"
@@ -63,6 +64,10 @@ NBinKokkos<DeviceType>::NBinKokkos(LAMMPS *lmp) : NBinStandard(lmp) {
 template<class DeviceType>
 void NBinKokkos<DeviceType>::bin_atoms_setup(int nall)
 {
+  if (mbins <= 0)
+    error->one(FLERR, "Kokkos neighbor list produced zero bins; "
+        "try setting neigh_modify binsize to a smaller value");
+
   if (mbins > (int)k_bins.view_device().extent(0)) {
     MemoryKokkos::realloc_kokkos(k_bins,"Neighbor::d_bins",mbins,atoms_per_bin);
     bins = k_bins.view<DeviceType>();

src/KOKKOS/npair_kokkos.cpp

@@ -16,6 +16,7 @@
 #include "atom_kokkos.h"
 #include "atom_masks.h"
 #include "domain_kokkos.h"
+#include "error.h"
 #include "update.h"
 #include "neighbor_kokkos.h"
 #include "nbin_kokkos.h"
@@ -234,6 +235,10 @@ void NPairKokkos<DeviceType,HALF,NEWTON,GHOST,TRI,SIZE>::build(NeighList *list_)
   data.special_flag[2] = special_flag[2];
   data.special_flag[3] = special_flag[3];
 
+  if (atoms_per_bin <= 0)
+    error->one(FLERR, "Kokkos neighbor list bin size produced zero atoms_per_bin; "
+        "try increasing neigh_modify binsize");
+
   data.h_resize()=1;
   while (data.h_resize()) {
     data.h_new_maxneighs() = list->maxneighs;

src/ML-METATOMIC/fix_metatomic.cpp

@@ -39,9 +39,11 @@
 #include "neigh_list.h"
 #include "neigh_request.h"
 #include "comm.h"
+#include "domain.h"
 
 #include <vector>
 #include <algorithm>
+#include <cmath>
 
 #include <metatomic/torch.hpp>
 #include <metatensor/torch.hpp>
@@ -337,16 +339,31 @@ void FixMetatomic::init() {
         this->system_adaptor->add_nl_request(cutoff, options);
     }
 
-    // HACK: Explicitly set the binsize for the neighbor list if there is no
-    // pair_style that would set it instead.
-    //
-    // Otherwise, the default binsize of box[0] is used, which crashes kokkos
-    // for large-ish boxes (~40A), and slow down the simulation for non-kokkos.
-    if (strcmp(force->pair_style, "none") == 0) {
+    // Check that neighbor list parameters are sufficient for this cutoff.
+    // Dense systems with large cutoffs can overflow the default one/page
+    // and crash the Kokkos NL builder with SIGFPE.
+    if (!neighbor->binsizeflag) {
+        // Keep the existing binsize hack for fix metatomic
         neighbor->binsize_user = 0.5 * mta_data->max_cutoff;
         neighbor->binsizeflag = 1;
     }
-    // END HACK
+
+    double volume = domain->xprd * domain->yprd * domain->zprd;
+    double density = (volume > 0) ? static_cast<double>(atom->natoms) / volume : 0.0;
+    double cutoff_with_skin = mta_data->max_cutoff + neighbor->skin;
+    int est_neighbors = static_cast<int>(
+        (4.0/3.0) * M_PI * pow(cutoff_with_skin, 3) * density * 2.0
+    );
+    if (est_neighbors > neighbor->oneatom) {
+        error->one(FLERR,
+            "The metatomic model cutoff ({:.4f}) with current system density "
+            "requires approximately {} neighbors per atom, but neigh_modify one "
+            "is only {}. Add 'neigh_modify one {} page {} binsize {:.4f}' "
+            "to your input script.",
+            mta_data->max_cutoff, est_neighbors, neighbor->oneatom,
+            est_neighbors, est_neighbors * 10,
+            0.5 * mta_data->max_cutoff);
+    }
 }
 
 void FixMetatomic::pick_device(c10::Device& device, const char* requested) {

src/ML-METATOMIC/pair_metatomic.cpp

@@ -38,6 +38,7 @@
 #endif
 
 #include <memory>
+#include <cmath>
 
 #include <metatensor/torch.hpp>
 #include <metatomic/torch.hpp>
@@ -559,6 +560,26 @@ void PairMetatomic::init_style() {
     auto request = neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_GHOST);
     request->set_cutoff(mta_data->max_cutoff);
 
+    // Check that neighbor list parameters are sufficient for this cutoff.
+    // Dense systems with large cutoffs can overflow the default one/page
+    // and crash the Kokkos NL builder with SIGFPE.
+    double volume = compute_volume(domain);
+    double density = (volume > 0) ? static_cast<double>(atom->natoms) / volume : 0.0;
+    double cutoff_with_skin = mta_data->max_cutoff + neighbor->skin;
+    int est_neighbors = static_cast<int>(
+        (4.0/3.0) * M_PI * pow(cutoff_with_skin, 3) * density * 2.0
+    );
+    if (est_neighbors > neighbor->oneatom) {
+        error->one(FLERR,
+            "The metatomic model cutoff ({:.4f}) with current system density "
+            "requires approximately {} neighbors per atom, but neigh_modify one "
+            "is only {}. Add 'neigh_modify one {} page {} binsize {:.4f}' "
+            "to your input script.",
+            mta_data->max_cutoff, est_neighbors, neighbor->oneatom,
+            est_neighbors, est_neighbors * 10,
+            0.5 * mta_data->max_cutoff);
+    }
+
     // Translate from the metatomic neighbor lists requests to LAMMPS neighbor
     // lists requests.
     auto requested_nl = mta_data->model->run_method("requested_neighbor_lists");