Progress towards supporting ghosts in neighbor list

examples/mliap/in.mliap.unified.hippynn.Al.ghostneigh

@@ -24,7 +24,7 @@ mass 1 26.981
 
 # choose potential
 
-pair_style	mliap unified ghostneigh ../../../hippynn/examples/mliap_unified_hippynn_Al.pkl
+pair_style	mliap unified ghostneigh ../../../hippynn/examples/mliap_unified_hippynn_Al_multiple_layers.pkl
 pair_coeff	* * Al
 
 # Setup output

src/ML-IAP/mliap_data.cpp

@@ -38,7 +38,6 @@ MLIAPData::MLIAPData(LAMMPS *lmp, int gradgradflag_in, int *map_in,
   jatoms(nullptr), jelems(nullptr), elems(nullptr), rij(nullptr),
   graddesc(nullptr), model(nullptr), descriptor(nullptr), list(nullptr)
 {
-  f = atom->f;
   gradgradflag = gradgradflag_in;
   map = map_in;
   model = model_in;
@@ -111,13 +110,10 @@ void MLIAPData::init()
 void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_in)
 {
   list = list_in;
+  f = atom->f;
   double **x = atom->x;
   int *type = atom->type;
 
-  nlocal = atom->nlocal;
-  nghost = atom->nghost;
-  ntotal = nlocal + nghost;
-
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
@@ -134,6 +130,7 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i
   // clear gradforce, elems arrays
 
   int nall = atom->nlocal + atom->nghost;
+  ntotal = nall;
   for (int i = 0; i < nall; i++) {
     for (int j = 0; j < size_gradforce; j++) {
       gradforce[i][j] = 0.0;
@@ -151,6 +148,11 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i
     nlistatoms_max = nlistatoms;
   }
   
+  if (list->ghost == 1)
+    nlistghosts = list->gnum;
+  else
+    nlistghosts = 0;
+
   // grow gamma arrays if necessary
 
   if (gradgradflag == 1) {

src/ML-IAP/mliap_data.h

@@ -57,9 +57,10 @@ class MLIAPData : protected Pointers {
   // data structures for mliap neighbor list
   // only neighbors strictly inside descriptor cutoff
 
-  int nlocal, nghost, ntotal;    // # of owned and ghost atoms on this proc, and sum of both
+  int ntotal;                    // total # of owned and ghost atoms on this proc
   int nlistatoms;                // current number of atoms in neighborlist
   int nlistatoms_max;            // allocated size of descriptor array
+  int nlistghosts;               // current number of ghost atoms in neighborlist
   int natomneigh_max;            // allocated size of atom neighbor arrays
   int *numneighs;                // neighbors count for each atom
   int *iatoms;                   // index of each atom

src/ML-IAP/mliap_unified.cpp

@@ -171,7 +171,7 @@ MLIAPBuildUnified_t LAMMPS_NS::build_unified(char *unified_fname, MLIAPData *dat
 
   PyObject *unified_module = PyImport_ImportModule("mliap_unifiedpy");
 
-  if (!unified_module) {
+  if (PyErr_Occurred()) {
     PyErr_Print();
     PyErr_Clear();
     PyGILState_Release(gstate);

src/ML-IAP/mliap_unifiedpy.pyx

@@ -52,11 +52,9 @@ cdef extern from "mliap_data.h" namespace "LAMMPS_NS":
         # data structures for mliap neighbor list
         # only neighbors strictly inside descriptor cutoff
 
-        int nlocal
+        int ntotal              # total # of owned and ghost atoms on this proc
-        int nghost
-        int ntotal
-        int * elems
         int nlistatoms          # current number of atoms in neighborlist
+        int nlistghosts         # current number of ghost atoms in neighborlist
         int * numneighs         # neighbors count for each atom
         int * iatoms            # index of each atom
         int * pair_i            # index of each i atom for each ij pair
@@ -65,6 +63,7 @@ cdef extern from "mliap_data.h" namespace "LAMMPS_NS":
         int npairs              # number of ij neighbor pairs
         int * jatoms            # index of each neighbor
         int * jelems            # element of each neighbor
+        int * elems             # element of each atom in or not in the neighborlist
         double ** rij           # distance vector of each neighbor
         # ----- write only -----
         double *** graddesc     # descriptor gradient w.r.t. each neighbor
@@ -75,7 +74,7 @@ cdef extern from "mliap_data.h" namespace "LAMMPS_NS":
 
 cdef extern from "mliap_unified.h" namespace "LAMMPS_NS":
     cdef cppclass MLIAPDummyDescriptor:
-        MLIAPDummyDescriptor(PyObject *, LAMMPS *)
+        MLIAPDummyDescriptor(PyObject *, LAMMPS *) except +
         int ndescriptors    # number of descriptors
         int nelements       # # of unique elements
         char **elements     # names of unique elements
@@ -88,15 +87,15 @@ cdef extern from "mliap_unified.h" namespace "LAMMPS_NS":
         void set_elements(char **, int)
 
     cdef cppclass MLIAPDummyModel:
-        MLIAPDummyModel(PyObject *, LAMMPS *, char * = NULL)
+        MLIAPDummyModel(PyObject *, LAMMPS *, char * = NULL) except +
         int ndescriptors    # number of descriptors
         int nparams         # number of parameters per element
         int nelements;      # # of unique elements
 
         void compute_gradients(MLIAPData *)
 
-    cdef void update_pair_energy(MLIAPData *, double *)
+    cdef void update_pair_energy(MLIAPData *, double *) except +
-    cdef void update_pair_forces(MLIAPData *, double *)
+    cdef void update_pair_forces(MLIAPData *, double *) except +
 
 
 # @property sans getter
@@ -209,14 +208,6 @@ cdef class MLIAPDataPy:
     # data structures for mliap neighbor list
     # only neighbors strictly inside descriptor cutoff
 
-    @property
-    def nlocal(self):
-        return self.data.nlocal
-
-    @property
-    def nghost(self):
-        return self.data.nghost
-
     @property
     def ntotal(self):
         return self.data.ntotal
@@ -231,6 +222,10 @@ cdef class MLIAPDataPy:
     def nlistatoms(self):
         return self.data.nlistatoms
     
+    @property
+    def nlistghosts(self):
+        return self.data.nlistghosts
+
     @property
     def numneighs(self):
         if self.data.numneighs is NULL:
@@ -318,19 +313,19 @@ cdef class MLIAPUnifiedInterface:
         self.unified_impl.compute_forces(data)
 
 
-cdef public void compute_gradients_python(unified_int, MLIAPData *data) with gil:
+cdef public void compute_gradients_python(unified_int, MLIAPData *data) except * with gil:
     pydata = MLIAPDataPy()
     pydata.data = data
     unified_int.compute_gradients(pydata)
 
 
-cdef public void compute_descriptors_python(unified_int, MLIAPData *data) with gil:
+cdef public void compute_descriptors_python(unified_int, MLIAPData *data) except * with gil:
     pydata = MLIAPDataPy()
     pydata.data = data
     unified_int.compute_descriptors(pydata)
 
 
-cdef public void compute_forces_python(unified_int, MLIAPData *data) with gil:
+cdef public void compute_forces_python(unified_int, MLIAPData *data) except * with gil:
     pydata = MLIAPDataPy()
     pydata.data = data
     unified_int.compute_forces(pydata)
@@ -361,8 +356,10 @@ cdef public object mliap_unified_connect(char *fname, MLIAPDummyModel * model,
     
     cdef int nelements = <int>len(unified.element_types)
     cdef char **elements = <char**>malloc(nelements * sizeof(char*))
+
     if not elements:
         raise MemoryError("failed to allocate memory for element names")
+
     cdef char *elem_name
     for i, elem in enumerate(unified.element_types):
         elem_name_bytes = elem.encode('UTF-8')

src/ML-IAP/pair_mliap.cpp

@@ -336,7 +336,7 @@ void PairMLIAP::init_style()
 
   // need a full neighbor list
 
-  if (ghostneigh)
+  if (ghostneigh == 1)
     neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_GHOST);
   else
     neighbor->add_request(this, NeighConst::REQ_FULL);