fix_dplr: make pppm_dplr optional

source/lmp/fix_dplr.cpp

@@ -175,11 +175,9 @@ void FixDPLR::get_valid_pairs(vector<pair<int, int> > &pairs) {
   int nall = nlocal + nghost;
   vector<int> dtype(nall);
   // get type
-  {
-    int *type = atom->type;
-    for (int ii = 0; ii < nall; ++ii) {
-      dtype[ii] = type[ii] - 1;
-    }
+  int *type = atom->type;
+  for (int ii = 0; ii < nall; ++ii) {
+    dtype[ii] = type[ii] - 1;
   }
 
   int **bondlist = neighbor->bondlist;
@@ -190,21 +188,51 @@ void FixDPLR::get_valid_pairs(vector<pair<int, int> > &pairs) {
     if (!binary_search(bond_type.begin(), bond_type.end(), bd_type)) {
       continue;
     }
-    if (binary_search(sel_type.begin(), sel_type.end(),
-                      dtype[bondlist[ii][0]]) &&
-        binary_search(dpl_type.begin(), dpl_type.end(),
-                      dtype[bondlist[ii][1]])) {
-      idx0 = bondlist[ii][0];
-      idx1 = bondlist[ii][1];
-    } else if (binary_search(sel_type.begin(), sel_type.end(),
-                             dtype[bondlist[ii][1]]) &&
-               binary_search(dpl_type.begin(), dpl_type.end(),
-                             dtype[bondlist[ii][0]])) {
-      idx0 = bondlist[ii][1];
-      idx1 = bondlist[ii][0];
+    std::vector<int>::iterator it =
+        find(sel_type.begin(), sel_type.end(), dtype[bondlist[ii][0]]);
+    if (it != sel_type.end()) {
+      int idx_type = distance(sel_type.begin(), it);
+      if (dtype[bondlist[ii][1]] == dpl_type[idx_type]) {
+        idx0 = bondlist[ii][0];
+        idx1 = bondlist[ii][1];
+      } else {
+        char str[300];
+        sprintf(str,
+                "Invalid pair: %d %d \n       A virtual atom of type %d is "
+                "expected, but the type of atom %d is "
+                "%d.\n       Please check your data file carefully.\n",
+                atom->tag[bondlist[ii][0]], atom->tag[bondlist[ii][1]],
+                dpl_type[idx_type] + 1, atom->tag[bondlist[ii][1]],
+                type[bondlist[ii][1]]);
+        error->all(FLERR, str);
+      }
     } else {
-      error->all(FLERR,
-                 "find a bonded pair the types of which are not associated");
+      it = find(sel_type.begin(), sel_type.end(), dtype[bondlist[ii][1]]);
+      if (it != sel_type.end()) {
+        int idx_type = distance(sel_type.begin(), it);
+        if (dtype[bondlist[ii][0]] == dpl_type[idx_type]) {
+          idx0 = bondlist[ii][1];
+          idx1 = bondlist[ii][0];
+        } else {
+          char str[300];
+          sprintf(str,
+                  "Invalid pair: %d %d \n       A virtual atom of type %d is "
+                  "expected, but the type of atom %d is %d.\n       Please "
+                  "check your data file carefully.\n",
+                  atom->tag[bondlist[ii][0]], atom->tag[bondlist[ii][1]],
+                  dpl_type[idx_type] + 1, atom->tag[bondlist[ii][0]],
+                  type[bondlist[ii][0]]);
+          error->all(FLERR, str);
+        }
+      } else {
+        char str[300];
+        sprintf(str,
+                "Invalid pair: %d %d \n       They are not expected to have "
+                "Wannier centroid.\n       Please check your data file "
+                "carefully.\n",
+                atom->tag[bondlist[ii][0]], atom->tag[bondlist[ii][1]]);
+        error->all(FLERR, str);
+      }
     }
     if (!(idx0 < nlocal && idx1 < nlocal)) {
       error->all(FLERR,
@@ -306,11 +334,6 @@ void FixDPLR::pre_force(int vflag) {
   //   }
   // }
 
-  // selected type
-  vector<int> dpl_type;
-  for (int ii = 0; ii < sel_type.size(); ++ii) {
-    dpl_type.push_back(type_asso[sel_type[ii]]);
-  }
   vector<int> sel_fwd, sel_bwd;
   int sel_nghost;
   deepmd_compat::select_by_type(sel_fwd, sel_bwd, sel_nghost, dcoord, dtype,
@@ -368,10 +391,6 @@ void FixDPLR::post_force(int vflag) {
   }
 
   PPPMDPLR *pppm_dplr = (PPPMDPLR *)force->kspace_match("pppm/dplr", 1);
-  if (!pppm_dplr) {
-    error->all(FLERR, "kspace_style pppm/dplr should be set before this fix\n");
-  }
-  const vector<double> &dfele_(pppm_dplr->get_fele());
   int nlocal = atom->nlocal;
   int nghost = atom->nghost;
   int nall = nlocal + nghost;
@@ -397,10 +416,13 @@ void FixDPLR::post_force(int vflag) {
         dcoord[ii * 3 + dd] = x[ii][dd] - domain->boxlo[dd];
       }
     }
-    assert(dfele_.size() == nlocal * 3);
     // revise force according to efield
-    for (int ii = 0; ii < nlocal * 3; ++ii) {
-      dfele[ii] = dfele_[ii];
+    if (pppm_dplr) {
+      const vector<double> &dfele_(pppm_dplr->get_fele());
+      assert(dfele_.size() == nlocal * 3);
+      for (int ii = 0; ii < nlocal * 3; ++ii) {
+        dfele[ii] += dfele_[ii];
+      }
     }
     // revise force and virial according to efield
     double *q = atom->q;