diff --git a/in.co_magnetomech b/in.co_magnetomech
index da9e5366f1..40a6d37a2c 100644
--- a/in.co_magnetomech
+++ b/in.co_magnetomech
@@ -63,6 +63,7 @@ pair_coeff * * eam/alloy ../Co_PurjaPun_2012.eam.alloy Co
 #pair_style pair/spin 4.0
 #type i and j | interaction type | cutoff | J1 (eV) | J2 (adim) | J3 (Ang) (for Exchange)
 pair_coeff * * pair/spin exchange 4.0 0.0446928 0.003496 1.4885 
+#pair_coeff * * pair/spin exchange 4.0 0.0 0.003496 1.4885 
 #type i and j | interaction type | cutoff | Int (eV) | [dx,dy,dz] (for DMI and ME)
 #pair_coeff * * dmi 2.6 0.001 0.0 0.0 1.0
 #pair_coeff * * me 2.6 0.01 1.0 1.0 1.0
@@ -76,8 +77,8 @@ neigh_modify every 1 check no delay 0
 
 #Magnetic field fix
 #Type | Intensity (T or eV) | Direction
-fix 1 all force/spin zeeman 10.0 0.0 0.0 1.0
-#fix 1 all force/spin anisotropy 0.001 0.0 0.0 1.0
+fix 1 all force/spin zeeman 0.1 0.0 0.0 1.0
+#fix 1 all force/spin anisotropy 0.01 0.0 0.0 1.0
 #fix 1 all force/spin anisotropy -0.1 0.0 0.0 1.0
 #fix 1 all force/spin anisotropy 0.1 0.0 1.0 0.0
 
@@ -85,7 +86,7 @@ fix 1 all force/spin zeeman 10.0 0.0 0.0 1.0
 #Temp | Alpha_trans | Alpha_long | Seed
 #fix 2 all langevin/spin 0.0 0.1 0.0 21
 #fix 2 all langevin/spin 1.0 0.1 0.0 21
-fix 2 all langevin/spin 0.0 0.0 0.0 21
+fix 2 all langevin/spin 0.1 0.01 0.0 21
 
 #Magnetic integration fix
 fix 3 all integration/spin mpi
@@ -95,6 +96,14 @@ fix 3 all integration/spin mpi
 compute mag all compute/spin
 fix outmag all ave/time 1 1 10 c_mag[1] c_mag[2] c_mag[3] c_mag[4] c_mag[5] c_mag[6] c_mag[7] file mag_Co_nodamp.dat format %20.16g
 
+#compute mechanical energy
+compute mech all pe
+fix outmech all ave/time 1 1 10 c_mech file mech_Co_nodamp.dat format %20.16g
+
+#compute kinetic energy
+compute kinetic all ke
+fix outke all ave/time 1 1 10 c_kinetic file kinetic_Co_nodamp.dat format %20.16g
+
 #Setting the timestep for the simulation (in ps)
 timestep	0.0001
 
@@ -106,6 +115,6 @@ timestep	0.0001
 dump 1 all custom 100 dump_spin_MM.lammpstrj type x y z spx spy spz
 
 #Running the simulations for N timesteps
-run 100000
+run 50000
 #run 1
 
diff --git a/in.cobalt b/in.cobalt
index ebc2017a5d..9fff327a51 100644
--- a/in.cobalt
+++ b/in.cobalt
@@ -59,7 +59,8 @@ set group all spin 1.72 1.0 0.0 0.0
 #Magnetic exchange interaction coefficient for bulk fcc Cobalt
 pair_style pair/spin 4.0
 #type i and j | interaction type | cutoff | J1 (eV) | J2 (adim) | J3 (Ang) (for Exchange)
-pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885 
+#pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885 
+pair_coeff * * exchange 4.0 0.0 0.003496 1.4885 
 #type i and j | interaction type | cutoff | Int (eV) | [dx,dy,dz] (for DMI and ME)
 #pair_coeff * * dmi 2.6 0.001 0.0 0.0 1.0
 #pair_coeff * * me 2.6 0.01 1.0 1.0 1.0
@@ -72,19 +73,19 @@ neigh_modify every 1 check no delay 0
 
 #Magnetic field fix
 #Type | Intensity (T or eV) | Direction
-fix 1 all force/spin zeeman 1.0 0.0 0.0 1.0
-#fix 1 all force/spin anisotropy 0.001 0.0 0.0 1.0
+#fix 1 all force/spin zeeman 1.0 0.0 0.0 1.0
+fix 1 all force/spin anisotropy 0.00 0.0 0.0 1.0
 #fix 1 all force/spin anisotropy -0.1 0.0 0.0 1.0
 #fix 1 all force/spin anisotropy 0.1 0.0 1.0 0.0
 
 #Fix Langevin spins (merging damping and temperature)
 #Temp | Alpha_trans | Alpha_long | Seed
 #fix 2 all langevin/spin 0.0 0.1 0.0 21
-fix 2 all langevin/spin 1.0 0.1 0.0 21
+fix 2 all langevin/spin 0.0 0.0 0.0 21
 #fix 2 all langevin/spin 0.0 0.0 0.0 21
 
 #Magnetic integration fix
-fix 3 all nve/spin mpi
+fix 3 all integration/spin mpi
 
 #compute real time, total magnetization, magnetic energy, and spin temperature
 #Iteration | Time | Mx | My | Mz | |M| | Em | Tm
@@ -102,6 +103,6 @@ timestep	0.0001
 dump 1 all custom 100 dump_spin_T100.lammpstrj type x y z spx spy spz
 
 #Running the simulations for N timesteps
-run 20
-#run 1
+#run 10000
+run 1
 
diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp
index 6018fc28f0..ccdaa13829 100644
--- a/src/SPIN/atom_vec_spin.cpp
+++ b/src/SPIN/atom_vec_spin.cpp
@@ -9,6 +9,12 @@
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
+
+------------------------------------------------------------------------- */
+
+/* ------------------------------------------------------------------------
+   Contributing authors: Julien Tranchida (SNL)
+                         Aidan Thompson (SNL)
 ------------------------------------------------------------------------- */
 
 #include <math.h>
