Merge pull request lammps#4301 from stanmoore1/reaxff_overflow_errors

Gracefully error out if integer overflow in ReaxFF or QEq

Author: akohlmey

src/KOKKOS/fix_qeq_reaxff_kokkos.cpp

@@ -301,7 +301,7 @@ void FixQEqReaxFFKokkos<DeviceType>::pre_force(int /*vflag*/)
 
 template<class DeviceType>
 KOKKOS_INLINE_FUNCTION
-void FixQEqReaxFFKokkos<DeviceType>::num_neigh_item(int ii, int &maxneigh) const
+void FixQEqReaxFFKokkos<DeviceType>::num_neigh_item(int ii, bigint &maxneigh) const
 {
   const int i = d_ilist[ii];
   maxneigh += d_numneigh[i];
@@ -316,13 +316,16 @@ void FixQEqReaxFFKokkos<DeviceType>::allocate_matrix()
 
   // determine the total space for the H matrix
 
-  m_cap = 0;
+  bigint m_cap_big = 0;
 
   // limit scope of functor to allow deallocation of views
   {
     FixQEqReaxFFKokkosNumNeighFunctor<DeviceType> neigh_functor(this);
-    Kokkos::parallel_reduce(nn,neigh_functor,m_cap);
+    Kokkos::parallel_reduce(nn,neigh_functor,m_cap_big);
   }
+  if (m_cap_big > MAXSMALLINT)
+    error->one(FLERR,"Too many neighbors in fix qeq/reaxff");
+  m_cap = m_cap_big;
 
   // deallocate first to reduce memory overhead
 

src/KOKKOS/fix_qeq_reaxff_kokkos.h

@@ -70,7 +70,7 @@ class FixQEqReaxFFKokkos : public FixQEqReaxFF, public KokkosBase {
   void pre_force(int) override;
 
   KOKKOS_INLINE_FUNCTION
-  void num_neigh_item(int, int&) const;
+  void num_neigh_item(int, bigint&) const;
 
   KOKKOS_INLINE_FUNCTION
   void operator()(TagQEqZero, const int&) const;
@@ -290,13 +290,13 @@ class FixQEqReaxFFKokkos : public FixQEqReaxFF, public KokkosBase {
 template <class DeviceType>
 struct FixQEqReaxFFKokkosNumNeighFunctor {
   typedef DeviceType device_type;
-  typedef int value_type;
+  typedef bigint value_type;
   FixQEqReaxFFKokkos<DeviceType> c;
   FixQEqReaxFFKokkosNumNeighFunctor(FixQEqReaxFFKokkos<DeviceType>* c_ptr):c(*c_ptr) {
     c.cleanup_copy();
   };
   KOKKOS_INLINE_FUNCTION
-  void operator()(const int ii, int &maxneigh) const {
+  void operator()(const int ii, bigint &maxneigh) const {
     c.num_neigh_item(ii, maxneigh);
   }
 };

src/KOKKOS/kokkos.cpp

@@ -638,10 +638,10 @@ void KokkosLMP::accelerator(int narg, char **arg)
    called by Finish
 ------------------------------------------------------------------------- */
 
-int KokkosLMP::neigh_count(int m)
+bigint KokkosLMP::neigh_count(int m)
 {
   int inum = 0;
-  int nneigh = 0;
+  bigint nneigh = 0;
 
   ArrayTypes<LMPHostType>::t_int_1d h_ilist;
   ArrayTypes<LMPHostType>::t_int_1d h_numneigh;

src/KOKKOS/kokkos.h

@@ -64,7 +64,7 @@ class KokkosLMP : protected Pointers {
   static void initialize(const Kokkos::InitializationSettings&, Error *);
   static void finalize();
   void accelerator(int, char **);
-  int neigh_count(int);
+  bigint neigh_count(int);
 
   template<class DeviceType>
   int need_dup(int qeq_flag = 0)

src/KOKKOS/pair_reaxff_kokkos.cpp

@@ -1504,10 +1504,18 @@ void PairReaxFFKokkos<DeviceType>::allocate_array()
   if (cut_hbsq > 0.0) {
     MemKK::realloc_kokkos(d_hb_first,"reaxff/kk:hb_first",nmax);
     MemKK::realloc_kokkos(d_hb_num,"reaxff/kk:hb_num",nmax);
+
+    if (((bigint) nmax*maxhb) > MAXSMALLINT)
+      error->one(FLERR,"Too many hydrogen bonds in pair reaxff");
+
     MemKK::realloc_kokkos(d_hb_list,"reaxff/kk:hb_list",nmax*maxhb);
   }
   MemKK::realloc_kokkos(d_bo_first,"reaxff/kk:bo_first",nmax);
   MemKK::realloc_kokkos(d_bo_num,"reaxff/kk:bo_num",nmax);
+
+  if (((bigint) nmax*maxbo) > MAXSMALLINT)
+    error->one(FLERR,"Too many bonds in pair reaxff");
+
   MemKK::realloc_kokkos(d_bo_list,"reaxff/kk:bo_list",nmax*maxbo);
 
   MemKK::realloc_kokkos(d_BO,"reaxff/kk:BO",nmax,maxbo);

src/QEQ/fix_qeq.cpp

@@ -237,8 +237,9 @@ void FixQEq::reallocate_storage()
 
 void FixQEq::allocate_matrix()
 {
-  int i,ii,inum,m;
+  int i,ii,inum;
   int *ilist, *numneigh;
+  bigint m;
 
   int mincap;
   double safezone;
@@ -261,7 +262,10 @@ void FixQEq::allocate_matrix()
     i = ilist[ii];
     m += numneigh[i];
   }
-  m_cap = MAX((int)(m * safezone), mincap * MIN_NBRS);
+  bigint m_cap_big = (bigint)MAX(m * safezone, mincap * MIN_NBRS);
+  if (m_cap_big > MAXSMALLINT)
+    error->one(FLERR,"Too many neighbors in fix qeq");
+  m_cap = m_cap_big;
 
   H.n = n_cap;
   H.m = m_cap;

src/REAXFF/fix_qeq_reaxff.cpp

@@ -338,7 +338,8 @@ void FixQEqReaxFF::reallocate_storage()
 
 void FixQEqReaxFF::allocate_matrix()
 {
-  int i,ii,m;
+  int i,ii;
+  bigint m;
 
   int mincap;
   double safezone;
@@ -360,7 +361,10 @@ void FixQEqReaxFF::allocate_matrix()
     i = ilist[ii];
     m += numneigh[i];
   }
-  m_cap = MAX((int)(m * safezone), mincap * REAX_MIN_NBRS);
+  bigint m_cap_big = (bigint)MAX(m * safezone, mincap * REAX_MIN_NBRS);
+  if (m_cap_big > MAXSMALLINT)
+    error->one(FLERR,"Too many neighbors in fix qeq/reaxff");
+  m_cap = m_cap_big;
 
   H.n = n_cap;
   H.m = m_cap;

src/REAXFF/reaxff_allocate.cpp

@@ -169,16 +169,23 @@ namespace ReaxFF {
   static int Reallocate_HBonds_List(reax_system *system, reax_list *hbonds)
   {
     int i, total_hbonds;
+    LAMMPS_NS::bigint total_hbonds_big;
 
     int mincap = system->mincap;
     double saferzone = system->saferzone;
 
-    total_hbonds = 0;
+    total_hbonds_big = 0;
     for (i = 0; i < system->n; ++i)
       if ((system->my_atoms[i].Hindex) >= 0) {
-        total_hbonds += system->my_atoms[i].num_hbonds;
+        total_hbonds_big += system->my_atoms[i].num_hbonds;
       }
-    total_hbonds = (int)(MAX(total_hbonds*saferzone, mincap*system->minhbonds));
+    total_hbonds_big = (LAMMPS_NS::bigint)(MAX(total_hbonds_big*saferzone, mincap*system->minhbonds));
+
+    auto error = system->error_ptr;
+    if (total_hbonds_big > MAXSMALLINT)
+      error->one(FLERR,"Too many hydrogen bonds in pair reaxff");
+
+    total_hbonds = total_hbonds_big;
 
     Delete_List(hbonds);
     Make_List(system->Hcap, total_hbonds, TYP_HBOND, hbonds);
@@ -190,17 +197,24 @@ namespace ReaxFF {
                                     reax_list *bonds, int *total_bonds, int *est_3body)
   {
     int i;
+    LAMMPS_NS::bigint total_bonds_big;
 
     int mincap = system->mincap;
     double safezone = system->safezone;
 
-    *total_bonds = 0;
+    total_bonds_big = 0;
     *est_3body = 0;
     for (i = 0; i < system->N; ++i) {
       *est_3body += SQR(system->my_atoms[i].num_bonds);
-      *total_bonds += system->my_atoms[i].num_bonds;
+      total_bonds_big += system->my_atoms[i].num_bonds;
     }
-    *total_bonds = (int)(MAX(*total_bonds * safezone, mincap*MIN_BONDS));
+    total_bonds_big = (LAMMPS_NS::bigint)(MAX(total_bonds_big * safezone, mincap*MIN_BONDS));
+
+    auto error = system->error_ptr;
+    if (total_bonds_big > MAXSMALLINT)
+      error->one(FLERR,"Too many bonds in pair reaxff");
+
+    *total_bonds = total_bonds_big;
 
     if (system->omp_active)
       for (i = 0; i < bonds->num_intrs; ++i)