diff --git a/doc/src/Commands_all.txt b/doc/src/Commands_all.txt
index 52c2e67e2e..58ca148555 100644
--- a/doc/src/Commands_all.txt
+++ b/doc/src/Commands_all.txt
@@ -83,7 +83,7 @@ An alphabetic list of all general LAMMPS commands.
 "molecule"_molecule.html,
 "ndx2group"_group2ndx.html,
 "neb"_neb.html,
-"neb_spin"_neb_spin.html,
+"neb/spin"_neb_spin.html,
 "neigh_modify"_neigh_modify.html,
 "neighbor"_neighbor.html,
 "newton"_newton.html,
diff --git a/doc/src/Eqs/fix_spin_cubic.jpg b/doc/src/Eqs/fix_spin_cubic.jpg
new file mode 100644
index 0000000000..871641111c
Binary files /dev/null and b/doc/src/Eqs/fix_spin_cubic.jpg differ
diff --git a/doc/src/Eqs/fix_spin_cubic.tex b/doc/src/Eqs/fix_spin_cubic.tex
new file mode 100644
index 0000000000..9c402294d7
--- /dev/null
+++ b/doc/src/Eqs/fix_spin_cubic.tex
@@ -0,0 +1,21 @@
+\documentclass[preview]{standalone}
+\usepackage{varwidth}
+\usepackage[utf8x]{inputenc}
+\usepackage{amsmath,amssymb,amsthm,bm}
+\begin{document}
+\begin{varwidth}{50in}
+  \begin{equation}
+    \bm{H}_{cubic} = -\sum_{{ i}=1}^{N} K_{1}
+    \Big[
+    \left(\vec{s}_{i} \cdot \vec{n1} \right)^2
+    \left(\vec{s}_{i} \cdot \vec{n2} \right)^2 +
+    \left(\vec{s}_{i} \cdot \vec{n2} \right)^2
+    \left(\vec{s}_{i} \cdot \vec{n3} \right)^2 +
+    \left(\vec{s}_{i} \cdot \vec{n1} \right)^2
+    \left(\vec{s}_{i} \cdot \vec{n3} \right)^2 \Big]
+    +K_{2}^{(c)} \left(\vec{s}_{i} \cdot \vec{n1} \right)^2
+    \left(\vec{s}_{i} \cdot \vec{n2} \right)^2
+    \left(\vec{s}_{i} \cdot \vec{n3} \right)^2 \nonumber
+  \end{equation}
+\end{varwidth}
+\end{document}
diff --git a/doc/src/Tools.txt b/doc/src/Tools.txt
index af9fd4298c..eb7b6d81b8 100644
--- a/doc/src/Tools.txt
+++ b/doc/src/Tools.txt
@@ -77,6 +77,7 @@ Post-processing tools :h3
 "python"_#pythontools,
 "reax"_#reax_tool,
 "smd"_#smd,
+"spin"_#spin,
 "xmgrace"_#xmgrace :tb(c=6,ea=c,a=l)
 
 Miscellaneous tools :h3
@@ -511,6 +512,20 @@ Ernst Mach Institute in Germany (georg.ganzenmueller at emi.fhg.de).
 
 :line
 
