From e9e9790d6e4a43e3d572323927afc8ee641ef316 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Thu, 15 Sep 2016 09:25:37 -0400
Subject: [PATCH] corrections for pair style lj/sf/dipole/sf and its /omp
 variant

(cherry picked from commit f0c8b2af28a58485f6795cf85c7d88eaafa3e52a)
---
 src/USER-MISC/pair_lj_sf_dipole_sf.cpp    | 35 ++++++++++++++++++-----
 src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp |  6 ++--
 2 files changed, 31 insertions(+), 10 deletions(-)

diff --git a/src/USER-MISC/pair_lj_sf_dipole_sf.cpp b/src/USER-MISC/pair_lj_sf_dipole_sf.cpp
index 072a6c9943..b88f72bac0 100644
--- a/src/USER-MISC/pair_lj_sf_dipole_sf.cpp
+++ b/src/USER-MISC/pair_lj_sf_dipole_sf.cpp
@@ -31,6 +31,8 @@
 
 using namespace LAMMPS_NS;
 
+static int warn_single = 0;
+
 /* ---------------------------------------------------------------------- */
 
 PairLJSFDipoleSF::PairLJSFDipoleSF(LAMMPS *lmp) : Pair(lmp)
@@ -87,7 +89,6 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
   double **torque = atom->torque;
   int *type = atom->type;
   int nlocal = atom->nlocal;
-// The global scaling parameters aren't used anymore
   double *special_coul = force->special_coul;
   double *special_lj = force->special_lj;
   int newton_pair = force->newton_pair;
@@ -355,7 +356,7 @@ void PairLJSFDipoleSF::settings(int narg, char **arg)
 
 void PairLJSFDipoleSF::coeff(int narg, char **arg)
 {
-  if (narg < 4 || narg > 7)
+  if (narg < 4 || narg > 8)
     error->all(FLERR,"Incorrect args for pair coefficients");
   if (!allocated) allocate();
 
@@ -366,13 +367,27 @@ void PairLJSFDipoleSF::coeff(int narg, char **arg)
   double epsilon_one = force->numeric(FLERR,arg[2]);
   double sigma_one = force->numeric(FLERR,arg[3]);
 
-  double scale_one = 1.0;
-  if (narg >= 5) scale_one = force->numeric(FLERR,arg[4]);
-
   double cut_lj_one = cut_lj_global;
   double cut_coul_one = cut_coul_global;
-  if (narg >= 6) cut_coul_one = cut_lj_one = force->numeric(FLERR,arg[5]);
-  if (narg == 7) cut_coul_one = force->numeric(FLERR,arg[6]);
+  double scale_one = 1.0;
+  int iarg = 4;
+
+  if ((narg > iarg) && (strcmp(arg[iarg],"scale") != 0)) {
+    cut_coul_one = cut_lj_one = force->numeric(FLERR,arg[iarg]);
+    ++iarg;
+  }
+  if ((narg > iarg) && (strcmp(arg[iarg],"scale") != 0)) {
+    cut_coul_one = force->numeric(FLERR,arg[iarg]);
+    ++iarg;
+  }
+  if (narg > iarg) {
+    if (strcmp(arg[iarg],"scale") == 0) {
+      scale_one = force->numeric(FLERR,arg[iarg+1]);
+      iarg += 2;
+    } else error->all(FLERR,"Incorrect args for pair coefficients");
+  }
+  if (iarg != narg)
+    error->all(FLERR,"Incorrect args for pair coefficients");
 
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
@@ -532,6 +547,12 @@ double PairLJSFDipoleSF::single(int i, int j, int itype, int jtype, double rsq,
   double *q = atom->q;
   double **mu = atom->mu;
 
+  if (!warn_single) {
+    warn_single = 1;
+    if (comm->me == 0) {
+      error->warning(FLERR,"Single method for lj/sf/dipole/sf does not compute forces");
+    }
+  }
   qtmp = q[i];
   xtmp = x[i][0];
   ytmp = x[i][1];
diff --git a/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp b/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp
index 1eadfe3d1b..29250b09f4 100644
--- a/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp
+++ b/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp
@@ -247,7 +247,7 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
 
         // total force
 
-        fq = factor_coul*qqrd2e;
+        fq = factor_coul*qqrd2e*scale[itype][jtype];
         fx = fq*forcecoulx + delx*forcelj;
         fy = fq*forcecouly + dely*forcelj;
         fz = fq*forcecoulz + delz*forcelj;
@@ -272,7 +272,7 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
 
         if (EFLAG) {
           if (rsq < cut_coulsq[itype][jtype]) {
-            ecoul = (1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype]));
+            ecoul = (1.0-sqrt(rsq/cut_coulsq[itype][jtype]));
             ecoul *= ecoul;
             ecoul *= qtmp * q[j] * rinv;
             if (mu[i].w > 0.0 && mu[j].w > 0.0)
@@ -281,7 +281,7 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
               ecoul += -q[j] * r3inv * pqfac * pidotr;
             if (mu[j].w > 0.0 && qtmp != 0.0)
               ecoul += qtmp * r3inv * qpfac * pjdotr;
-            ecoul *= factor_coul*qqrd2e;
+            ecoul *= factor_coul*qqrd2e*scale[itype][jtype];
           } else ecoul = 0.0;
 
           if (rsq < cut_ljsq[itype][jtype]) {