diff --git a/src/.gitignore b/src/.gitignore
index c26eaaba30..a93abe03f4 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -437,6 +437,8 @@
 /pair_coul_cut_global.h
 /pair_coul_streitz.cpp
 /pair_coul_streitz.h
+/pair_coul_ctip.cpp
+/pair_coul_ctip.h
 /pair_gauss.cpp
 /pair_gauss.h
 /pair_lj96_cut.cpp
@@ -976,6 +978,8 @@
 /fix_pour.h
 /fix_qeq_comb.cpp
 /fix_qeq_comb.h
+/fix_qeq_ctip.cpp
+/fix_qeq_ctip.h
 /fix_qeq_fire.cpp
 /fix_qeq_fire.h
 /fix_qeq_reaxff.cpp
diff --git a/src/EXTRA-PAIR/pair_coul_ctip.cpp b/src/EXTRA-PAIR/pair_coul_ctip.cpp
index e9f47a17a0..0291886987 100644
--- a/src/EXTRA-PAIR/pair_coul_ctip.cpp
+++ b/src/EXTRA-PAIR/pair_coul_ctip.cpp
@@ -20,6 +20,7 @@
 #include "atom.h"
 #include "comm.h"
 #include "error.h"
+#include "ewald_const.h"
 #include "force.h"
 #include "math_const.h"
 #include "memory.h"
@@ -32,21 +33,15 @@
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
+using namespace EwaldConst;
 
-#define EWALD_F 1.12837917
-#define EWALD_P 0.3275911
-#define A1 0.254829592
-#define A2 -0.284496736
-#define A3 1.421413741
-#define A4 -1.453152027
-#define A5 1.061405429
-#define DELTA 4
+static constexpr int DELTA = 4;
 
 /* ---------------------------------------------------------------------- */
 
 PairCoulCTIP::PairCoulCTIP(LAMMPS *lmp) : Pair(lmp)
 {
-params= nullptr;
+  params = nullptr;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -87,8 +82,10 @@ void PairCoulCTIP::compute(int eflag, int vflag)
   double cut_coulcu, cut_coul4, alphacu, erfcd_cut, t_cut, erfcc_cut;
   int elt1, elt2;
 
-  double shield[nelements][nelements], shieldcu[nelements][nelements], reffc[nelements][nelements], reffcsq[nelements][nelements], reffc4[nelements][nelements], reffc7[nelements][nelements];
-  double s2d_shift[nelements][nelements], f_shift[nelements][nelements], e_shift[nelements][nelements], self_factor[nelements][nelements];
+  double shield[nelements][nelements], shieldcu[nelements][nelements], reffc[nelements][nelements],
+      reffcsq[nelements][nelements], reffc4[nelements][nelements], reffc7[nelements][nelements];
+  double s2d_shift[nelements][nelements], f_shift[nelements][nelements],
+      e_shift[nelements][nelements], self_factor[nelements][nelements];
 
   ecoul = 0.0;
   selfion = 0.0;
@@ -123,9 +120,14 @@ void PairCoulCTIP::compute(int eflag, int vflag)
       reffcsq[elt1][elt2] = reffc[elt1][elt2] * reffc[elt1][elt2];
       reffc4[elt1][elt2] = reffcsq[elt1][elt2] * reffcsq[elt1][elt2];
       reffc7[elt1][elt2] = reffc4[elt1][elt2] * reffcsq[elt1][elt2] * reffc[elt1][elt2];
-      s2d_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coulcu + 4.0 * alpha / MY_PIS * erfcd_cut / cut_coulsq + 4.0 * alphacu / MY_PIS * erfcd_cut - 2 / cut_coulcu + 4 * cut_coul4 / reffc7[elt1][elt2] - 2 * cut_coul / reffc4[elt1][elt2];
-      f_shift[elt1][elt2] = erfcc_cut / cut_coulsq + 2.0 * alpha / MY_PIS * erfcd_cut / cut_coul - 1 / cut_coulsq + cut_coulsq / reffc4[elt1][elt2];
-      e_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coul + 2.0 * alpha / MY_PIS * erfcd_cut - 2 / cut_coul + 1 / reffc[elt1][elt2] + cut_coulcu / reffc4[elt1][elt2] + s2d_shift[elt1][elt2] * cut_coulsq * 0.5;
+      s2d_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coulcu +
+          4.0 * alpha / MY_PIS * erfcd_cut / cut_coulsq + 4.0 * alphacu / MY_PIS * erfcd_cut -
+          2 / cut_coulcu + 4 * cut_coul4 / reffc7[elt1][elt2] - 2 * cut_coul / reffc4[elt1][elt2];
+      f_shift[elt1][elt2] = erfcc_cut / cut_coulsq + 2.0 * alpha / MY_PIS * erfcd_cut / cut_coul -
+          1 / cut_coulsq + cut_coulsq / reffc4[elt1][elt2];
+      e_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coul + 2.0 * alpha / MY_PIS * erfcd_cut -
+          2 / cut_coul + 1 / reffc[elt1][elt2] + cut_coulcu / reffc4[elt1][elt2] +
+          s2d_shift[elt1][elt2] * cut_coulsq * 0.5;
       self_factor[elt1][elt2] = -(e_shift[elt1][elt2] * 0.5 + alpha / MY_PIS) * qqrd2e;
     }
   }