+spin tool :h4,link(spin)
+
+The spin sub-directory contains a C file interpolate.c which can
+be compiled and used to perform a cubic polynomial interpolation of 
+the MEP following a GNEB calculation.
+
+See the README file in tools/spin/interpolate_gneb for more details.
+
+This tool was written by the SPIN package author, Julien
+Tranchida at Sandia National Labs (jtranch at sandia.gov, and by Aleksei 
+Ivanov, at University of Iceland (ali5 at hi.is).
+
+:line
+
 vim tool :h4,link(vim)
 
 The files in the tools/vim directory are add-ons to the VIM editor
diff --git a/doc/src/commands_list.txt b/doc/src/commands_list.txt
index cf716df9ac..e224dc8096 100644
--- a/doc/src/commands_list.txt
+++ b/doc/src/commands_list.txt
@@ -62,12 +62,12 @@ Commands :h1
    mass
    message
    min_modify
-   min_spin
+   min/spin
    min_style
    minimize
    molecule
    neb
-   neb_spin
+   neb/spin
    neigh_modify
    neighbor
    newton
diff --git a/doc/src/fix_precession_spin.txt b/doc/src/fix_precession_spin.txt
index 05814931ea..708b2bd7aa 100644
--- a/doc/src/fix_precession_spin.txt
+++ b/doc/src/fix_precession_spin.txt
@@ -14,19 +14,23 @@ fix ID group precession/spin style args :pre
 
 ID, group are documented in "fix"_fix.html command :ulb,l
 precession/spin = style name of this fix command :l
-style = {zeeman} or {anisotropy} :l
+style = {zeeman} or {anisotropy} or {cubic} :l
   {zeeman} args = H x y z
     H = intensity of the magnetic field (in Tesla)
     x y z = vector direction of the field
   {anisotropy} args = K x y z
     K = intensity of the magnetic anisotropy (in eV)
     x y z = vector direction of the anisotropy :pre
+  {cubic} args = K1 K2c n1x n1y n1x n2x n2y n2z n3x n3y n3z 
+    K1 and K2c = intensity of the magnetic anisotropy (in eV)
+    n1x to n3z = three direction vectors of the cubic anisotropy :pre
 :ule
 
 [Examples:]
 
 fix 1 all precession/spin zeeman 0.1 0.0 0.0 1.0
-fix 1 all precession/spin anisotropy 0.001 0.0 0.0 1.0
+fix 1 3 precession/spin anisotropy 0.001 0.0 0.0 1.0
+fix 1 iron precession/spin cubic 0.001 0.0005 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0
 fix 1 all precession/spin zeeman 0.1 0.0 0.0 1.0 anisotropy 0.001 0.0 0.0 1.0 :pre
 
 [Description:]
@@ -50,10 +54,29 @@ for the magnetic spins in the defined group:
 with n defining the direction of the anisotropy, and K (in eV) its intensity.
 If K>0, an easy axis is defined, and if K<0, an easy plane is defined.
 
-In both cases, the choice of (x y z) imposes the vector direction for the force.
-Only the direction of the vector is important; it's length is ignored.
+Style {cubic} is used to simulate a cubic anisotropy, with three
+possible easy axis for the magnetic spins in the defined group: 
 
-Both styles can be combined within one single command line.
+:c,image(Eqs/fix_spin_cubic.jpg)
+
+with K1 and K2c (in eV) the intensity coefficients and 
+n1, n2 and n3 defining the three anisotropic directions
+defined by the command (from n1x to n3z). 
+For n1 = (100), n2 = (010), and n3 = (001), K1 < 0 defines an 
+iron type anisotropy (easy axis along the (001)-type cube
+edges), and K1 > 0 defines a nickel type anisotropy (easy axis
+along the (111)-type cube diagonals). 
+K2^c > 0 also defines easy axis along the (111)-type cube
+diagonals.
+See chapter 2 of "(Skomski)"_#Skomski1 for more details on cubic
+anisotropies.
+
+In all cases, the choice of (x y z) only imposes the vector
+directions for the forces. Only the direction of the vector is 
+important; it's length is ignored (the entered vectors are
+normalized).
+
+Those styles can be combined within one single command line.
 
 :line
 
@@ -85,3 +108,9 @@ package"_Build_package.html doc page for more info.
 "atom_style spin"_atom_style.html
 
 [Default:] none
+
+:line
+
+:link(Skomski1)
+[(Skomski)] Skomski, R. (2008). Simple models of magnetism.
+Oxford University Press.
diff --git a/doc/src/neb_spin.txt b/doc/src/neb_spin.txt
index 70c0fe2ab5..7dbd924cd2 100644
--- a/doc/src/neb_spin.txt
+++ b/doc/src/neb_spin.txt
@@ -6,7 +6,7 @@
 
 :line
 
-neb command :h3
+neb/spin command :h3
 
 [Syntax:]
 
