diff --git a/doc/src/JPG/zeeman_langevin.jpg b/doc/src/JPG/zeeman_langevin.jpg
index c42c49b6c0..53ce77d2a8 100644
Binary files a/doc/src/JPG/zeeman_langevin.jpg and b/doc/src/JPG/zeeman_langevin.jpg differ
diff --git a/src/SPIN/pair_spin_dipole_cut.cpp b/src/SPIN/pair_spin_dipole_cut.cpp
index a393fe7021..a7372b480d 100644
--- a/src/SPIN/pair_spin_dipole_cut.cpp
+++ b/src/SPIN/pair_spin_dipole_cut.cpp
@@ -252,7 +252,7 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
       if (eflag) {
         if (rsq <= local_cut2) {
           evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
-          evdwl *= hbar;
+          evdwl *= 0.5*hbar;
         }
       } else evdwl = 0.0;
 
diff --git a/src/SPIN/pair_spin_dmi.cpp b/src/SPIN/pair_spin_dmi.cpp
index af68c0d256..04c2dc408d 100644
--- a/src/SPIN/pair_spin_dmi.cpp
+++ b/src/SPIN/pair_spin_dmi.cpp
@@ -319,7 +319,6 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3])
   // if interaction applies to type ii,
   // locflag = 1 and compute pair interaction
 
-  //i = ilist[ii];
   if (locflag == 1) {
 
     xi[0] = x[ii][0];
@@ -373,9 +372,9 @@ void PairSpinDmi::compute_dmi(int i, int j, double eij[3], double fmi[3], double
   dmiy = eij[2]*v_dmx[itype][jtype] - eij[0]*v_dmz[itype][jtype];
   dmiz = eij[0]*v_dmy[itype][jtype] - eij[1]*v_dmx[itype][jtype];
 
-  fmi[0] -= 2.0*(dmiy*spj[2] - dmiz*spj[1]);
-  fmi[1] -= 2.0*(dmiz*spj[0] - dmix*spj[2]);
-  fmi[2] -= 2.0*(dmix*spj[1] - dmiy*spj[0]);
+  fmi[0] -= (dmiy*spj[2] - dmiz*spj[1]);
+  fmi[1] -= (dmiz*spj[0] - dmix*spj[2]);
+  fmi[2] -= (dmix*spj[1] - dmiy*spj[0]);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/SPIN/pair_spin_exchange.cpp b/src/SPIN/pair_spin_exchange.cpp
index 2aa8b99f32..6eacb04ee3 100644
--- a/src/SPIN/pair_spin_exchange.cpp
+++ b/src/SPIN/pair_spin_exchange.cpp
@@ -242,7 +242,8 @@ void PairSpinExchange::compute(int eflag, int vflag)
 
       if (eflag) {
         evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
-        evdwl *= hbar;
+        evdwl *= 0.5*hbar;
+        // evdwl *= hbar;
       } else evdwl = 0.0;
 
       if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
diff --git a/src/SPIN/pair_spin_magelec.cpp b/src/SPIN/pair_spin_magelec.cpp
index 664f8df56c..fabad4ae4d 100644
--- a/src/SPIN/pair_spin_magelec.cpp
+++ b/src/SPIN/pair_spin_magelec.cpp
@@ -251,7 +251,7 @@ void PairSpinMagelec::compute(int eflag, int vflag)
 
       if (eflag) {
         evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
-        evdwl *= hbar;
+        evdwl *= 0.5*hbar;
       } else evdwl = 0.0;
 
       if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
@@ -362,17 +362,17 @@ void PairSpinMagelec::compute_magelec(int i, int j, double eij[3], double fmi[3]
   vy = v_mey[itype][jtype];
   vz = v_mez[itype][jtype];
 
-  meix = vy*eij[2] - vz*eij[1];
-  meiy = vz*eij[0] - vx*eij[2];
-  meiz = vx*eij[1] - vy*eij[0];
+  meix = (vy*eij[2] - vz*eij[1]);
+  meiy = (vz*eij[0] - vx*eij[2]);
+  meiz = (vx*eij[1] - vy*eij[0]);
 
   meix *= ME[itype][jtype];
   meiy *= ME[itype][jtype];
   meiz *= ME[itype][jtype];
 
-  fmi[0] += spj[1]*meiz - spj[2]*meiy;
-  fmi[1] += spj[2]*meix - spj[0]*meiz;
-  fmi[2] += spj[0]*meiy - spj[1]*meix;
+  fmi[0] += (spj[1]*meiz - spj[2]*meiy);
+  fmi[1] += (spj[2]*meix - spj[0]*meiz);
+  fmi[2] += (spj[0]*meiy - spj[1]*meix);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -391,17 +391,17 @@ void PairSpinMagelec::compute_magelec_mech(int i, int j, double fi[3], double sp
   vy = v_mey[itype][jtype];
   vz = v_mez[itype][jtype];
 
-  meix = spi[1]*spi[2] - spi[2]*spj[1];
-  meiy = spi[2]*spi[0] - spi[0]*spj[2];
-  meiz = spi[0]*spi[1] - spi[1]*spj[0];
+  meix = (spi[1]*spi[2] - spi[2]*spj[1]);
+  meiy = (spi[2]*spi[0] - spi[0]*spj[2]);
+  meiz = (spi[0]*spi[1] - spi[1]*spj[0]);
 
   meix *= ME_mech[itype][jtype];
   meiy *= ME_mech[itype][jtype];
   meiz *= ME_mech[itype][jtype];
 
-  fi[0] += meiy*vz - meiz*vy;
-  fi[1] += meiz*vx - meix*vz;
-  fi[2] += meix*vy - meiy*vx;
+  fi[0] += (meiy*vz - meiz*vy);
+  fi[1] += (meiz*vx - meix*vz);
+  fi[2] += (meix*vy - meiy*vx);
 
 }
 
diff --git a/src/SPIN/pair_spin_neel.cpp b/src/SPIN/pair_spin_neel.cpp
index a3114497a6..4a5d453de2 100644
--- a/src/SPIN/pair_spin_neel.cpp
+++ b/src/SPIN/pair_spin_neel.cpp
@@ -260,7 +260,7 @@ void PairSpinNeel::compute(int eflag, int vflag)
 
       if (eflag) {
         evdwl = (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
-        evdwl *= hbar;
+        evdwl *= 0.5*hbar;
       } else evdwl = 0.0;
 
       if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,