From cfc8137d2d98c2790b4ec39a8ba1948d882d0a42 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Thu, 8 Sep 2022 22:38:27 -0400
Subject: [PATCH] add check to atom style dielectric code so it will only run
 with dielectric pair styles

---
 src/DIELECTRIC/atom_vec_dielectric.cpp | 40 ++++++++++++++++++++++++++
 src/DIELECTRIC/atom_vec_dielectric.h   |  1 +
 src/pair_hybrid.h                      |  1 +
 3 files changed, 42 insertions(+)

diff --git a/src/DIELECTRIC/atom_vec_dielectric.cpp b/src/DIELECTRIC/atom_vec_dielectric.cpp
index 518faa6e94..8fb4744fa6 100644
--- a/src/DIELECTRIC/atom_vec_dielectric.cpp
+++ b/src/DIELECTRIC/atom_vec_dielectric.cpp
@@ -15,6 +15,10 @@
 
 #include "atom.h"
 #include "citeme.h"
+#include "error.h"
+#include "force.h"
+#include "pair.h"
+#include "pair_hybrid.h"
 
 #include <cmath>
 
@@ -87,6 +91,42 @@ AtomVecDielectric::AtomVecDielectric(LAMMPS *_lmp) : AtomVec(_lmp)
   bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void AtomVecDielectric::init()
+{
+  AtomVec::init();
+
+  // since atom style dielectric modifies the charge q, it will produce incorrect results
+  // with pair styles using coulomb without dielectric support.
+
+  std::string pair_style(force->pair_style);
+  if ((pair_style != "none") && (pair_style != "zero") && !utils::strmatch(force->pair_style,"/dielectric")) {
+    bool mismatch = false;
+    if (utils::strmatch(force->pair_style,"^reaxff")) mismatch = true;
+    if (utils::strmatch(force->pair_style,"^comb")) mismatch = true;
+    if (utils::strmatch(force->pair_style,"coul")) mismatch = true;
+    if (utils::strmatch(force->pair_style,"tip4p")) mismatch = true;
+    if (utils::strmatch(force->pair_style,"dipole")) mismatch = true;
+
+    if (utils::strmatch(force->pair_style,"^hybrid")) {
+      auto hybrid = dynamic_cast<PairHybrid *>(force->pair);
+      if (hybrid) {
+        for (int i = 0; i < hybrid->nstyles; i++) {
+          if (utils::strmatch(hybrid->keywords[i],"^reaxff")) mismatch = true;
+          if (utils::strmatch(hybrid->keywords[i],"^comb")) mismatch = true;
+          if (utils::strmatch(hybrid->keywords[i],"coul")) mismatch = true;
+          if (utils::strmatch(hybrid->keywords[i],"tip4p")) mismatch = true;
+          if (utils::strmatch(hybrid->keywords[i],"dipole")) mismatch = true;
+        }
+      }
+    }
+    if (mismatch)
+      error->all(FLERR, "Pair style {} is not compatible with atom style {}",
+                 pair_style, atom->get_style());
+  }
+}
+
 /* ----------------------------------------------------------------------
    set local copies of all grow ptrs used by this class, except defaults
    needed in replicate when 2 atom classes exist and it calls pack_restart()
diff --git a/src/DIELECTRIC/atom_vec_dielectric.h b/src/DIELECTRIC/atom_vec_dielectric.h
index 91c3b84124..77b6001896 100644
--- a/src/DIELECTRIC/atom_vec_dielectric.h
+++ b/src/DIELECTRIC/atom_vec_dielectric.h
@@ -31,6 +31,7 @@ class AtomVecDielectric : public AtomVec {
  public:
   AtomVecDielectric(class LAMMPS *);
 
+  void init() override;
   void grow_pointers() override;
   void create_atom_post(int) override;
   void data_atom_post(int) override;
diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h
index 8c7782bc4f..a326f0d80c 100644
--- a/src/pair_hybrid.h
+++ b/src/pair_hybrid.h
@@ -25,6 +25,7 @@ PairStyle(hybrid,PairHybrid);
 namespace LAMMPS_NS {
 
 class PairHybrid : public Pair {
+  friend class AtomVecDielectric;
   friend class ComputeSpin;
   friend class FixGPU;
   friend class FixIntel;