diff --git a/examples/SPIN/iron/in.spin.iron_cubic b/examples/SPIN/iron/in.spin.iron_cubic
new file mode 100644
index 0000000000..d4703a2959
--- /dev/null
+++ b/examples/SPIN/iron/in.spin.iron_cubic
@@ -0,0 +1,60 @@
+# bcc iron in a 3d periodic box
+
+clear 
+units 		metal
+atom_style 	spin
+
+dimension 	3
+boundary 	p p p
+
+# necessary for the serial algorithm (sametag)
+atom_modify 	map array 
+
+lattice 	bcc 2.8665
+region 		box block 0.0 5.0 0.0 5.0 0.0 5.0
+create_box 	1 box
+create_atoms 	1 box
+
+# setting mass, mag. moments, and interactions for bcc iron
+
+mass		1 55.845
+set 		group all spin 2.2 -1.0 0.0 0.0
+velocity 	all create 100 4928459 rot yes dist gaussian
+
+pair_style 	hybrid/overlay eam/alloy spin/exchange 3.5
+pair_coeff 	* * eam/alloy Fe_Mishin2006.eam.alloy Fe
+pair_coeff 	* * spin/exchange exchange 3.4 0.02726 0.2171 1.841
+neighbor 	0.1 bin
+neigh_modify 	every 10 check yes delay 20
+
+fix 		1 all precession/spin cubic 0.001 0.0005 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0
+fix_modify 	1 energy yes
+fix 		2 all langevin/spin 0.0 0.0 21
+
+fix 		3 all nve/spin lattice yes
+timestep	0.0001
+
+# compute and output options
+
+compute 	out_mag    all spin
+compute 	out_pe     all pe
+compute 	out_ke     all ke
+compute 	out_temp   all temp
+
+variable 	magx      equal c_out_mag[1]
+variable 	magy      equal c_out_mag[2]
+variable 	magz      equal c_out_mag[3]
+variable 	magnorm   equal c_out_mag[4]
+variable 	emag      equal c_out_mag[5]
+variable 	tmag      equal c_out_mag[6]
+
+thermo_style    custom step time v_magx v_magy v_magz v_magnorm v_tmag v_emag pe etotal
+thermo          50
+
+compute 	outsp all property/atom spx spy spz sp fmx fmy fmz
+dump 		100 all custom 1 dump_iron.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3]
+
+run 		2000
+# min_style 	spin
+# min_modify	alpha_damp 1.0 discrete_factor 10
+# minimize 	1.0e-16 1.0e-16 10000 10000
diff --git a/examples/SPIN/nickel/in.spin.nickel_cubic b/examples/SPIN/nickel/in.spin.nickel_cubic
new file mode 100644
index 0000000000..3c97b284ae
--- /dev/null
+++ b/examples/SPIN/nickel/in.spin.nickel_cubic
@@ -0,0 +1,60 @@
+# fcc nickel in a 3d periodic box
+
+clear 
+units 		metal
+atom_style 	spin
+
+dimension 	3
+boundary 	p p p
+
+# necessary for the serial algorithm (sametag)
+atom_modify 	map array 
+
+lattice 	fcc 3.524
+region 		box block 0.0 5.0 0.0 5.0 0.0 5.0
+create_box 	1 box
+create_atoms 	1 box
+
+# setting mass, mag. moments, and interactions for cobalt
+
+mass		1 58.69
+
+set 		group all spin/random 31 0.63
+#set 		group all spin 0.63 0.0 0.0 1.0 
+velocity 	all create 100 4928459 rot yes dist gaussian
+
+pair_style 	hybrid/overlay eam/alloy spin/exchange 4.0
+pair_coeff 	* * eam/alloy Ni99.eam.alloy Ni
+pair_coeff 	* * spin/exchange exchange 4.0 0.50 0.2280246862 1.229983475
+
+neighbor 	0.1 bin
+neigh_modify 	every 10 check yes delay 20
+
+fix 		1 all precession/spin cubic -0.0001 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 &
+				      zeeman 0.0 0.0 0.0 1.0
+fix_modify 	1 energy yes
+fix 		2 all langevin/spin 0.0 0.0 21
+
+fix 		3 all nve/spin lattice yes
+timestep	0.0001
+
+# compute and output options
+
+compute 	out_mag    all spin
+compute 	out_pe     all pe
+compute 	out_ke     all ke
+compute 	out_temp   all temp
+
+variable 	magz      equal c_out_mag[3]
+variable 	magnorm   equal c_out_mag[4]
+variable 	emag      equal c_out_mag[5]
+variable 	tmag      equal c_out_mag[6]
+
+thermo_style    custom step time v_magnorm v_emag temp v_tmag etotal
+thermo          50
+
+compute 	outsp all property/atom spx spy spz sp fmx fmy fmz
+dump 		50 all custom 1 dump.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3] c_outsp[4] c_outsp[5] c_outsp[6] c_outsp[7]
+
+run 		2000
+
diff --git a/src/.gitignore b/src/.gitignore
index 9e10059ddb..764bc8c2f3 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -156,6 +156,8 @@
 /fix_nve_spin.h
 /fix_precession_spin.cpp
 /fix_precession_spin.h
