diff --git a/doc/src/Eqs/fix_spin_zeeman.jpg b/doc/src/Eqs/fix_spin_zeeman.jpg
new file mode 100644
index 0000000000..85cd362e2f
Binary files /dev/null and b/doc/src/Eqs/fix_spin_zeeman.jpg differ
diff --git a/doc/src/Eqs/force_spin_zeeman.tex b/doc/src/Eqs/fix_spin_zeeman.tex
similarity index 70%
rename from doc/src/Eqs/force_spin_zeeman.tex
rename to doc/src/Eqs/fix_spin_zeeman.tex
index d623ef85ae..79dd49526a 100644
--- a/doc/src/Eqs/force_spin_zeeman.tex
+++ b/doc/src/Eqs/fix_spin_zeeman.tex
@@ -5,7 +5,8 @@
 \begin{document}
 \begin{varwidth}{50in}
   \begin{equation}
-    \bm{H}_{zeeman} = -\mu_{B}\mu_0\sum_{i=0}^{N}g_{i} \vec{s}_{i} \cdot \vec{H}_{ext} \nonumber 
+    \bm{H}_{Zeeman} = -g \sum_{i=0}^{N}\mu_{i}\, 
+    \vec{s}_{i} \cdot\vec{B}_{ext} \nonumber
   \end{equation}
 \end{varwidth}
 \end{document}
diff --git a/doc/src/Eqs/force_spin_zeeman.jpg b/doc/src/Eqs/force_spin_zeeman.jpg
deleted file mode 100644
index 14fb5500cf..0000000000
Binary files a/doc/src/Eqs/force_spin_zeeman.jpg and /dev/null differ
diff --git a/doc/src/JPG/zeeman_langevin.jpg b/doc/src/JPG/zeeman_langevin.jpg
new file mode 100644
index 0000000000..c42c49b6c0
Binary files /dev/null and b/doc/src/JPG/zeeman_langevin.jpg differ
diff --git a/doc/src/fix_precession_spin.txt b/doc/src/fix_precession_spin.txt
index 040a3086d3..152f717155 100644
--- a/doc/src/fix_precession_spin.txt
+++ b/doc/src/fix_precession_spin.txt
@@ -41,10 +41,37 @@ Style {zeeman} is used for the simulation of the interaction
 between the magnetic spins in the defined group and an external
 magnetic field:
 
-:c,image(Eqs/force_spin_zeeman.jpg)
+:c,image(Eqs/fix_spin_zeeman.jpg)
 
-with mu0 the vacuum permeability, muB the Bohr magneton (muB = 5.788 eV/T
-in metal units).
+with:
+
+Bext the external magnetic field (in T) :ulb,l
+g the Lande factor (hard-coded as g=2.0) :l
+si the unitary vector describing the orientation of spin i :l
+mui the atomic moment of spin i given as a multiple of the
+Bohr magneton muB (for example, mui ~ 2.2 in bulk iron).  :l
+:ule
+
+The field value in Tesla is multiplied by the gyromagnetic
+ratio, g*muB/hbar, converting it into a precession frequency in
+rad.THz (in metal units and with muB = 5.788 eV/T).
+
+As a comparison, the figure below displays the simulation of a
+single spin (of norm mui = 1.0) submitted to an external
+magnetic field of |Bext| = 10.0 Tesla (and oriented along the z
+axis). 
+The upper plot shows the average magnetization along the
+external magnetic field axis and the lower plot the Zeeman
+energy, both as a function of temperature.
+The reference result is provided by the plot of the Langevin 
+function for the same parameters.
+
+:c,image(JPG/zeeman_langevin.jpg)
+
+The temperature effects are accounted for by connecting the spin
+i to a thermal bath using a Langevin thermostat (see 
+"fix_langevin_spin"_fix_langevin_spin.html for the definition of 
+this thermostat).
 
 Style {anisotropy} is used to simulate an easy axis or an easy plane
 for the magnetic spins in the defined group:
diff --git a/src/SPIN/fix_precession_spin.cpp b/src/SPIN/fix_precession_spin.cpp
index 3296b28228..3b8817704d 100644
--- a/src/SPIN/fix_precession_spin.cpp
+++ b/src/SPIN/fix_precession_spin.cpp
@@ -173,9 +173,9 @@ int FixPrecessionSpin::setmask()
 
 void FixPrecessionSpin::init()
 {
-  const double hbar = force->hplanck/MY_2PI;    // eV/(rad.THz)
-  const double mub = 5.78901e-5;                // in eV/T
-  const double gyro = mub/hbar;                 // in rad.THz/T
+  const double hbar = force->hplanck/MY_2PI;	// eV/(rad.THz)
+  const double mub = 5.78901e-5;		// in eV/T
+  const double gyro = 2.0*mub/hbar;		// in rad.THz/T
 
   // convert field quantities to rad.THz
 
@@ -240,6 +240,7 @@ void FixPrecessionSpin::post_force(int /* vflag */)
   }
 
   int *mask = atom->mask;
+  double *emag = atom->emag;
   double **fm = atom->fm;
   double **sp = atom->sp;
   const int nlocal = atom->nlocal;
@@ -257,7 +258,7 @@ void FixPrecessionSpin::post_force(int /* vflag */)
 
       if (zeeman_flag) {          // compute Zeeman interaction
         compute_zeeman(i,fmi);
-        epreci -= hbar*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
+        epreci -= 2.0*hbar*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
       }
 
       if (aniso_flag) {           // compute magnetic anisotropy
@@ -271,6 +272,7 @@ void FixPrecessionSpin::post_force(int /* vflag */)
       }
 
       eprec += epreci;