@@ -143,9 +145,9 @@ void PairCoulCTIP::compute(int eflag, int vflag)
     jnum = numneigh[i];
     iparam_i = elem1param[itype];
 
-    selfion = self(&params[iparam_i],qtmp);
+    selfion = self(&params[iparam_i], qtmp);
 
-    if (evflag) ev_tally(i,i,nlocal,0,0.0,selfion,0.0,0.0,0.0,0.0);
+    if (evflag) ev_tally(i, i, nlocal, 0, 0.0, selfion, 0.0, 0.0, 0.0, 0.0);
 
     if (eflag) {
       double e_self = self_factor[iparam_i][iparam_i] * qtmp * qtmp;
@@ -166,7 +168,7 @@ void PairCoulCTIP::compute(int eflag, int vflag)
         jtype = map[type[j]];
         jparam_j = elem1param[jtype];
         r = sqrt(rsq);
-        reff = cbrt(rsq * r + 1/shieldcu[iparam_i][jparam_j]);
+        reff = cbrt(rsq * r + 1 / shieldcu[iparam_i][jparam_j]);
         reffsq = reff * reff;
         reff4 = reffsq * reffsq;
         prefactor = qqrd2e * qtmp * q[j] / r;
@@ -174,7 +176,10 @@ void PairCoulCTIP::compute(int eflag, int vflag)
         t = 1.0 / (1.0 + EWALD_P * alpha * r);
         erfcc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * erfcd;
 
-        forcecoul = prefactor * (erfcc / r + 2.0 * alpha / MY_PIS * erfcd + rsq * r / reff4 - 1 / r - r * f_shift[iparam_i][jparam_j] + r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul)) * r;
+        forcecoul = prefactor *
+            (erfcc / r + 2.0 * alpha / MY_PIS * erfcd + rsq * r / reff4 - 1 / r -
+             r * f_shift[iparam_i][jparam_j] + r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul)) *
+            r;
         if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
         fpair = forcecoul / rsq;
 
@@ -188,7 +193,10 @@ void PairCoulCTIP::compute(int eflag, int vflag)
         }
 
         if (eflag) {
-          ecoul = prefactor * (erfcc + r / reff - 1 - r * erfcc_cut / cut_coul - r / reffc[iparam_i][jparam_j] + r / cut_coul + r * f_shift[iparam_i][jparam_j] * (r - cut_coul) - r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul) * (r - cut_coul)* 0.5);
+          ecoul = prefactor *
+              (erfcc + r / reff - 1 - r * erfcc_cut / cut_coul - r / reffc[iparam_i][jparam_j] +
+               r / cut_coul + r * f_shift[iparam_i][jparam_j] * (r - cut_coul) -
+               r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul) * (r - cut_coul) * 0.5);
           if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
         } else
           ecoul = 0.0;