+/fix_setforce_spin.cpp
+/fix_setforce_spin.h
 /min_spin.cpp
 /min_spin.h
 /neb_spin.cpp
diff --git a/src/SPIN/fix_precession_spin.cpp b/src/SPIN/fix_precession_spin.cpp
index 9ee691c227..1ea134d367 100644
--- a/src/SPIN/fix_precession_spin.cpp
+++ b/src/SPIN/fix_precession_spin.cpp
@@ -25,8 +25,8 @@
 #include <cstdio>
 #include <cstdlib>
 #include <cstring>
-
 #include "atom.h"
+#include "error.h"
 #include "domain.h"
 #include "error.h"
 #include "fix_precession_spin.h"
@@ -35,6 +35,7 @@
 #include "math_const.h"
 #include "memory.h"
 #include "modify.h"
+#include "neigh_list.h"
 #include "respa.h"
 #include "update.h"
 #include "variable.h"
@@ -71,8 +72,12 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lm
   Ka = 0.0;
   nax = nay = naz = 0.0;
   Kax = Kay = Kaz = 0.0;
+  k1c = k2c = 0.0;
+  nc1x = nc1y = nc1z = 0.0;
+  nc2x = nc2y = nc2z = 0.0;
+  nc3x = nc3y = nc3z = 0.0;
 
-  zeeman_flag = aniso_flag = 0;
+  zeeman_flag = aniso_flag = cubic_flag = 0;
 
   int iarg = 3;
   while (iarg < narg) {
@@ -92,14 +97,61 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lm
       nay = force->numeric(FLERR,arg[iarg+3]);
       naz = force->numeric(FLERR,arg[iarg+4]);
       iarg += 5;
+    } else if (strcmp(arg[iarg],"cubic") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix precession/spin command");
+      cubic_flag = 1;
+      k1c = force->numeric(FLERR,arg[iarg+1]);
+      k2c = force->numeric(FLERR,arg[iarg+2]);
+      nc1x = force->numeric(FLERR,arg[iarg+3]);
+      nc1y = force->numeric(FLERR,arg[iarg+4]);
+      nc1z = force->numeric(FLERR,arg[iarg+5]);
+      nc2x = force->numeric(FLERR,arg[iarg+6]);
+      nc2y = force->numeric(FLERR,arg[iarg+7]);
+      nc2z = force->numeric(FLERR,arg[iarg+8]);
+      nc3x = force->numeric(FLERR,arg[iarg+9]);
+      nc3y = force->numeric(FLERR,arg[iarg+10]);
+      nc3z = force->numeric(FLERR,arg[iarg+11]);
+      iarg += 12;
     } else error->all(FLERR,"Illegal precession/spin command");
   }
 
+  // normalize vectors
+
+  double inorm;
+  if (zeeman_flag) {
+    inorm = 1.0/sqrt(nhx*nhx + nhy*nhy + nhz*nhz);
+    nhx *= inorm;
+    nhy *= inorm;
+    nhz *= inorm;
+  }
+
+  if (aniso_flag) {
+    inorm = 1.0/sqrt(nax*nax + nay*nay + naz*naz);
+    nax *= inorm;
+    nay *= inorm;
+    naz *= inorm;
+  }
+  
+  if (cubic_flag) {
+    inorm = 1.0/sqrt(nc1x*nc1x + nc1y*nc1y + nc1z*nc1z);
+    nc1x *= inorm;
+    nc1y *= inorm;
+    nc1z *= inorm;
+    inorm = 1.0/sqrt(nc2x*nc2x + nc2y*nc2y + nc2z*nc2z);
+    nc2x *= inorm;
+    nc2y *= inorm;
+    nc2z *= inorm;
+    inorm = 1.0/sqrt(nc3x*nc3x + nc3y*nc3y + nc3z*nc3z);
+    nc3x *= inorm;
+    nc3y *= inorm;
+    nc3z *= inorm;
+  }
+
   degree2rad = MY_PI/180.0;
   time_origin = update->ntimestep;
 
   eflag = 0;
-  emag = 0.0;
+  eprec = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -130,8 +182,12 @@ void FixPrecessionSpin::init()
   const double mub = 5.78901e-5;                // in eV/T
   const double gyro = mub/hbar;                 // in rad.THz/T
 