+      emag[i] += epreci;
       fm[i][0] += fmi[0];
       fm[i][1] += fmi[1];
       fm[i][2] += fmi[2];
@@ -295,9 +297,9 @@ void FixPrecessionSpin::compute_single_precession(int i, double spi[3], double f
 void FixPrecessionSpin::compute_zeeman(int i, double fmi[3])
 {
   double **sp = atom->sp;
-  fmi[0] += sp[i][3]*hx;
-  fmi[1] += sp[i][3]*hy;
-  fmi[2] += sp[i][3]*hz;
+  fmi[0] += 0.5*sp[i][3]*hx;
+  fmi[1] += 0.5*sp[i][3]*hy;
+  fmi[2] += 0.5*sp[i][3]*hz;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -305,9 +307,9 @@ void FixPrecessionSpin::compute_zeeman(int i, double fmi[3])
 void FixPrecessionSpin::compute_anisotropy(double spi[3], double fmi[3])
 {
   double scalar = nax*spi[0] + nay*spi[1] + naz*spi[2];
-  fmi[0] += scalar*Kax;
-  fmi[1] += scalar*Kay;
-  fmi[2] += scalar*Kaz;
+  fmi[0] += 0.5*scalar*Kax;
+  fmi[1] += 0.5*scalar*Kay;
+  fmi[2] += 0.5*scalar*Kaz;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -316,52 +318,10 @@ double FixPrecessionSpin::compute_anisotropy_energy(double spi[3])
 {
   double energy = 0.0;
   double scalar = nax*spi[0] + nay*spi[1] + naz*spi[2];
-  energy = Ka*scalar*scalar;
+  energy = 2.0*Ka*scalar*scalar;
   return energy; 
 }
 
-/* ---------------------------------------------------------------------- */
-
-void FixPrecessionSpin::post_force_respa(int vflag, int ilevel, int /*iloop*/)
-{
-  if (ilevel == ilevel_respa) post_force(vflag);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixPrecessionSpin::set_magneticprecession()
-{
-  if (zeeman_flag) {
-    hx = H_field*nhx;
-    hy = H_field*nhy;
-    hz = H_field*nhz;
-  }
-  if (aniso_flag) {
-    Kax = 2.0*Kah*nax;
-    Kay = 2.0*Kah*nay;
-    Kaz = 2.0*Kah*naz;
-  }
-}
-
-/* ----------------------------------------------------------------------
-   compute cubic aniso energy of spin i
-------------------------------------------------------------------------- */
-
-double FixPrecessionSpin::compute_cubic_energy(double spi[3])
-{
-  double energy = 0.0;
-  double skx,sky,skz;
-
-  skx = spi[0]*nc1x+spi[1]*nc1y+spi[2]*nc1z;
-  sky = spi[0]*nc2x+spi[1]*nc2y+spi[2]*nc2z;
-  skz = spi[0]*nc3x+spi[1]*nc3y+spi[2]*nc3z;
-
-  energy = k1c*(skx*skx*sky*sky + sky*sky*skz*skz + skx*skx*skz*skz);
-  energy += k2c*skx*skx*sky*sky*skz*skz;
-
-  return energy;
-}
-
 /* ----------------------------------------------------------------------
    compute cubic anisotropy interaction for spin i
 ------------------------------------------------------------------------- */
@@ -396,9 +356,44 @@ void FixPrecessionSpin::compute_cubic(double spi[3], double fmi[3])
   sixy = k2ch*(nc1y*six1 + nc2y*six2 + nc3y*six3);
   sixz = k2ch*(nc1z*six1 + nc2z*six2 + nc3z*six3);
   
-  fmi[0] += fourx + sixx;
-  fmi[1] += foury + sixy;
-  fmi[2] += fourz + sixz;
+  fmi[0] += 0.5*(fourx + sixx);
+  fmi[1] += 0.5*(foury + sixy);
+  fmi[2] += 0.5*(fourz + sixz);
+}
+
+/* ----------------------------------------------------------------------
+   compute cubic aniso energy of spin i
+------------------------------------------------------------------------- */
+
+double FixPrecessionSpin::compute_cubic_energy(double spi[3])
+{
+  double energy = 0.0;
+  double skx,sky,skz;
+
+  skx = spi[0]*nc1x+spi[1]*nc1y+spi[2]*nc1z;
+  sky = spi[0]*nc2x+spi[1]*nc2y+spi[2]*nc2z;
+  skz = spi[0]*nc3x+spi[1]*nc3y+spi[2]*nc3z;
+
+  energy = k1c*(skx*skx*sky*sky + sky*sky*skz*skz + skx*skx*skz*skz);
+  energy += k2c*skx*skx*sky*sky*skz*skz;
+
+  return 2.0*energy;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPrecessionSpin::set_magneticprecession()
+{
+  if (zeeman_flag) {
+    hx = H_field*nhx;
+    hy = H_field*nhy;
+    hz = H_field*nhz;
+  }
+  if (aniso_flag) {
+    Kax = 2.0*Kah*nax;
+    Kay = 2.0*Kah*nay;
+    Kaz = 2.0*Kah*naz;
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -422,3 +417,10 @@ void FixPrecessionSpin::min_post_force(int vflag)
 {
   post_force(vflag);
 }
+
+/* ---------------------------------------------------------------------- */
+
+void FixPrecessionSpin::post_force_respa(int vflag, int ilevel, int /*iloop*/)
+{
+  if (ilevel == ilevel_respa) post_force(vflag);
+}