@@ -30,16 +36,16 @@ using namespace LAMMPS_NS;
 AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp)
 {
   molecular = 0;
-  mass_type = 1; //check why
+  mass_type = 1; // check why
 
   comm_x_only = 0;
-  comm_f_only = 1;
+  comm_f_only = 0;
 
   size_forward = 7;
   size_reverse = 6;
   size_border = 11;
   size_velocity = 3;
-  size_data_atom = 9; //to check later
+  size_data_atom = 10; // to check later
   size_data_vel = 4;
   xcol_data = 4;
  
@@ -67,10 +73,11 @@ void AtomVecSpin::grow(int n)
   type = memory->grow(atom->type,nmax,"atom:type");
   mask = memory->grow(atom->mask,nmax,"atom:mask");
   image = memory->grow(atom->image,nmax,"atom:image");
+  // allocating mech. quantities
   x = memory->grow(atom->x,nmax,3,"atom:x");
   v = memory->grow(atom->v,nmax,3,"atom:v");
   f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
-  //Allocating mag. quantities
+  // allocating mag. quantities
   mumag = memory->grow(atom->mumag,nmax,"atom:mumag");
   sp = memory->grow(atom->sp,nmax,4,"atom:sp");
   fm = memory->grow(atom->fm,nmax*comm->nthreads,3,"atom:fm");
@@ -330,6 +337,7 @@ int AtomVecSpin::pack_reverse(int n, int first, double *buf)
     buf[m++] = fm[i][1];
     buf[m++] = fm[i][2];
   }
+
   return m;
 }
 
@@ -934,9 +942,9 @@ bigint AtomVecSpin::memory_usage()
   if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
   if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
   
-  if (atom->memcheck("fm")) bytes += memory->usage(fm,nmax*comm->nthreads,3);
   if (atom->memcheck("mumag")) bytes += memory->usage(mumag,nmax);
   if (atom->memcheck("sp")) bytes += memory->usage(sp,nmax,4);
+  if (atom->memcheck("fm")) bytes += memory->usage(fm,nmax*comm->nthreads,3);
 
   return bytes;
 }
diff --git a/src/SPIN/atom_vec_spin.h b/src/SPIN/atom_vec_spin.h
index f4b13926f0..aedc3efb3e 100644
--- a/src/SPIN/atom_vec_spin.h
+++ b/src/SPIN/atom_vec_spin.h
@@ -58,7 +58,7 @@ class AtomVecSpin : public AtomVec {
   int write_data_hybrid(FILE *, double *); 
   bigint memory_usage();
 
-  //Test force clear
+  // clear mag. and mech. forces
   void force_clear(int, size_t); 
 
 
@@ -66,8 +66,8 @@ class AtomVecSpin : public AtomVec {
   tagint *tag;
   int *type,*mask;
   imageint *image;
-  double **x,**v,**f;     //MD quantities: position, velocity and force
-  double *mumag,**sp, **fm;         //Magnetic quantities: mu, spin direction, magnetic force
+  double **x,**v,**f;       // lattice quantities
+  double *mumag,**sp, **fm; // spin quantities
 
 };
 
diff --git a/src/SPIN/compute_spin.cpp b/src/SPIN/compute_spin.cpp
index ed69a76229..140d6ce8ce 100644
--- a/src/SPIN/compute_spin.cpp
+++ b/src/SPIN/compute_spin.cpp
@@ -11,6 +11,11 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+/* ------------------------------------------------------------------------
+   Contributing authors: Julien Tranchida (SNL)
+                         Aidan Thompson (SNL)
+------------------------------------------------------------------------- */
+
 #include <mpi.h>
 #include <string.h>
 #include "atom.h"
@@ -86,9 +91,9 @@ void ComputeSpin::compute_vector()
   double tx,ty,tz;
 	
   int nlocal = atom->nlocal;
-    
+
   // compute total magnetization and magnetic energy 
-  // compute spin temperature; See Nurdin et al., Phys. Rev. E 61, 2000
+  // compute spin temperature (Nurdin et al., Phys. Rev. E 61, 2000)
   for (i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
       if (atom->mumag_flag && atom->sp_flag) {
@@ -108,7 +113,7 @@ void ComputeSpin::compute_vector()
       }
       else error->all(FLERR,"Compute spin/compute declared magnetic quantities (sp and mumag flags)");
   }
-  
+ 
   MPI_Allreduce(mag,magtot,4,MPI_DOUBLE,MPI_SUM,world);
   MPI_Allreduce(&magenergy,&magenergytot,1,MPI_DOUBLE,MPI_SUM,world);
   MPI_Allreduce(&tempnum,&tempnumtot,1,MPI_DOUBLE,MPI_SUM,world);
diff --git a/src/SPIN/fix_force_spin.cpp b/src/SPIN/fix_force_spin.cpp
index daa34e9df6..9791d622a0 100644
--- a/src/SPIN/fix_force_spin.cpp
+++ b/src/SPIN/fix_force_spin.cpp
@@ -11,6 +11,11 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+/* ------------------------------------------------------------------------
+   Contributing authors: Julien Tranchida (SNL)
+                         Aidan Thompson (SNL)
+------------------------------------------------------------------------- */
+
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -41,7 +46,7 @@ enum{CONSTANT,EQUAL};
 FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
 {
 	
-  if (narg < 7) error->all(FLERR,"Illegal fix spin command");
+  if (narg < 7) error->all(FLERR,"Illegal force/spin command");
   // 7 arguments for a force/spin fix command:
   //(fix  ID  group  force/spin  magnitude (T or eV)  style (zeeman or anisotropy)  direction (3 cartesian coordinates) 
   
@@ -65,7 +70,7 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a
   zeeman_flag = aniso_flag = 0;
  
   if (strcmp(arg[3],"zeeman") == 0) {
-	  if (narg != 8) error->all(FLERR,"Illegal fix zeeman command");
+	  if (narg != 8) error->all(FLERR,"Illegal force/spin command");
 	  style = ZEEMAN;
           zeeman_flag = 1;
 	  H_field = force->numeric(FLERR,arg[4]);
@@ -74,14 +79,14 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a
 	  Hz = force->numeric(FLERR,arg[7]);	
 	  magfieldstyle = CONSTANT; 	  	  
   } else if (strcmp(arg[3],"anisotropy") == 0) {
-	  if (narg != 8) error->all(FLERR,"Illegal fix anisotropy command");
+	  if (narg != 8) error->all(FLERR,"Illegal force/spin command");
 	  style = ANISOTROPY;
           aniso_flag = 1;
 	  Ka = force->numeric(FLERR,arg[4]);
 	  Kax = force->numeric(FLERR,arg[5]);
 	  Kay = force->numeric(FLERR,arg[6]);
 	  Kaz = force->numeric(FLERR,arg[7]);	  
-  } else error->all(FLERR,"Illegal fix force/spin command");
+  } else error->all(FLERR,"Illegal force/spin command");
  
   degree2rad = MY_PI/180.0;
   time_origin = update->ntimestep;
@@ -115,12 +120,12 @@ int FixForceSpin::setmask()
 
 void FixForceSpin::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 = mub/hbar; // in rad.THz/T  
 
-  H_field *= gyro; //in rad.THz
-  Ka /= hbar; //in rad.THz
+  H_field *= gyro; // in rad.THz
+  Ka /= hbar; // in rad.THz
 
   if (strstr(update->integrate_style,"respa")) {
     ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
@@ -131,15 +136,15 @@ void FixForceSpin::init()
   if (magstr) {
   magvar = input->variable->find(magstr);
   if (magvar < 0) 
-        error->all(FLERR,"Variable name for fix magnetic field does not exist");
+        error->all(FLERR,"Illegal force/spin command");
   if (!input->variable->equalstyle(magvar))
-        error->all(FLERR,"Variable for fix magnetic field is invalid style");
+        error->all(FLERR,"Illegal force/spin command");
   }
   
   varflag = CONSTANT;
   if (magfieldstyle != CONSTANT) varflag = EQUAL;
  
-  // set magnetic field components once and for all
+  // set magnetic field components
   if (varflag == CONSTANT) set_magneticforce();
   
   memory->create(spi,3,"forcespin:spi"); 
@@ -168,26 +173,23 @@ void FixForceSpin::post_force(int vflag)
   if (varflag != CONSTANT) {
     modify->clearstep_compute();
     modify->addstep_compute(update->ntimestep + 1);
-    set_magneticforce(); //Update value of the mag. field if time-dependent
+    set_magneticforce(); // update mag. field if time-dep.
   }
 
-//  double **x = atom->x;
   double **sp = atom->sp; 
   double *mumag = atom->mumag;
   double **fm = atom->fm; 
   const int nlocal = atom->nlocal;  
   double scalar;
-  
+
   eflag = 0;
   emag = 0.0;
 
   for (int i = 0; i < nlocal; i++) {
     fmi[0] = fmi[1] = fmi[2] = 0.0;
-    //if (style == ZEEMAN) {
     if (zeeman_flag) {
       compute_zeeman(i,fmi);
     }  
-    //if (style == ANISOTROPY) {
     if (aniso_flag) {
       spi[0] = sp[i][0];
       spi[1] = sp[i][1];
@@ -197,7 +199,9 @@ void FixForceSpin::post_force(int vflag)
     fm[i][0] += fmi[0];
     fm[i][1] += fmi[1];
     fm[i][2] += fmi[2];
-  }  
+
+  } 
+
 }
 
 
@@ -226,7 +230,6 @@ void FixForceSpin::post_force_respa(int vflag, int ilevel, int iloop)
 }
 
 /* ---------------------------------------------------------------------- */