@@ -205,17 +213,17 @@ void PairCoulCTIP::compute(int eflag, int vflag)
 
 double PairCoulCTIP::self(Param *param, double qi)
 {
- double s1=param->chi, s2=param->eta, s3=param->qmin, s4=param->qmax, s5=param->omega;
+  double s1 = param->chi, s2 = param->eta, s3 = param->qmin, s4 = param->qmax, s5 = param->omega;
 
- if ( qi < s3 ) {
-         return qi*((s1-2*s3*s5)+qi*(0.50*s2+s5))+s3*s3*s5;
- } else if ( qi < s4 ) {
-         return qi*(s1+qi*(0.50*s2));
- } else {
-         return qi*((s1-2*s4*s5)+qi*(0.50*s2+s5))+s4*s4*s5;
- }
+  if (qi < s3) {
+    return qi * ((s1 - 2 * s3 * s5) + qi * (0.50 * s2 + s5)) + s3 * s3 * s5;
+  } else if (qi < s4) {
+    return qi * (s1 + qi * (0.50 * s2));
+  } else {
+    return qi * ((s1 - 2 * s4 * s5) + qi * (0.50 * s2 + s5)) + s4 * s4 * s5;
+  }
 
- return 0.0;
+  return 0.0;
 }
 
 /* ----------------------------------------------------------------------
@@ -233,16 +241,16 @@ void PairCoulCTIP::allocate()
 
   memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
 
-  memory->create(qeq_x,n+1,"pair:qeq_x");
-  memory->create(qeq_j,n+1,"pair:qeq_j");
-  memory->create(qeq_g,n+1,"pair:qeq_g");
-  memory->create(qeq_z,n+1,"pair:qeq_z");
-  memory->create(qeq_c,n+1,"pair:qeq_c");
-  memory->create(qeq_q1,n+1,"pair:qeq_q1");
-  memory->create(qeq_q2,n+1,"pair:qeq_q2");
-  memory->create(qeq_w,n+1,"pair:qeq_w");
+  memory->create(qeq_x, n + 1, "pair:qeq_x");
+  memory->create(qeq_j, n + 1, "pair:qeq_j");
+  memory->create(qeq_g, n + 1, "pair:qeq_g");
+  memory->create(qeq_z, n + 1, "pair:qeq_z");
+  memory->create(qeq_c, n + 1, "pair:qeq_c");
+  memory->create(qeq_q1, n + 1, "pair:qeq_q1");
+  memory->create(qeq_q2, n + 1, "pair:qeq_q2");
+  memory->create(qeq_w, n + 1, "pair:qeq_w");
 
-  map = new int[n+1];
+  map = new int[n + 1];
 }
 
 /* ----------------------------------------------------------------------
@@ -265,13 +273,12 @@ void PairCoulCTIP::coeff(int narg, char **arg)
 {
   if (!allocated) allocate();
 
-  map_element2type(narg-3, arg+3);
+  map_element2type(narg - 3, arg + 3);
 
   // read potential file and initialize potential parameters
 
   read_file(arg[2]);
   setup_params();
-
 }
 
 /* ----------------------------------------------------------------------
@@ -285,7 +292,6 @@ void PairCoulCTIP::init_style()
   neighbor->add_request(this);
 
   cut_coulsq = cut_coul * cut_coul;
-
 }
 
 /* ----------------------------------------------------------------------
@@ -309,7 +315,7 @@ void PairCoulCTIP::read_file(char *file)
   // open file on proc 0
   if (comm->me == 0) {
     PotentialFileReader reader(lmp, file, "coul/ctip");
-    char * line;
+    char *line;
 
     while ((line = reader.next_line(NPARAMS_PER_LINE))) {
       try {
@@ -327,16 +333,14 @@ void PairCoulCTIP::read_file(char *file)
 
         if (nparams == maxparam) {
           maxparam += DELTA;
-          params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
-                                              "pair:params");
+          params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
 
           // make certain all addional allocated storage is initialized
           // to avoid false positives when checking with valgrind
 
-          memset(params + nparams, 0, DELTA*sizeof(Param));
+          memset(params + nparams, 0, DELTA * sizeof(Param));
         }
 
-
         params[nparams].ielement = ielement;
         params[nparams].chi = values.next_double();
         params[nparams].eta = values.next_double();
@@ -353,9 +357,10 @@ void PairCoulCTIP::read_file(char *file)
 
       // parameter sanity check
 
-      if (params[nparams].eta < 0.0 || params[nparams].zeta < 0.0 ||
-          params[nparams].zcore < 0.0 || params[nparams].gamma < 0.0 || params[nparams].qmin > params[nparams].qmax || params[nparams].omega < 0.0 )
-        error->one(FLERR,"Illegal coul/ctip parameter");
+      if (params[nparams].eta < 0.0 || params[nparams].zeta < 0.0 || params[nparams].zcore < 0.0 ||
+          params[nparams].gamma < 0.0 || params[nparams].qmin > params[nparams].qmax ||
+          params[nparams].omega < 0.0)
+        error->one(FLERR, "Illegal coul/ctip parameter");
 
       nparams++;
     }
@@ -365,35 +370,34 @@ void PairCoulCTIP::read_file(char *file)
   MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
 
   if (comm->me != 0) {
-    params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
+    params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
   }
 
-  MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
+  MPI_Bcast(params, maxparam * sizeof(Param), MPI_BYTE, 0, world);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void PairCoulCTIP::setup_params()
 {
-  int i,m,n;
+  int i, m, n;
 
   // set elem1param
 
   memory->destroy(elem1param);
-  memory->create(elem1param,nelements,"pair:elem1param");
+  memory->create(elem1param, nelements, "pair:elem1param");
 
   for (i = 0; i < nelements; i++) {
     n = -1;
     for (m = 0; m < nparams; m++) {
       if (i == params[m].ielement) {
-        if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
+        if (n >= 0) error->all(FLERR, "Potential file has duplicate entry");
         n = m;
       }
     }
-    if (n < 0) error->all(FLERR,"Potential file is missing an entry");
+    if (n < 0) error->all(FLERR, "Potential file is missing an entry");
     elem1param[i] = n;
   }
-
 }
 
 /* ----------------------------------------------------------------------
@@ -465,60 +469,76 @@ void *PairCoulCTIP::extract(const char *str, int &dim)
     dim = 0;
     return (void *) &cut_coul;
   }
-  if (strcmp(str,"chi") == 0 && qeq_x) {
+  if (strcmp(str, "chi") == 0 && qeq_x) {
     dim = 1;
     for (int i = 1; i <= atom->ntypes; i++)
-      if (map[i] >= 0) qeq_x[i] = params[map[i]].chi;
-      else qeq_x[i] = 0.0;
+      if (map[i] >= 0)
+        qeq_x[i] = params[map[i]].chi;
+      else
+        qeq_x[i] = 0.0;
     return (void *) qeq_x;
   }
-  if (strcmp(str,"eta") == 0 && qeq_j) {
+  if (strcmp(str, "eta") == 0 && qeq_j) {
     dim = 1;
     for (int i = 1; i <= atom->ntypes; i++)
-      if (map[i] >= 0) qeq_j[i] = params[map[i]].eta;
-      else qeq_j[i] = 0.0;
+      if (map[i] >= 0)
+        qeq_j[i] = params[map[i]].eta;
+      else
+        qeq_j[i] = 0.0;
     return (void *) qeq_j;
   }
-  if (strcmp(str,"gamma") == 0 && qeq_g) {
+  if (strcmp(str, "gamma") == 0 && qeq_g) {
     dim = 1;
     for (int i = 1; i <= atom->ntypes; i++)
-      if (map[i] >= 0) qeq_g[i] = params[map[i]].gamma;
-      else qeq_g[i] = 0.0;
+      if (map[i] >= 0)
+        qeq_g[i] = params[map[i]].gamma;
+      else
+        qeq_g[i] = 0.0;
     return (void *) qeq_g;
   }
-  if (strcmp(str,"zeta") == 0 && qeq_z) {
+  if (strcmp(str, "zeta") == 0 && qeq_z) {
     dim = 1;
     for (int i = 1; i <= atom->ntypes; i++)
-      if (map[i] >= 0) qeq_z[i] = params[map[i]].zeta;
-      else qeq_z[i] = 0.0;
+      if (map[i] >= 0)
+        qeq_z[i] = params[map[i]].zeta;
+      else
+        qeq_z[i] = 0.0;
     return (void *) qeq_z;
   }
-  if (strcmp(str,"zcore") == 0 && qeq_c) {
+  if (strcmp(str, "zcore") == 0 && qeq_c) {
     dim = 1;
     for (int i = 1; i <= atom->ntypes; i++)
-      if (map[i] >= 0) qeq_c[i] = params[map[i]].zcore;
-      else qeq_c[i] = 0.0;
+      if (map[i] >= 0)
+        qeq_c[i] = params[map[i]].zcore;
+      else
+        qeq_c[i] = 0.0;
     return (void *) qeq_c;
   }
-  if (strcmp(str,"qmin") == 0 && qeq_q1) {
+  if (strcmp(str, "qmin") == 0 && qeq_q1) {
     dim = 1;
     for (int i = 1; i <= atom->ntypes; i++)
-      if (map[i] >= 0) qeq_q1[i] = params[map[i]].qmin;
-      else qeq_q1[i] = 0.0;
+      if (map[i] >= 0)
+        qeq_q1[i] = params[map[i]].qmin;
+      else
+        qeq_q1[i] = 0.0;
     return (void *) qeq_q1;
   }
-  if (strcmp(str,"qmax") == 0 && qeq_q2) {
+  if (strcmp(str, "qmax") == 0 && qeq_q2) {
     dim = 1;
     for (int i = 1; i <= atom->ntypes; i++)
-      if (map[i] >= 0) qeq_q2[i] = params[map[i]].qmax;
-      else qeq_q2[i] = 0.0;
+      if (map[i] >= 0)
+        qeq_q2[i] = params[map[i]].qmax;
+      else
+        qeq_q2[i] = 0.0;
     return (void *) qeq_q2;
   }
-  if (strcmp(str,"omega") == 0 && qeq_w) {
+  if (strcmp(str, "omega") == 0 && qeq_w) {
     dim = 1;
     for (int i = 1; i <= atom->ntypes; i++)
-      if (map[i] >= 0) qeq_w[i] = params[map[i]].omega;
-      else qeq_w[i] = 0.0;
+      if (map[i] >= 0)
+        qeq_w[i] = params[map[i]].omega;
+      else
+        qeq_w[i] = 0.0;
     return (void *) qeq_w;
   }
   return nullptr;