-  H_field *= gyro;                              // in rad.THz
-  Ka /= hbar;                                   // in rad.THz
+  // convert field quantities to rad.THz
+
+  H_field *= gyro;
+  Kah = Ka/hbar;
+  k1ch = k1c/hbar;
+  k2ch = k2c/hbar;
 
   if (strstr(update->integrate_style,"respa")) {
     ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
@@ -185,53 +241,60 @@ void FixPrecessionSpin::post_force(int /* vflag */)
   if (varflag != CONSTANT) {
     modify->clearstep_compute();
     modify->addstep_compute(update->ntimestep + 1);
-    set_magneticprecession();                   // update mag. field if time-dep.
+    set_magneticprecession();		// update mag. field if time-dep.
   }
 
-  double **sp = atom->sp;
+  int *mask = atom->mask;
   double **fm = atom->fm;
-  double spi[3], fmi[3];
+  double **sp = atom->sp;
   const int nlocal = atom->nlocal;
-
+  double spi[3], fmi[3], epreci;
+  
   eflag = 0;
-  emag = 0.0;
-
+  eprec = 0.0;
   for (int i = 0; i < nlocal; i++) {
-    spi[0] = sp[i][0];
-    spi[1] = sp[i][1];
-    spi[2] = sp[i][2];
-    fmi[0] = fmi[1] = fmi[2] = 0.0;
+    if (mask[i] & groupbit) {
+      epreci = 0.0;
+      spi[0] = sp[i][0];
+      spi[1] = sp[i][1];
+      spi[2] = sp[i][2];
+      fmi[0] = fmi[1] = fmi[2] = 0.0;
 
-    if (zeeman_flag) {          // compute Zeeman interaction
-      compute_zeeman(i,fmi);
-      emag -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
+      if (zeeman_flag) {          // compute Zeeman interaction
+        compute_zeeman(i,fmi);
+        epreci -= hbar*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
+      }
+
+      if (aniso_flag) {           // compute magnetic anisotropy
+        compute_anisotropy(spi,fmi);
+        epreci -= compute_anisotropy_energy(spi);
+      }
+
+      if (cubic_flag) {		// compute cubic anisotropy
+	compute_cubic(spi,fmi);
+	epreci -= compute_cubic_energy(spi);
+      }
+
+      eprec += epreci;
+      fm[i][0] += fmi[0];
+      fm[i][1] += fmi[1];
+      fm[i][2] += fmi[2];
     }
-
-    if (aniso_flag) {           // compute magnetic anisotropy
-      compute_anisotropy(spi,fmi);
-      emag -= 0.5*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
-    }
-
-    fm[i][0] += fmi[0];
-    fm[i][1] += fmi[1];
-    fm[i][2] += fmi[2];
   }
-  emag *= hbar;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixPrecessionSpin::compute_single_precession(int i, double spi[3], double fmi[3])
 {
-  if (zeeman_flag) {
-    compute_zeeman(i,fmi);
-  }
-  if (aniso_flag) {
-    compute_anisotropy(spi,fmi);
+  int *mask = atom->mask;
+  if (mask[i] & groupbit) {
+    if (zeeman_flag) compute_zeeman(i,fmi);
+    if (aniso_flag) compute_anisotropy(spi,fmi);
+    if (cubic_flag) compute_cubic(spi,fmi);
   }
 }
 
-
 /* ---------------------------------------------------------------------- */
 
 void FixPrecessionSpin::compute_zeeman(int i, double fmi[3])
@@ -254,6 +317,16 @@ void FixPrecessionSpin::compute_anisotropy(double spi[3], double fmi[3])
 
 /* ---------------------------------------------------------------------- */
 
+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;
+  return energy; 
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixPrecessionSpin::post_force_respa(int vflag, int ilevel, int /*iloop*/)
 {
   if (ilevel == ilevel_respa) post_force(vflag);
@@ -264,17 +337,75 @@ void FixPrecessionSpin::post_force_respa(int vflag, int ilevel, int /*iloop*/)
 void FixPrecessionSpin::set_magneticprecession()
 {
   if (zeeman_flag) {
-          hx = H_field*nhx;
-          hy = H_field*nhy;
-          hz = H_field*nhz;
+    hx = H_field*nhx;
+    hy = H_field*nhy;
+    hz = H_field*nhz;
   }
   if (aniso_flag) {
-          Kax = 2.0*Ka*nax;
-          Kay = 2.0*Ka*nay;
-          Kaz = 2.0*Ka*naz;
+    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
+------------------------------------------------------------------------- */
+
+void FixPrecessionSpin::compute_cubic(double spi[3], double fmi[3])
+{
+  double skx,sky,skz,skx2,sky2,skz2;
+  double four1,four2,four3,fourx,foury,fourz;
+  double six1,six2,six3,sixx,sixy,sixz;
+
+  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;
+
+  skx2 = skx*skx;
+  sky2 = sky*sky;
+  skz2 = skz*skz;
+
+  four1 = 2.0*skx*(sky2+skz2);
+  four2 = 2.0*sky*(skx2+skz2);
+  four3 = 2.0*skz*(skx2+sky2);
+
+  fourx = k1ch*(nc1x*four1 + nc2x*four2 + nc3x*four3);
+  foury = k1ch*(nc1y*four1 + nc2y*four2 + nc3y*four3);
+  fourz = k1ch*(nc1z*four1 + nc2z*four2 + nc3z*four3);
+
+  six1 = 2.0*skx*sky2*skz2;
+  six2 = 2.0*sky*skx2*skz2;
+  six3 = 2.0*skz*skx2*sky2;
+  
+  sixx = k2ch*(nc1x*six1 + nc2x*six2 + nc3x*six3);
+  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;
+}
+
 /* ----------------------------------------------------------------------
    potential energy in magnetic field
 ------------------------------------------------------------------------- */
@@ -284,10 +415,10 @@ double FixPrecessionSpin::compute_scalar()
   // only sum across procs one time
 
   if (eflag == 0) {
-    MPI_Allreduce(&emag,&emag_all,1,MPI_DOUBLE,MPI_SUM,world);
+    MPI_Allreduce(&eprec,&eprec_all,1,MPI_DOUBLE,MPI_SUM,world);
     eflag = 1;
   }
-  return emag_all;
+  return eprec_all;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/SPIN/fix_precession_spin.h b/src/SPIN/fix_precession_spin.h
index 1db4d32ae9..0037784a48 100644
--- a/src/SPIN/fix_precession_spin.h
+++ b/src/SPIN/fix_precession_spin.h
@@ -39,11 +39,20 @@ class FixPrecessionSpin : public Fix {
   void min_post_force(int);
   double compute_scalar();
 
-  int zeeman_flag, aniso_flag;
+  int zeeman_flag, aniso_flag, cubic_flag;
   void compute_single_precession(int, double *, double *);
   void compute_zeeman(int, double *);
-  void compute_anisotropy(double *, double *);
+  
+  // uniaxial aniso calculations
 
+  void compute_anisotropy(double *, double *);
+  double compute_anisotropy_energy(double *);
+
+  // cubic aniso calculations
+
+  void compute_cubic(double *, double *);
+  double compute_cubic_energy(double *);
+ 
  protected:
   int style; 			// style of the magnetic precession
 
@@ -52,7 +61,7 @@ class FixPrecessionSpin : public Fix {
   int ilevel_respa;
   int time_origin;
   int eflag;
-  double emag, emag_all;
+  double eprec, eprec_all;
 
   int varflag;
   int magfieldstyle;
@@ -67,10 +76,19 @@ class FixPrecessionSpin : public Fix {
 
   // magnetic anisotropy intensity and direction
 
-  double Ka;
+  double Ka;			// aniso const. in eV
+  double Kah;			// aniso const. in rad.THz
   double nax, nay, naz;
   double Kax, Kay, Kaz; 	// temp. force variables
 
+  // cubic anisotropy intensity
+
+  double k1c,k2c;		// cubic const. in eV
+  double k1ch,k2ch;		// cubic const. in rad.THz
+  double nc1x,nc1y,nc1z;
+  double nc2x,nc2y,nc2z;
+  double nc3x,nc3y,nc3z;
+
   void set_magneticprecession();
 
 };
diff --git a/tools/README b/tools/README
index 6bf7d6b878..54f8d86898 100644
--- a/tools/README
+++ b/tools/README
@@ -39,6 +39,7 @@ pymol_asphere	       convert LAMMPS output of ellipsoids to PyMol format
 python		       Python scripts for post-processing LAMMPS output
 reax	       	       Tools for analyzing output of ReaxFF simulations
 smd                    convert Smooth Mach Dynamics triangles to VTK
+spin                   perform a cubic polynomial interpolation of a GNEB MEP
 vim		       add-ons to VIM editor for editing LAMMPS input scripts
 xmgrace                a collection of scripts to generate xmgrace plots