diff --git a/doc/src/read_data.txt b/doc/src/read_data.txt
index 62112ea886..5fb91ba797 100644
--- a/doc/src/read_data.txt
+++ b/doc/src/read_data.txt
@@ -565,6 +565,7 @@ molecular: atom-ID molecule-ID atom-type x y z
 peri: atom-ID atom-type volume density x y z
 smd: atom-ID atom-type molecule volume mass kernel-radius contact-radius x y z
 sphere: atom-ID atom-type diameter density x y z
+spin: atom-ID atom-type sp x y z spx spy spz
 template: atom-ID molecule-ID template-index template-atom atom-type x y z
 tri: atom-ID molecule-ID atom-type triangleflag density x y z
 wavepacket: atom-ID atom-type charge spin eradius etag cs_re cs_im x y z
@@ -595,6 +596,8 @@ mux,muy,muz = components of dipole moment of atom (dipole units)
 q = charge on atom (charge units)
 rho = density (need units) for SPH particles
 spin = electron spin (+1/-1), 0 = nuclei, 2 = fixed-core, 3 = pseudo-cores (i.e. ECP)
+sp = norm of magnetic spin of atom (in number of Bohr magnetons)
+spx,spy,spz = components of magnetic spin of atom (adim normalized vector)
 template-atom = which atom within a template molecule the atom is
 template-index = which molecule within the molecule template the atom is part of
 theta = internal temperature of a DPD particle
diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp
index 5b4173509f..8cc38dbc15 100644
--- a/src/SPIN/atom_vec_spin.cpp
+++ b/src/SPIN/atom_vec_spin.cpp
@@ -263,10 +263,10 @@ int AtomVecSpin::pack_comm_hybrid(int n, int *list, double *buf)
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
+    buf[m++] = sp[j][3];
     buf[m++] = sp[j][0];
     buf[m++] = sp[j][1];
     buf[m++] = sp[j][2];
-    buf[m++] = sp[j][3];
   }
   return m;
 }
@@ -321,10 +321,10 @@ int AtomVecSpin::unpack_comm_hybrid(int n, int first, double *buf)
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
+    sp[i][3] = buf[m++];
     sp[i][0] = buf[m++];
     sp[i][1] = buf[m++];
     sp[i][2] = buf[m++];
-    sp[i][3] = buf[m++];
   }
   return m;
 }
@@ -524,10 +524,10 @@ int AtomVecSpin::pack_border_hybrid(int n, int *list, double *buf)
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
+    buf[m++] = sp[j][3];
     buf[m++] = sp[j][0];
     buf[m++] = sp[j][1];
     buf[m++] = sp[j][2];
-    buf[m++] = sp[j][3];
   }
 
   return m;
@@ -603,10 +603,10 @@ int AtomVecSpin::unpack_border_hybrid(int n, int first, double *buf)
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
+    sp[i][3] = buf[m++];
     sp[i][0] = buf[m++];
     sp[i][1] = buf[m++];
     sp[i][2] = buf[m++];
-    sp[i][3] = buf[m++];
   }
 
   return m;
@@ -894,10 +894,10 @@ void AtomVecSpin::pack_data(double **buf)
 
 int AtomVecSpin::pack_data_hybrid(int i, double *buf)
 {
-  buf[0] = sp[i][0];
-  buf[1] = sp[i][1];
-  buf[2] = sp[i][2];
-  buf[3] = sp[i][3];
+  buf[0] = sp[i][3];
+  buf[1] = sp[i][0];
+  buf[2] = sp[i][1];
+  buf[3] = sp[i][2];
   return 4;
 }
 
@@ -924,7 +924,7 @@ void AtomVecSpin::write_data(FILE *fp, int n, double **buf)
 
 int AtomVecSpin::write_data_hybrid(FILE *fp, double *buf)
 {
-  fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e",buf[0],buf[1],buf[2],buf[3],buf[4]);
+  fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e",buf[0],buf[1],buf[2],buf[3]);
   return 4;
 }
 
diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp
index 595ddb0cc2..b992f17cc3 100644
--- a/src/SPIN/fix_nve_spin.cpp
+++ b/src/SPIN/fix_nve_spin.cpp
@@ -312,29 +312,37 @@ void FixNVESpin::initial_integrate(int /*vflag*/)
       comm->forward_comm();
       int i = stack_foot[j];
       while (i >= 0) {
-        ComputeInteractionsSpin(i);
-        AdvanceSingleSpin(i);
-        i = forward_stacks[i];
+        if (mask[i] & groupbit) {
+          ComputeInteractionsSpin(i);
+          AdvanceSingleSpin(i);
+          i = forward_stacks[i];
+	}
       }
     }
     for (int j = nsectors-1; j >= 0; j--) {     // advance quarter s for nlocal
       comm->forward_comm();
       int i = stack_head[j];
       while (i >= 0) {
-        ComputeInteractionsSpin(i);
-        AdvanceSingleSpin(i);
-        i = backward_stacks[i];
+        if (mask[i] & groupbit) {
+          ComputeInteractionsSpin(i);
+          AdvanceSingleSpin(i);
+          i = backward_stacks[i];
+	}
       }
     }
   } else if (sector_flag == 0) {                // serial seq. update
     comm->forward_comm();                       // comm. positions of ghost atoms
     for (int i = 0; i < nlocal; i++){           // advance quarter s for nlocal
-      ComputeInteractionsSpin(i);
-      AdvanceSingleSpin(i);
+      if (mask[i] & groupbit) {
+        ComputeInteractionsSpin(i);
+        AdvanceSingleSpin(i);
+      }
     }
     for (int i = nlocal-1; i >= 0; i--){        // advance quarter s for nlocal
-      ComputeInteractionsSpin(i);
-      AdvanceSingleSpin(i);
+      if (mask[i] & groupbit) {
+        ComputeInteractionsSpin(i);
+        AdvanceSingleSpin(i);
+      }
     }
   } else error->all(FLERR,"Illegal fix NVE/spin command");
 
@@ -357,29 +365,37 @@ void FixNVESpin::initial_integrate(int /*vflag*/)
       comm->forward_comm();
       int i = stack_foot[j];
       while (i >= 0) {
-        ComputeInteractionsSpin(i);
-        AdvanceSingleSpin(i);
-        i = forward_stacks[i];
+        if (mask[i] & groupbit) {
+          ComputeInteractionsSpin(i);
+          AdvanceSingleSpin(i);
+          i = forward_stacks[i];
+	}
       }
     }
     for (int j = nsectors-1; j >= 0; j--) {     // advance quarter s for nlocal
       comm->forward_comm();
       int i = stack_head[j];
       while (i >= 0) {
-        ComputeInteractionsSpin(i);
-        AdvanceSingleSpin(i);
-        i = backward_stacks[i];
+        if (mask[i] & groupbit) {
+          ComputeInteractionsSpin(i);
+          AdvanceSingleSpin(i);
+          i = backward_stacks[i];
+	}
       }
     }
   } else if (sector_flag == 0) {                // serial seq. update
     comm->forward_comm();                       // comm. positions of ghost atoms
     for (int i = 0; i < nlocal; i++){           // advance quarter s for nlocal-1
-      ComputeInteractionsSpin(i);
-      AdvanceSingleSpin(i);
+      if (mask[i] & groupbit) {
+        ComputeInteractionsSpin(i);
+        AdvanceSingleSpin(i);
+      }
     }
     for (int i = nlocal-1; i >= 0; i--){        // advance quarter s for nlocal-1
-      ComputeInteractionsSpin(i);
-      AdvanceSingleSpin(i);
+      if (mask[i] & groupbit) {
+        ComputeInteractionsSpin(i);
+        AdvanceSingleSpin(i);
+      }
     }
   } else error->all(FLERR,"Illegal fix NVE/spin command");
 
@@ -613,11 +629,11 @@ void FixNVESpin::AdvanceSingleSpin(int i)
 
   // renormalization (check if necessary)
 
-  //msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
-  //scale = 1.0/sqrt(msq);
-  //sp[i][0] *= scale;
-  //sp[i][1] *= scale;
-  //sp[i][2] *= scale;
+  // msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
+  // scale = 1.0/sqrt(msq);
+  // sp[i][0] *= scale;
+  // sp[i][1] *= scale;
+  // sp[i][2] *= scale;
 
   // comm. sp[i] to atoms with same tag (for serial algo)
 
diff --git a/src/SPIN/min_spin.cpp b/src/SPIN/min_spin.cpp
index 2bddc110e7..dda414ff97 100644
--- a/src/SPIN/min_spin.cpp
+++ b/src/SPIN/min_spin.cpp
@@ -287,12 +287,6 @@ void MinSpin::advance_spins(double dts)
 
     // renormalization (check if necessary)
 
-    msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
-    scale = 1.0/sqrt(msq);
-    sp[i][0] *= scale;
-    sp[i][1] *= scale;
-    sp[i][2] *= scale;
-
     // no comm. to atoms with same tag
     // because no need for simplecticity
   }