-//No acceleration for magnetic EOM, only a "magnetic force" 
 void FixForceSpin::set_magneticforce()
 {
   if (style == ZEEMAN) {
diff --git a/src/SPIN/fix_force_spin.h b/src/SPIN/fix_force_spin.h
index a422abad2c..a31e2f888c 100644
--- a/src/SPIN/fix_force_spin.h
+++ b/src/SPIN/fix_force_spin.h
@@ -42,9 +42,9 @@ class FixForceSpin : public Fix {
   void compute_anisotropy(int, double *, double *);
 
  protected:
-  int style;
+  int style; // style of the magnetic force
   
-  double xmag, ymag, zmag; //Magnetic force
+  double xmag, ymag, zmag; // temp. force variables
   double degree2rad;
   int ilevel_respa;
   int time_origin;
@@ -56,13 +56,16 @@ class FixForceSpin : public Fix {
   int magvar;
   char *magstr;
    
-  double H_field; //Zeeman field intensity and direction
+  // zeeman field intensity and direction
+  double H_field; 
   double Hx, Hy, Hz;
   
-  double Ka; //Magnetic anisotropy intensity and direction
+  // magnetic anisotropy intensity and direction
+  double Ka; 
   double Kax, Kay, Kaz;
 
-  double *spi, *fmi; //Temp var. in compute
+  // temp. spin variables
+  double *spi, *fmi;
   
   void set_magneticforce();
  
@@ -75,18 +78,10 @@ class FixForceSpin : public Fix {
 
 /* ERROR/WARNING messages:
 
-E: Illegal ... command
+E: Illegal force/spin command
 
 Self-explanatory.  Check the input script syntax and compare to the
 documentation for the command.  You can use -echo screen as a
 command-line option when running LAMMPS to see the offending line.
 
-E: Variable name for fix force/spin does not exist
-
-Self-explanatory.
-
-E: Variable for fix force/spin is invalid style
-
-Only equal-style variables can be used.
-
 */
diff --git a/src/SPIN/fix_integration_spin.cpp b/src/SPIN/fix_integration_spin.cpp
index d05f8b39ba..b09aabac05 100644
--- a/src/SPIN/fix_integration_spin.cpp
+++ b/src/SPIN/fix_integration_spin.cpp
@@ -11,6 +11,11 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+/* ------------------------------------------------------------------------
+   Contributing authors: Julien Tranchida (SNL)
+                         Aidan Thompson (SNL)
+------------------------------------------------------------------------- */
+
 #include <math.h>
 #include <stdio.h>
 #include <string.h>
@@ -45,7 +50,7 @@ enum{NONE,SPIN};
 /* ---------------------------------------------------------------------- */
 
 FixIntegrationSpin::FixIntegrationSpin(LAMMPS *lmp, int narg, char **arg) :
-  FixNVE(lmp, narg, arg)
+  Fix(lmp, narg, arg)
 {
 	
   if (narg != 4) error->all(FLERR,"Illegal fix integration/spin command");	
@@ -68,9 +73,11 @@ FixIntegrationSpin::FixIntegrationSpin(LAMMPS *lmp, int narg, char **arg) :
   // error checks
   if (extra == SPIN && !atom->mumag_flag)
     error->all(FLERR,"Fix integration/spin requires spin attribute mumag");
-
+  magpair_flag = 0;
   exch_flag = dmi_flag = me_flag = 0;
+  magforce_flag = 0;
   zeeman_flag = aniso_flag = 0;
+  maglangevin_flag = 0;
   tdamp_flag = temp_flag = 0;
 
   lockpairspin = NULL;
@@ -96,11 +103,25 @@ FixIntegrationSpin::~FixIntegrationSpin(){
 
 /* ---------------------------------------------------------------------- */
 
+int FixIntegrationSpin::setmask()
+{
+  int mask = 0;
+  mask |= INITIAL_INTEGRATE;
+  mask |= FINAL_INTEGRATE;
+  return mask;  
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixIntegrationSpin::init()
 {
-  FixNVE::init();       
+  //FixNVE::init();       
   
+  // set timesteps
+  dtv = update->dt;
+  dtf = 0.5 * update->dt * force->ftm2v;
   dts = update->dt;
+
   memory->create(xi,3,"integrations:xi");
   memory->create(sec,3,"integrations:sec");
   memory->create(rsec,3,"integrations:rsec");
@@ -113,6 +134,7 @@ void FixIntegrationSpin::init()
 
 
   if (strstr(force->pair_style,"pair/spin")) {
+    magpair_flag = 1;
     lockpairspin = (PairSpin *) force->pair;
   } else if (strstr(force->pair_style,"hybrid/overlay")) {
     PairHybrid *lockhybrid = (PairHybrid *) force->pair;
@@ -120,35 +142,48 @@ void FixIntegrationSpin::init()
     int ipair;
     for (ipair = 0; ipair < nhybrid_styles; ipair++) {
       if (strstr(lockhybrid->keywords[ipair],"pair/spin")) {
+        magpair_flag = 1;
 	lockpairspin = (PairSpin *) lockhybrid->styles[ipair];
       }
     }
   } else error->all(FLERR,"Illegal fix integration/spin command");
 
-  // check errors, and handle simple hybrid (not overlay)
+  // check errors, and handle simple hybrid (not overlay), and no pair/spin interaction
 
   int iforce;
-  for (iforce = 0; iforce < modify->nfix; iforce++)
-    if (strstr(modify->fix[iforce]->style,"force/spin")) break;
-  lockforcespin = (FixForceSpin *) modify->fix[iforce]; 
+  for (iforce = 0; iforce < modify->nfix; iforce++) {
+    if (strstr(modify->fix[iforce]->style,"force/spin")) {
+      magforce_flag = 1;
+      lockforcespin = (FixForceSpin *) modify->fix[iforce]; 
+    }
+  }
 
-  for (iforce = 0; iforce < modify->nfix; iforce++)
-    if (strstr(modify->fix[iforce]->style,"langevin/spin")) break;
-  locklangevinspin = (FixLangevinSpin *) modify->fix[iforce]; 
+  for (iforce = 0; iforce < modify->nfix; iforce++) {
+    if (strstr(modify->fix[iforce]->style,"langevin/spin")) {
+      maglangevin_flag = 1;
+      locklangevinspin = (FixLangevinSpin *) modify->fix[iforce]; 
+    } 
+  }
 
-  // set flags for the different magnetic interactions
-  if (lockpairspin->exch_flag == 1) exch_flag = 1;
-  if (lockpairspin->dmi_flag == 1) dmi_flag = 1;
-  if (lockpairspin->me_flag == 1) me_flag = 1;
+  // set flags for the different magnetic interactionsi
+  if (magpair_flag) {
+    if (lockpairspin->exch_flag == 1) exch_flag = 1;
+    if (lockpairspin->dmi_flag == 1) dmi_flag = 1;
+    if (lockpairspin->me_flag == 1) me_flag = 1;
+  }
 
-  if (lockforcespin->zeeman_flag == 1) zeeman_flag = 1;
-  if (lockforcespin->aniso_flag == 1) aniso_flag = 1;
+  if (magforce_flag) { 
+    if (lockforcespin->zeeman_flag == 1) zeeman_flag = 1;
+    if (lockforcespin->aniso_flag == 1) aniso_flag = 1;
+  }
 
-  if (locklangevinspin->tdamp_flag == 1) tdamp_flag = 1;
-  if (locklangevinspin->temp_flag == 1) temp_flag = 1;
+  if (maglangevin_flag) {
+   if (locklangevinspin->tdamp_flag == 1) tdamp_flag = 1;
+   if (locklangevinspin->temp_flag == 1) temp_flag = 1;
+  }
 
   // perform the sectoring if mpi integration
-  if (mpi_flag == 1) sectoring();
+  if (mpi_flag) sectoring();
 
 }
 
@@ -171,11 +206,15 @@ void FixIntegrationSpin::initial_integrate(int vflag)
   int *type = atom->type;
   int *mask = atom->mask;  
 
+  printf("before mechmag, spin 1: [%g,%g,%g] \n",f[1][0],f[1][1],f[1][2]);
+
   // compute and add magneto-mech. force 
-  if (exch_flag) {
-      lockpairspin->compute_magnetomech(0,vflag);
+  if (magpair_flag == 1) { 
+    if (exch_flag) lockpairspin->compute_magnetomech(0,vflag);
   }
 
+  printf("after mechmag, spin 1: [%g,%g,%g] \n",f[1][0],f[1][1],f[1][2]);
+
   // update half v all particles
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
@@ -263,7 +302,7 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
   const int newton_pair = force->newton_pair;
 
   // add test here
-  if (exch_flag == 1 || dmi_flag == 1 || me_flag == 1 ) { 
+  if (magpair_flag) { 
     inum = lockpairspin->list->inum;
     ilist = lockpairspin->list->ilist;
     numneigh = lockpairspin->list->numneigh;
@@ -286,9 +325,6 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
   // force computation for spin i
   i = ilist[ii];
   
-  // clear atom i 
-  fm[i][0] = fm[i][1] = fm[i][2] = 0.0;
-  
   spi[0] = sp[i][0];
   spi[1] = sp[i][1];
   spi[2] = sp[i][2];
@@ -316,50 +352,54 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
     itype = type[ii];
     jtype = type[j];
 
-    if (exch_flag) {
-      cut_ex_2 = (lockpairspin->cut_spin_exchange[itype][jtype])*(lockpairspin->cut_spin_exchange[itype][jtype]);
-      if (rsq <= cut_ex_2) {
-        lockpairspin->compute_exchange(i,j,rsq,fmi,fmj,spi,spj);
-      }  
-    }
-
-    if (dmi_flag) {
-      cut_dmi_2 = (lockpairspin->cut_spin_dmi[itype][jtype])*(lockpairspin->cut_spin_dmi[itype][jtype]);
-      if (rsq <= cut_dmi_2) {
-        lockpairspin->compute_dmi(i,j,fmi,fmj,spi,spj);
-      }  
-    }
-
-    if (me_flag) {
-      cut_me_2 = (lockpairspin->cut_spin_me[itype][jtype])*(lockpairspin->cut_spin_me[itype][jtype]);
-      if (rsq <= cut_me_2) {
-        lockpairspin->compute_me(i,j,fmi,fmj,spi,spj);
-      }  
-    }
-
-  }
-
-  // post force
-  if (zeeman_flag) {
-    lockforcespin->compute_zeeman(i,fmi);
-  }
-
-  if (aniso_flag) {
-    spi[0] = sp[i][0];
-    spi[1] = sp[i][1];                                       
-    spi[2] = sp[i][2]; 
-    lockforcespin->compute_anisotropy(i,spi,fmi);
-  }
-
-  if (tdamp_flag) {
-    locklangevinspin->add_tdamping(spi,fmi);   
-  }
-
-  if (temp_flag) {
-    locklangevinspin->add_temperature(fmi);
+    if (magpair_flag) { // mag. pair inter.
+      double temp_cut;
+      if (exch_flag) { // exchange
+        temp_cut = lockpairspin->cut_spin_exchange[itype][jtype];
+	cut_ex_2 = temp_cut*temp_cut;
+        if (rsq <= cut_ex_2) { 
+          lockpairspin->compute_exchange(i,j,rsq,fmi,fmj,spi,spj);
+        }  
+      }
+      if (dmi_flag) { // dmi
+	temp_cut = lockpairspin->cut_spin_dmi[itype][jtype];
+	cut_dmi_2 = temp_cut*temp_cut;
+	if (rsq <= cut_dmi_2) {
+	  lockpairspin->compute_dmi(i,j,fmi,fmj,spi,spj);
+	}
+      }
+      if (me_flag) { // me
+	temp_cut = lockpairspin->cut_spin_me[itype][jtype];
+	cut_me_2 = temp_cut*temp_cut;
+	if (rsq <= cut_me_2) {
+	  lockpairspin->compute_me(i,j,fmi,fmj,spi,spj);
+	}
+      }
+    }  
   } 
- 
-  // replace the magnetic force by its new value
+
+  if (magforce_flag) { // mag. forces
+    if (zeeman_flag) { // zeeman
+      lockforcespin->compute_zeeman(i,fmi);
+    }
+    if (aniso_flag) { // aniso
+      spi[0] = sp[i][0];
+      spi[1] = sp[i][1];                                       
+      spi[2] = sp[i][2];
+      lockforcespin->compute_anisotropy(i,spi,fmi);
+    }
+  }
+
+  if (maglangevin_flag) { // mag. langevin
+    if (tdamp_flag) { // trans. damping
+      locklangevinspin->add_tdamping(spi,fmi);   
+    }
+    if (temp_flag) { // temp.
+      locklangevinspin->add_temperature(fmi);
+    } 
+  }
+
+  // replace the mag. force i by its new value
   fm[i][0] = fmi[0];
   fm[i][1] = fmi[1];
   fm[i][2] = fmi[2];
@@ -397,7 +437,7 @@ void FixIntegrationSpin::sectoring()
   nsectors = sec[0]*sec[1]*sec[2];
 
   if (mpi_flag == 1 && nsectors != 8) 
-    error->all(FLERR,"System too small for sectoring operation");
+    error->all(FLERR,"Illegal sectoring operation");
 
   rsec[0] = rsx;
   rsec[1] = rsy;
@@ -406,17 +446,6 @@ void FixIntegrationSpin::sectoring()
   if (sec[1] == 2) rsec[1] = rsy/2.0;
   if (sec[2] == 2) rsec[2] = rsz/2.0;
 
-  /*
-  if (2.0 * rv >= rsx && sec[0] >= 2)
-    error->all(FLERR,"Illegal number of sectors"); 
-
-  if (2.0 * rv >= rsy && sec[1] >= 2)
-    error->all(FLERR,"Illegal number of sectors"); 
-
-  if (2.0 * rv >= rsz && sec[2] >= 2)
-    error->all(FLERR,"Illegal number of sectors"); 
-*/
-
 }
 
 /* ---------------------------------------------------------------------- */
@@ -501,9 +530,9 @@ void FixIntegrationSpin::final_integrate()
   int *type = atom->type;
   int *mask = atom->mask; 
 
-  // compute and add exchange magneto-mech. force 
-  if (exch_flag) {
-      lockpairspin->compute_magnetomech(0,0);
+  // compute and add magneto-mech. force 
+  if (magpair_flag == 1) { 
+    if (exch_flag) lockpairspin->compute_magnetomech(0,0);
   }
 
   // update half v for all particles
diff --git a/src/SPIN/fix_integration_spin.h b/src/SPIN/fix_integration_spin.h
index ff85dedc03..96f0035919 100644
--- a/src/SPIN/fix_integration_spin.h
+++ b/src/SPIN/fix_integration_spin.h
@@ -17,43 +17,57 @@ FixStyle(integration/spin,FixIntegrationSpin)
 
 #else
 
-#ifndef LMP_FIX_NVE_SPIN_H
-#define LMP_FIX_NVE_SPIN_H
+#ifndef LMP_FIX_INTEGRATION_SPIN_H
+#define LMP_FIX_INTEGRATION_SPIN_H
 
-#include "fix_nve.h"
+#include "fix.h"
 
 namespace LAMMPS_NS {
 
-class FixIntegrationSpin : public FixNVE {
+class FixIntegrationSpin : public Fix {
 	
  public:
   FixIntegrationSpin(class LAMMPS *, int, char **);
   virtual ~FixIntegrationSpin();
+  int setmask();
   void init();
   virtual void initial_integrate(int);
-  void AdvanceSingleSpin(int, double, double **, double **);
   virtual void final_integrate();
-  void ComputeInteractionsSpin(int);   
 
-  void sectoring();
+  // compute and advance single spin
+  void ComputeInteractionsSpin(int);   
+  void AdvanceSingleSpin(int, double, double **, double **);
+
+  // sectoring operations
+  void sectoring(); 
   int coords2sector(double *);
 
  protected:
   int extra, mpi_flag;
-  double dts;
+
+  // vel., force, and spin timesteps
+  double dtv,dtf,dts;
   
+  // mag. interaction flags
+  int magpair_flag;
   int exch_flag, dmi_flag, me_flag;
+  int magforce_flag;
   int zeeman_flag, aniso_flag;
+  int maglangevin_flag;
   int tdamp_flag, temp_flag;
 
+  // pointers to mag. interaction classes
   class PairHybrid *lockhybrid; 
   class PairSpin *lockpairspin;
   class FixForceSpin *lockforcespin;
   class FixLangevinSpin *locklangevinspin; 
 
-  double *spi, *spj, *fmi, *fmj; //Temp var.
+  // temporary variables
   double *xi;
+  double *spi, *spj;
+  double *fmi, *fmj; 
  
+  // sectoring variables
   int nsectors;
   int *sec, *seci;
   double *rsec;
@@ -67,27 +81,19 @@ class FixIntegrationSpin : public FixNVE {
 
 /* ERROR/WARNING messages:
 
-E: Illegal ... command
+E: Illegal fix integration/spin command
 
-Self-explanatory.  Check the input script syntax and compare to the
+Self-explanatory.  Check the input script syntax and compare to the 
 documentation for the command.  You can use -echo screen as a
 command-line option when running LAMMPS to see the offending line.
 
-E: Fix nve/sphere requires atom style sphere
+E: Fix integration/spin requires spin attribute mumag
 
-Self-explanatory.
+An atom/spin style with this attribute is needed.
 
-E: Fix nve/sphere update dipole requires atom attribute mu
-
-An atom style with this attribute is needed.
-
-E: Fix nve/sphere requires extended particles
-
-This fix can only be used for particles of a finite size.
- 
-E: Fix nve/sphere dlm must be used with update dipole
- 
-The DLM algorithm can only be used in conjunction with update dipole.
+E: Illegal sectoring operation
 
+The number of processes does not match the size of the system. 
+See the documentation of the sectoring method.
 
 */
diff --git a/src/SPIN/fix_langevin_spin.cpp b/src/SPIN/fix_langevin_spin.cpp
index 53caa7391e..4f59f8f4d7 100644
--- a/src/SPIN/fix_langevin_spin.cpp
+++ b/src/SPIN/fix_langevin_spin.cpp
@@ -12,8 +12,8 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing authors: Carolyn Phillips (U Mich), reservoir energy tally
-                         Aidan Thompson (SNL) GJF formulation
+   Contributing authors: Julien Tranchida (SNL)
+                         Aidan Thompson (SNL) 
 ------------------------------------------------------------------------- */
 
 #include <mpi.h>
@@ -51,7 +51,7 @@ using namespace MathConst;
 FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg), id_temp(NULL), random(NULL)
 {
-  if (narg != 7) error->all(FLERR,"Illegal fix langevin/spin command");
+  if (narg != 7) error->all(FLERR,"Illegal langevin/spin command");
 
   dynamic_group_allow = 1;
   scalar_flag = 1;
@@ -65,7 +65,7 @@ FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) :
   seed = force->inumeric(FLERR,arg[6]);
 
   if (alpha_t < 0.0) {
-    error->all(FLERR,"Illegal fix/langevin/spin command");
+    error->all(FLERR,"Illegal langevin/spin command");
   } else if (alpha_t == 0.0) {
     tdamp_flag = 0;
   } else {
@@ -73,7 +73,7 @@ FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) :
   }
 
   if (alpha_l < 0.0) {
-    error->all(FLERR,"Illegal fix/langevin/spin command");
+    error->all(FLERR,"Illegal langevin/spin command");
   } else if (alpha_l == 0.0) {
     ldamp_flag = 0;
   } else {
@@ -81,7 +81,7 @@ FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) :
   }
   
   if (temp < 0.0) {
-    error->all(FLERR,"Illegal fix/langevin/spin command");
+    error->all(FLERR,"Illegal langevin/spin command");
   } else if (temp == 0.0) {
     temp_flag = 0;
   } else {
@@ -201,8 +201,9 @@ void FixLangevinSpin::add_tdamping(double *spi, double *fmi)
   double cpx = fmi[1]*spi[2] - fmi[2]*spi[1];
   double cpy = fmi[2]*spi[0] - fmi[0]*spi[2];
   double cpz = fmi[0]*spi[1] - fmi[1]*spi[0];
-		
-  fmi[0] -= alpha_t*cpx;//Taking the damping value away
+	
+  // taking away the transverse damping 
+  fmi[0] -= alpha_t*cpx;
   fmi[1] -= alpha_t*cpy;
   fmi[2] -= alpha_t*cpz;
 }
@@ -212,7 +213,8 @@ void FixLangevinSpin::add_temperature(double *fmi)
 {
 //#define GAUSSIAN_R
 #if defined GAUSSIAN_R
-  double rx = sigma*random->gaussian();//Drawing random dist
+  // drawing gausian random dist
+  double rx = sigma*random->gaussian();
   double ry = sigma*random->gaussian();
   double rz = sigma*random->gaussian();
 #else 
@@ -221,11 +223,13 @@ void FixLangevinSpin::add_temperature(double *fmi)
   double rz = sigma*(random->uniform() - 0.5);
 #endif
 
-  fmi[0] += rx;//Adding random field
+  // adding the random field 
+  fmi[0] += rx; 
   fmi[1] += ry;
   fmi[2] += rz;
-                
-  fmi[0] *= Gil_factor;//Multiplying by Gilbert's prefactor 
+               
+  // adding Gilbert's prefactor 
+  fmi[0] *= Gil_factor; 
   fmi[1] *= Gil_factor; 
   fmi[2] *= Gil_factor; 
 
diff --git a/src/SPIN/fix_langevin_spin.h b/src/SPIN/fix_langevin_spin.h
index a6b9bc5df3..2fe6c6fab5 100644
--- a/src/SPIN/fix_langevin_spin.h
+++ b/src/SPIN/fix_langevin_spin.h
@@ -33,15 +33,19 @@ class FixLangevinSpin : public Fix {
   void setup(int);
   virtual void post_force(int);
   void post_force_respa(int, int, int);
+  // add transverse damping and temperature
   void add_tdamping(double *, double *);
   void add_temperature(double *);
+  // associated flags
   int tdamp_flag, ldamp_flag, temp_flag;
 
  protected:
   double *spi, *fmi;
-  double alpha_t, alpha_l; //Transverse and long. damping value
-  double dts,temp,D,sigma;//timestep, temp, noise  
-  double Gil_factor;//Gilbert's prefactor 
+  // transverse and longitudinal damping coeff.
+  double alpha_t, alpha_l; 
+  // timestep, temperature, noise intensity
+  double dts,temp,D,sigma;   
+  double Gil_factor; 
  
   char *id_temp;
   class Compute *temperature;
@@ -60,7 +64,7 @@ class FixLangevinSpin : public Fix {
 
 /* ERROR/WARNING messages:
 
-E: Illegal ... command
+E: Illegal langevin/spin command
 
 Self-explanatory.  Check the input script syntax and compare to the
 documentation for the command.  You can use -echo screen as a
diff --git a/src/SPIN/pair_spin.cpp b/src/SPIN/pair_spin.cpp
index 689e3f8844..b62fd867de 100755
--- a/src/SPIN/pair_spin.cpp
+++ b/src/SPIN/pair_spin.cpp
@@ -11,6 +11,11 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+/* ------------------------------------------------------------------------
+   Contributing authors: Julien Tranchida (SNL)
+                         Aidan Thompson (SNL)
+------------------------------------------------------------------------- */
+
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
@@ -48,9 +53,10 @@ PairSpin::~PairSpin()
     memory->destroy(setflag);
     
     memory->destroy(cut_spin_exchange);
-    memory->destroy(J_1);
-    memory->destroy(J_2);
-    memory->destroy(J_2);  
+    memory->destroy(J1_mag);
+    memory->destroy(J1_mech);
+    memory->destroy(J2);
+    memory->destroy(J3);  
     
     memory->destroy(cut_spin_dmi);
     memory->destroy(DM);
@@ -70,6 +76,7 @@ PairSpin::~PairSpin()
     memory->destroy(fj);
     memory->destroy(fmi);
     memory->destroy(fmj);
+    memory->destroy(del);
 
     memory->destroy(cutsq);
   }
@@ -85,7 +92,7 @@ void PairSpin::compute_magnetomech(int eflag, int vflag)
   double xtmp,ytmp,ztmp;
   double fix,fiy,fiz,fjx,fjy,fjz;
   double cut_ex_2,cut_dmi_2,cut_me_2;
-  double rsq,rd,delx,dely,delz;
+  double rsq,rd;
   int *ilist,*jlist,*numneigh,**firstneigh;  
 
   evdwl = ecoul = 0.0;
@@ -119,21 +126,21 @@ void PairSpin::compute_magnetomech(int eflag, int vflag)
     spi[1] = sp[i][1]; 
     spi[2] = sp[i][2]; 
 
-    //Loop on Neighbors
+    // loop on neighbors
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;
       spj[0] = sp[j][0]; 
       spj[1] = sp[j][1]; 
       spj[2] = sp[j][2]; 
-
       fi[0] = fi[1] = fi[2] = 0.0;
       fj[0] = fj[1] = fj[2] = 0.0;
-     
-      delx = xtmp - x[j][0];
-      dely = ytmp - x[j][1];
-      delz = ztmp - x[j][2];
-      rsq = delx*delx + dely*dely + delz*delz;  //square or inter-atomic distance
+      del[0] = del[1] = del[2] = 0.0;
+
+      del[0] = x[j][0] - xtmp;
+      del[1] = x[j][1] - ytmp;
+      del[2] = x[j][2] - ztmp;
+      rsq = del[0]*del[0] + del[1]*del[1] + del[2]*del[2];  //square or inter-atomic distance
       itype = type[i];
       jtype = type[j];
 
@@ -141,7 +148,7 @@ void PairSpin::compute_magnetomech(int eflag, int vflag)
       if (exch_flag) {
         cut_ex_2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
         if (rsq <= cut_ex_2) {
-          compute_exchange_mech(i,j,rsq,fi,fj,spi,spj);   
+          compute_exchange_mech(i,j,rsq,del,fi,fj,spi,spj);   
         }
       }
 
@@ -204,7 +211,7 @@ void PairSpin::compute(int eflag, int vflag)
     spi[1] = sp[i][1]; 
     spi[2] = sp[i][2]; 
    
-    //Loop on Neighbors
+    // loop on neighbors
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;
@@ -218,25 +225,27 @@ void PairSpin::compute(int eflag, int vflag)
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
-      rsq = delx*delx + dely*dely + delz*delz;  //square or inter-atomic distance
+      // square or inter-atomic distance
+      rsq = delx*delx + dely*dely + delz*delz; 
+
       itype = type[i];
       jtype = type[j];
 
-      //Exchange interaction
+      // exchange interaction
       if (exch_flag) {
         cut_ex_2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
         if (rsq <= cut_ex_2) {
           compute_exchange(i,j,rsq,fmi,fmj,spi,spj);   
         }
       }
-      //DM interaction
+      // dm interaction
       if (dmi_flag){
         cut_dmi_2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];
         if (rsq <= cut_dmi_2){
           compute_dmi(i,j,fmi,fmj,spi,spj);
         } 
       }
-      //ME interaction
+      // me interaction
       if (me_flag){
         cut_me_2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype];
         if (rsq <= cut_me_2){
@@ -249,12 +258,12 @@ void PairSpin::compute(int eflag, int vflag)
       fm[i][2] += fmi[2];
 
       if (newton_pair || j < nlocal) {
-          fm[j][0] += fmj[0];	 
-          fm[j][1] += fmj[1];	  	  
-          fm[j][2] += fmj[2];
-          }
+        fm[j][0] += fmj[0];	 
+        fm[j][1] += fmj[1];	  	  
+        fm[j][2] += fmj[2];
       }
-  }
+    }
+  }  
 
   if (vflag_fdotr) virial_fdotr_compute();
 }
@@ -264,14 +273,13 @@ void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi,  double *
 {
   int *type = atom->type;  
   int itype, jtype;
-  //double dmix,dmiy,dmiz;	
   double Jex, ra;
   itype = type[i];
   jtype = type[j];
           
-  ra = rsq/J_3[itype][jtype]/J_3[itype][jtype]; 
-  Jex = 4.0*J_1[itype][jtype]*ra;
-  Jex *= (1.0-J_2[itype][jtype]*ra);
+  ra = rsq/J3[itype][jtype]/J3[itype][jtype]; 
+  Jex = 4.0*J1_mag[itype][jtype]*ra;
+  Jex *= (1.0-J2[itype][jtype]*ra);
   Jex *= exp(-ra);
 
   fmi[0] += Jex*spj[0];
@@ -285,29 +293,29 @@ void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi,  double *
 }
 
 /* ---------------------------------------------------------------------- */
-void PairSpin::compute_exchange_mech(int i, int j, double rsq, double *fmi,  double *fmj, double *spi, double *spj)
+void PairSpin::compute_exchange_mech(int i, int j, double rsq, double *del, double *fi,  double *fj, double *spi, double *spj)
 {
   int *type = atom->type;  
   int itype, jtype;
-  //double dmix,dmiy,dmiz;	
-  double Jex_mech, ra, rr;
+  double Jex, Jex_mech, ra, rr;
   itype = type[i];
   jtype = type[j];
-          
-  ra = rsq/J_3[itype][jtype]/J_3[itype][jtype]; 
-  rr = sqrt(rsq)/J_3[itype][jtype]/J_3[itype][jtype];
-  Jex_mech = 1.0-2.0*J_2[itype][jtype]*ra;
-  Jex_mech -= ra*(1.0-J_2[itype][jtype]*ra);
-  Jex_mech *= 8.0*J_1[itype][jtype]*rr*exp(-ra);
+  Jex = J1_mech[itype][jtype];        
 
+  ra = rsq/J3[itype][jtype]/J3[itype][jtype]; 
+  rr = sqrt(rsq)/J3[itype][jtype]/J3[itype][jtype];
 
-  fi[0] += Jex_mech*spi[0]*spj[0];
-  fi[1] += Jex_mech*spi[1]*spj[1];
-  fi[2] += Jex_mech*spi[2]*spj[2];
+  Jex_mech = 1.0-ra-J2[itype][jtype]*ra*(2.0-ra);
+  Jex_mech *= 8.0*Jex*rr*exp(-ra);
+  Jex_mech *= (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]); 
+
+  fi[0] += Jex_mech*del[0];
+  fi[1] += Jex_mech*del[1];
+  fi[2] += Jex_mech*del[2];
           
-  fj[0] += Jex_mech*spi[0]*spj[0];
-  fj[1] += Jex_mech*spi[1]*spj[1];
-  fj[2] += Jex_mech*spi[2]*spj[2];
+  fj[0] -= Jex_mech*del[0];
+  fj[1] -= Jex_mech*del[1];
+  fj[2] -= Jex_mech*del[2];
   
 }
 
@@ -386,9 +394,10 @@ void PairSpin::allocate()
       setflag[i][j] = 0;
       
   memory->create(cut_spin_exchange,n+1,n+1,"pair:cut_spin_exchange");
-  memory->create(J_1,n+1,n+1,"pair:J_1");
-  memory->create(J_2,n+1,n+1,"pair:J_2");  
-  memory->create(J_3,n+1,n+1,"pair:J_3");
+  memory->create(J1_mag,n+1,n+1,"pair:J1_mag");
+  memory->create(J1_mech,n+1,n+1,"pair:J1_mech");
+  memory->create(J2,n+1,n+1,"pair:J2");  
+  memory->create(J3,n+1,n+1,"pair:J3");
 
   memory->create(cut_spin_dmi,n+1,n+1,"pair:cut_spin_dmi");
   memory->create(DM,n+1,n+1,"pair:DM");
@@ -404,10 +413,11 @@ void PairSpin::allocate()
  
   memory->create(spi,3,"pair:spi");
   memory->create(spj,3,"pair:spj");
-  memory->create(fi,3,"pair:fmi");
-  memory->create(fj,3,"pair:fmj");
+  memory->create(fi,3,"pair:fi");
+  memory->create(fj,3,"pair:fj");
   memory->create(fmi,3,"pair:fmi");
   memory->create(fmj,3,"pair:fmj");
+  memory->create(del,3,"pair:del");
  
   memory->create(cutsq,n+1,n+1,"pair:cutsq");  
   
@@ -461,17 +471,18 @@ void PairSpin::coeff(int narg, char **arg)
     force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
     
     const double rij = force->numeric(FLERR,arg[3]);
-    const double J1 = (force->numeric(FLERR,arg[4]))/hbar;
-    const double J2 = force->numeric(FLERR,arg[5]);  
-    const double J3 = force->numeric(FLERR,arg[6]); 
+    const double j1 = (force->numeric(FLERR,arg[4]));
+    const double j2 = force->numeric(FLERR,arg[5]);  
+    const double j3 = force->numeric(FLERR,arg[6]); 
   
     int count = 0;
     for (int i = ilo; i <= ihi; i++) {
       for (int j = MAX(jlo,i); j <= jhi; j++) {
         cut_spin_exchange[i][j] = rij;   
-        J_1[i][j] = J1;
-        J_2[i][j] = J2;
-        J_3[i][j] = J3;
+        J1_mag[i][j] = j1/hbar;
+        J1_mech[i][j] = j1;
+        J2[i][j] = j2;
+        J3[i][j] = j3;
         setflag[i][j] = 1;
         count++;
       }
@@ -592,9 +603,10 @@ void PairSpin::write_restart(FILE *fp)
       fwrite(&setflag[i][j],sizeof(int),1,fp);
       if (setflag[i][j]) {
         if (exch_flag){
-          fwrite(&J_1[i][j],sizeof(double),1,fp);
-          fwrite(&J_2[i][j],sizeof(double),1,fp);
-          fwrite(&J_3[i][j],sizeof(double),1,fp);
+          fwrite(&J1_mag[i][j],sizeof(double),1,fp);
+          fwrite(&J1_mech[i][j],sizeof(double),1,fp);
+          fwrite(&J2[i][j],sizeof(double),1,fp);
+          fwrite(&J3[i][j],sizeof(double),1,fp);
           fwrite(&cut_spin_exchange[i][j],sizeof(double),1,fp);
         }
         if (dmi_flag) {
@@ -633,14 +645,16 @@ void PairSpin::read_restart(FILE *fp)
       MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
       if (setflag[i][j]) {
         if (me == 0) {
-          fread(&J_1[i][j],sizeof(double),1,fp);
-          fread(&J_2[i][j],sizeof(double),1,fp);
-          fread(&J_2[i][j],sizeof(double),1,fp);
+          fread(&J1_mag[i][j],sizeof(double),1,fp);
+          fread(&J1_mech[i][j],sizeof(double),1,fp);
+          fread(&J2[i][j],sizeof(double),1,fp);
+          fread(&J2[i][j],sizeof(double),1,fp);
           fread(&cut_spin_exchange[i][j],sizeof(double),1,fp);
         }
-        MPI_Bcast(&J_1[i][j],1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&J_2[i][j],1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&J_3[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&J1_mag[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&J1_mech[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&J2[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&J3[i][j],1,MPI_DOUBLE,0,world);
         MPI_Bcast(&cut_spin_exchange[i][j],1,MPI_DOUBLE,0,world);
       }
     }
diff --git a/src/SPIN/pair_spin.h b/src/SPIN/pair_spin.h
index 36469d5acf..7ef78f2767 100755
--- a/src/SPIN/pair_spin.h
+++ b/src/SPIN/pair_spin.h
@@ -41,7 +41,7 @@ class PairSpin : public Pair {
   void read_restart_settings(FILE *);
   
   void compute_exchange(int, int, double, double *, double *,double *, double *);
-  void compute_exchange_mech(int, int, double, double *, double *,double *, double *);
+  void compute_exchange_mech(int, int, double, double *, double *, double *,double *, double *);
   void compute_dmi(int, int, double *, double *, double *, double *);
   void compute_me(int, int, double *, double *, double *, double *);  
  
@@ -56,18 +56,19 @@ class PairSpin : public Pair {
   double **cut_spin_me;       // cutoff distance me 
 
  protected: 
-  double **J_1, **J_2, **J_3; // exchange coeffs Jij
-                              // J1 in eV, J2 adim and J3 in Ang
-  double **DM;
-  double **v_dmx, **v_dmy, **v_dmz;// DMI coeffs
-                                   // DM int. in eV, v direction
+  double **J1_mag, **J1_mech; // exchange coeffs Jij
+  double **J2, **J3; // J1 in eV, J2 adim, J3 in Ang
 
-  double **ME;                     // ME in eV
-  double **v_mex, **v_mey, **v_mez;// ME direction
+  double **DM;                     // dmi coeff in eV 
+  double **v_dmx, **v_dmy, **v_dmz;// dmi direction
 
-  double *spi, *spj;
+  double **ME;                     // me coeff in eV
+  double **v_mex, **v_mey, **v_mez;// me direction
+
+  double *spi, *spj; // temp. spin vals. in compute
   double *fi, *fj;   // temp. mech. forces  in compute
   double *fmi, *fmj; // temp. mag. forces in compute
+  double *del;       // inter atomic distances
 
   void allocate();
 };
@@ -79,7 +80,7 @@ class PairSpin : public Pair {
 
 /* ERROR/WARNING messages:
 
-E: Incorrect args in pair_style command
+E: Incorrect args in pair_spin command
 
 Self-explanatory.