diff --git a/examples/SPIN/gneb_bfo/final.iron_spin b/examples/SPIN/gneb_bfo/final.iron_spin
index a921287ccb..aa1cbae770 100644
--- a/examples/SPIN/gneb_bfo/final.iron_spin
+++ b/examples/SPIN/gneb_bfo/final.iron_spin
@@ -1,36 +1,36 @@
 32
-2.2000000000000002e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 1.4332499999999999e+00 1.4332499999999999e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 2.8664999999999998e+00 0.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 4.2997499999999995e+00 1.4332499999999999e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 5.7329999999999997e+00 0.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 7.1662499999999998e+00 1.4332499999999999e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 8.5994999999999990e+00 0.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 1.0032750000000000e+01 1.4332499999999999e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 0.0000000000000000e+00 2.8664999999999998e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 2.8664999999999998e+00 2.8664999999999998e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 5.7329999999999997e+00 2.8664999999999998e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 8.5994999999999990e+00 2.8664999999999998e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 1.4332499999999999e+00 4.2997499999999995e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 4.2997499999999995e+00 4.2997499999999995e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 7.1662499999999998e+00 4.2997499999999995e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 1.0032750000000000e+01 4.2997499999999995e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 0.0000000000000000e+00 5.7329999999999997e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 1.4332499999999999e+00 7.1662499999999998e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 2.8664999999999998e+00 5.7329999999999997e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 4.2997499999999995e+00 7.1662499999999998e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 5.7329999999999997e+00 5.7329999999999997e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 7.1662499999999998e+00 7.1662499999999998e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 8.5994999999999990e+00 5.7329999999999997e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 1.0032750000000000e+01 7.1662499999999998e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 0.0000000000000000e+00 8.5994999999999990e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 2.8664999999999998e+00 8.5994999999999990e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 5.7329999999999997e+00 8.5994999999999990e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 8.5994999999999990e+00 8.5994999999999990e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 1.4332499999999999e+00 1.0032750000000000e+01 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 4.2997499999999995e+00 1.0032750000000000e+01 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 7.1662499999999998e+00 1.0032750000000000e+01 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2.2000000000000002e+00 1.0032750000000000e+01 1.0032750000000000e+01 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+1 2.2000000000000002e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+2 2.2000000000000002e+00 1.4332499999999999e+00 1.4332499999999999e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+3 2.2000000000000002e+00 2.8664999999999998e+00 0.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+4 2.2000000000000002e+00 4.2997499999999995e+00 1.4332499999999999e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+5 2.2000000000000002e+00 5.7329999999999997e+00 0.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+6 2.2000000000000002e+00 7.1662499999999998e+00 1.4332499999999999e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+7 2.2000000000000002e+00 8.5994999999999990e+00 0.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+8 2.2000000000000002e+00 1.0032750000000000e+01 1.4332499999999999e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+9 2.2000000000000002e+00 0.0000000000000000e+00 2.8664999999999998e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+10 2.2000000000000002e+00 2.8664999999999998e+00 2.8664999999999998e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+11 2.2000000000000002e+00 5.7329999999999997e+00 2.8664999999999998e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+12 2.2000000000000002e+00 8.5994999999999990e+00 2.8664999999999998e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+13 2.2000000000000002e+00 1.4332499999999999e+00 4.2997499999999995e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+14 2.2000000000000002e+00 4.2997499999999995e+00 4.2997499999999995e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+15 2.2000000000000002e+00 7.1662499999999998e+00 4.2997499999999995e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+16 2.2000000000000002e+00 1.0032750000000000e+01 4.2997499999999995e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+17 2.2000000000000002e+00 0.0000000000000000e+00 5.7329999999999997e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+18 2.2000000000000002e+00 1.4332499999999999e+00 7.1662499999999998e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+19 2.2000000000000002e+00 2.8664999999999998e+00 5.7329999999999997e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+20 2.2000000000000002e+00 4.2997499999999995e+00 7.1662499999999998e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+21 2.2000000000000002e+00 5.7329999999999997e+00 5.7329999999999997e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+22 2.2000000000000002e+00 7.1662499999999998e+00 7.1662499999999998e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+23 2.2000000000000002e+00 8.5994999999999990e+00 5.7329999999999997e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+24 2.2000000000000002e+00 1.0032750000000000e+01 7.1662499999999998e+00 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+25 2.2000000000000002e+00 0.0000000000000000e+00 8.5994999999999990e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+26 2.2000000000000002e+00 2.8664999999999998e+00 8.5994999999999990e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+27 2.2000000000000002e+00 5.7329999999999997e+00 8.5994999999999990e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+28 2.2000000000000002e+00 8.5994999999999990e+00 8.5994999999999990e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+29 2.2000000000000002e+00 1.4332499999999999e+00 1.0032750000000000e+01 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+30 2.2000000000000002e+00 4.2997499999999995e+00 1.0032750000000000e+01 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+31 2.2000000000000002e+00 7.1662499999999998e+00 1.0032750000000000e+01 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+32 2.2000000000000002e+00 1.0032750000000000e+01 1.0032750000000000e+01 1.4332499999999999e+00 -1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
 
 
 
diff --git a/examples/SPIN/gneb_bfo/in.neb.hop1 b/examples/SPIN/gneb_bfo/in.neb.hop1
index 4e203afd82..6697330faa 100644
--- a/examples/SPIN/gneb_bfo/in.neb.hop1
+++ b/examples/SPIN/gneb_bfo/in.neb.hop1
@@ -1,4 +1,5 @@
 # 2d NEB surface simulation, hop from surface to become adatom
+print "Test 1"
 
 dimension	2
 boundary	p s p
diff --git a/examples/SPIN/gneb_bfo/in.neb.spin_iron b/examples/SPIN/gneb_bfo/in.neb.spin_iron
index 6b99cd4639..2cf3a8653e 100644
--- a/examples/SPIN/gneb_bfo/in.neb.spin_iron
+++ b/examples/SPIN/gneb_bfo/in.neb.spin_iron
@@ -1,12 +1,12 @@
 # bcc iron in a 3d periodic box
+print "Test 1"
 
-clear 
 units 		metal
-atom_style 	spin
-
 dimension 	3
 boundary 	p p f
 
+atom_style 	spin
+
 # necessary for the serial algorithm (sametag)
 atom_modify 	map array 
 
@@ -20,7 +20,7 @@ read_data        ../examples/SPIN/gneb_bfo/initial.iron_spin
 # 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
+#set 		group all spin 2.2 -1.0 0.0 0.0
 
 pair_style 	spin/exchange 3.5
 pair_coeff 	* * exchange 3.4 0.02726 0.2171 1.841
@@ -28,73 +28,13 @@ pair_coeff 	* * 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 zeeman 0.0 0.0 0.0 1.0
-fix 		2 all langevin/spin 0.0 0.0 21
-fix		2 nebatoms neb_spin 1.0 
+fix 		1 all precession/spin zeeman 0.001 0.0 0.0 1.0 anisotropy 0.01 1.0 0.0 0.0 
+#fix 		2 all langevin/spin 0.0 0.0 21
+fix		2 all neb/spin 1.0 
 #parallel ideal
 
-min_style	quickmin
-neb		0.0 0.1 1000 1000 100 final ../examples/SPIN/gneb_bfo/final.iron_spin
+timestep	0.0001
 
-
-#############################################################
-
-# 2d NEB surface simulation, hop from surface to become adatom
-
-
-variable	u uloop 20
-
-# create geometry with flat surface
-
-lattice		hex 0.9
-region		box block 0 20 0 10 -0.25 0.25
-
-#create_box	3 box
-#create_atoms	1 box
-#mass		* 1.0
-#write_data      initial.hop1
-
-read_data        ../examples/SPIN/gneb_bfo/initial.hop1
-
-# LJ potentials
-
-pair_style	lj/cut 2.5
-pair_coeff	* * 1.0 1.0 2.5
-pair_modify	shift yes
-
-# initial minimization to relax surface
-
-minimize	1.0e-6 1.0e-4 1000 10000
-reset_timestep	0
-
-# define groups
-
-region	        1 block INF INF INF 1.25 INF INF
-group		lower region 1
-group		mobile subtract all lower
-set		group lower type 2
-
-timestep	0.05
-
-# group of NEB atoms - either block or single atom ID 412
-
-region		surround block 10 18 17 20 0 0 units box
-group		nebatoms region surround
-#group		nebatoms id 412
-set		group nebatoms type 3
-group		nonneb subtract all nebatoms
-
-fix		1 lower setforce 0.0 0.0 0.0
-fix		2 nebatoms neb 1.0 parallel ideal
-fix		3 all enforce2d
-
-thermo		100
-
-#dump		1 nebatoms atom 10 dump.neb.$u
-#dump		2 nonneb atom 10 dump.nonneb.$u
-
-# run NEB for 2000 steps or to force tolerance
-
-min_style	quickmin
-
-neb		0.0 0.1 1000 1000 100 final ../examples/SPIN/gneb_bfo/final.hop1
+#min_style	quickmin
+neb/spin		0.0 0.1 1 1 1 final ../examples/SPIN/gneb_bfo/final.iron_spin
+#neb/spin		0.0 0.1 1000 1000 100 final ../examples/SPIN/gneb_bfo/final.iron_spin
diff --git a/examples/SPIN/gneb_bfo/in.spin.iron b/examples/SPIN/gneb_bfo/in.spin.iron
index 6d291c633b..0c84845f5f 100644
--- a/examples/SPIN/gneb_bfo/in.spin.iron
+++ b/examples/SPIN/gneb_bfo/in.spin.iron
@@ -19,7 +19,7 @@ create_atoms 	1 box
 
 mass		1 55.845
 
-set 		group all spin 2.2 -1.0 0.0 0.0
+set 		group all spin 2.2 1.0 0.0 0.0
 #velocity 	all create 100 4928459 rot yes dist gaussian
 
 pair_style 	spin/exchange 3.5
@@ -53,4 +53,4 @@ 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 		0
-write_data	final.iron_spin
+write_data	initial.iron_spin
diff --git a/src/REPLICA/fix_neb_spin.cpp b/src/REPLICA/fix_neb_spin.cpp
index a96af45267..150c37a03e 100644
--- a/src/REPLICA/fix_neb_spin.cpp
+++ b/src/REPLICA/fix_neb_spin.cpp
@@ -11,11 +11,6 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-/* ----------------------------------------------------------------------
-   Contributing author for: Emile Maras (CEA, France)
-     new options for inter-replica forces, first/last replica treatment
-------------------------------------------------------------------------- */
-
 #include <mpi.h>
 #include <cmath>
 #include <cstdlib>
@@ -392,6 +387,7 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
   nlen = 0.0;
   double tlen = 0.0;
   double gradnextlen = 0.0;
+  double delndots, delpdots;
 
   dotgrad = gradlen = dotpath = dottangrad = 0.0;
 
@@ -403,12 +399,20 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
 
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
+	
 	// tangent vector
 	delspxp = sp[i][0] - spprev[i][0];
 	delspyp = sp[i][1] - spprev[i][1];
 	delspzp = sp[i][2] - spprev[i][2];
-        domain->minimum_image(delspxp,delspyp,delspzp); // check what it does
-	delsqp = delspxp*delspxp+delspyp*delspyp+delspzp*delspzp;
+        
+	// project delp vector on tangent space
+	delpdots = delspxp*sp[i][0]+delspyp*sp[i][1]+delspzp*sp[i][2];
+	delspxp -= delpdots*sp[i][0];
+	delspyp -= delpdots*sp[i][1];
+	delspzp -= delpdots*sp[i][2];
+        
+        // adjust distance if pbc	
+	//domain->minimum_image(delspxp,delspyp,delspzp);
 
 	// calc. geodesic length
 	spi[0]=sp[i][0];
@@ -419,7 +423,7 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
 	spj[2]=spprev[i][2];
 	templen = geodesic_distance(spi, spj);
 	plen += templen*templen;
-	dottangrad += delxp*fm[i][0]+ delyp*fm[i][1]+delzp*fm[i][2];
+	dottangrad += delspxp*fm[i][0]+ delspyp*fm[i][1]+delspzp*fm[i][2];
         gradlen += fm[i][0]*fm[i][0] + fm[i][1]*fm[i][1] + fm[i][2]*fm[i][2];
 
 	//plen += delxp*delxp + delyp*delyp + delzp*delzp;
@@ -430,14 +434,14 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
 	// (unless FreeEnd option)
         if (FreeEndFinal||FreeEndFinalWithRespToEIni) {
 	  error->all(FLERR,"Free End option not yet active");
-	  tangent[i][0]=delspxp;
-	  tangent[i][1]=delspyp;
-	  tangent[i][2]=delspzp;
-	  // if needed, tlen has to be modified
-          tlen += tangent[i][0]*tangent[i][0] +
-            tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2];
-          dot += fm[i][0]*tangent[i][0] + fm[i][1]*tangent[i][1] +
-            fm[i][2]*tangent[i][2];
+	  //tangent[i][0]=delspxp;
+	  //tangent[i][1]=delspyp;
+	  //tangent[i][2]=delspzp;
+	  //// if needed, tlen has to be modified
+          //tlen += tangent[i][0]*tangent[i][0] +
+          //  tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2];
+          //dot += fm[i][0]*tangent[i][0] + fm[i][1]*tangent[i][1] +
+          //  fm[i][2]*tangent[i][2];
 	}
        
        
@@ -469,9 +473,16 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
 	delspxn = spnext[i][0]- sp[i][0];
 	delspyn = spnext[i][1]- sp[i][1];
 	delspzn = spnext[i][2]- sp[i][2];
-        domain->minimum_image(delspxn,delspyn,delspzn); // check what it does
-	delsqn = delspxn*delspxn+delspyn*delspyn+delspzn*delspzn;
 
+        // project deln vector on tangent space
+	delndots = delspxn*sp[i][0]+delspyn*sp[i][1]+delspzn*sp[i][2];
+	delspxn -= delndots*sp[i][0];
+	delspyn -= delndots*sp[i][1];
+	delspzn -= delndots*sp[i][2];
+
+        // adjust del. if pbc	
+        //domain->minimum_image(delspxn,delspyn,delspzn);
+	
 	// calc. geodesic length
 	spi[0]=sp[i][0];
 	spi[1]=sp[i][1];
@@ -485,14 +496,14 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
         gradlen += fm[i][0]*fm[i][0] + fm[i][1]*fm[i][1] + fm[i][2]*fm[i][2];
         if (FreeEndIni) {
 	  error->all(FLERR,"Free End option not yet active");
-          tangent[i][0]=delxn;
-          tangent[i][1]=delyn;
-          tangent[i][2]=delzn;
-	  // if needed, tlen has to be modified
-          tlen += tangent[i][0]*tangent[i][0] +
-            tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2];
-          dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] +
-            f[i][2]*tangent[i][2];
+          //tangent[i][0]=delxn;
+          //tangent[i][1]=delyn;
+          //tangent[i][2]=delzn;
+	  //// if needed, tlen has to be modified
+          //tlen += tangent[i][0]*tangent[i][0] +
+          //  tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2];
+          //dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] +
+          //  f[i][2]*tangent[i][2];
         }
 
 	//delxn = xnext[i][0] - x[i][0];
@@ -528,10 +539,20 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
 
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
-        delspxp = spx[i][0] - spxprev[i][0];
-        delspyp = spx[i][1] - spxprev[i][1];
-        delspzp = spx[i][2] - spxprev[i][2];
-        domain->minimum_image(delspxp,delspyp,delspzp);
+        
+	// calc. delp vector
+	delspxp = sp[i][0] - spprev[i][0];
+        delspyp = sp[i][1] - spprev[i][1];
+        delspzp = sp[i][2] - spprev[i][2];
+        
+	// project delp vector on tangent space
+	delndots = delspxp*sp[i][0]+delspyp*sp[i][1]+delspzp*sp[i][2];
+	delspxp -= delpdots*sp[i][0];
+	delspyp -= delpdots*sp[i][1];
+	delspzp -= delpdots*sp[i][2];
+       
+        // adjust distance if pbc	
+	//domain->minimum_image(delspxp,delspyp,delspzp);
 
 	// calc. geodesic length
 	spi[0]=sp[i][0];
@@ -543,10 +564,19 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
 	templen = geodesic_distance(spi, spj);
 	plen += templen*templen;
 
-        delspxn = spxnext[i][0] - spx[i][0];
-        delspyn = spxnext[i][1] - spx[i][1];
-        delspzn = spxnext[i][2] - spx[i][2];
-        domain->minimum_image(delspxn,delspyn,delspzn);
+	// calc. deln vector
+        delspxn = spnext[i][0] - sp[i][0];
+        delspyn = spnext[i][1] - sp[i][1];
+        delspzn = spnext[i][2] - sp[i][2];
+        
+	// project deln vector on tangent space
+	delndots = delspxn*sp[i][0]+delspyn*sp[i][1]+delspzn*sp[i][2];
+	delspxn -= delndots*sp[i][0];
+	delspyn -= delndots*sp[i][1];
+	delspzn -= delndots*sp[i][2];
+        
+        // adjust distance if pbc	
+	//domain->minimum_image(delspxn,delspyn,delspzn);
 
         if (vnext > veng && veng > vprev) {
           tangent[i][0] = delspxn;
@@ -593,63 +623,11 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
         dotgrad += fm[i][0]*fnext[i][0] + fm[i][1]*fnext[i][1] +
           fm[i][2]*fnext[i][2];
 
-	// is this a perpandicular spring force?
+        // no Perpendicular nudging force option active yet
+        // see fix_neb for example
         if (kspringPerp != 0.0) 
-	  error->all(FLERR,"Perpendicular nudging force not yet active");
-        //springF[i][0] = kspringPerp*(delxn-delxp);
-        //springF[i][1] = kspringPerp*(delyn-delyp);
-        //springF[i][2] = kspringPerp*(delzn-delzp);
+	  error->all(FLERR,"NEB_spin Perpendicular nudging force not yet active");
 
-        //delxp = x[i][0] - xprev[i][0];
-        //delyp = x[i][1] - xprev[i][1];
-        //delzp = x[i][2] - xprev[i][2];
-        //domain->minimum_image(delxp,delyp,delzp);
-        //plen += delxp*delxp + delyp*delyp + delzp*delzp;
-
-        //delxn = xnext[i][0] - x[i][0];
-        //delyn = xnext[i][1] - x[i][1];
-        //delzn = xnext[i][2] - x[i][2];
-        //domain->minimum_image(delxn,delyn,delzn);
-
-        //if (vnext > veng && veng > vprev) {
-        //  tangent[i][0] = delxn;
-        //  tangent[i][1] = delyn;
-        //  tangent[i][2] = delzn;
-        //} else if (vnext < veng && veng < vprev) {
-        //  tangent[i][0] = delxp;
-        //  tangent[i][1] = delyp;
-        //  tangent[i][2] = delzp;
-        //} else {
-        //  if (vnext > vprev) {
-        //    tangent[i][0] = vmax*delxn + vmin*delxp;
-        //    tangent[i][1] = vmax*delyn + vmin*delyp;
-        //    tangent[i][2] = vmax*delzn + vmin*delzp;
-        //  } else if (vnext < vprev) {
-        //    tangent[i][0] = vmin*delxn + vmax*delxp;
-        //    tangent[i][1] = vmin*delyn + vmax*delyp;
-        //    tangent[i][2] = vmin*delzn + vmax*delzp;
-        //  } else { // vnext == vprev, e.g. for potentials that do not compute an energy
-        //    tangent[i][0] = delxn + delxp;
-        //    tangent[i][1] = delyn + delyp;
-        //    tangent[i][2] = delzn + delzp;
-        //  }
-        //}
-
-        //nlen += delxn*delxn + delyn*delyn + delzn*delzn;
-        //tlen += tangent[i][0]*tangent[i][0] +
-        //  tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2];
-        //gradlen += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2];
-        //dotpath += delxp*delxn + delyp*delyn + delzp*delzn;
-        //dottangrad += tangent[i][0]*f[i][0] +
-        //  tangent[i][1]*f[i][1] + tangent[i][2]*f[i][2];
-        //gradnextlen += fnext[i][0]*fnext[i][0] +
-        //  fnext[i][1]*fnext[i][1] +fnext[i][2] * fnext[i][2];
-        //dotgrad += f[i][0]*fnext[i][0] + f[i][1]*fnext[i][1] +
-        //  f[i][2]*fnext[i][2];
-
-        //springF[i][0] = kspringPerp*(delxn-delxp);
-        //springF[i][1] = kspringPerp*(delyn-delyp);
-        //springF[i][2] = kspringPerp*(delzn-delzp);
       }
   }
 
@@ -675,6 +653,24 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
   dottangrad = bufout[6];
   dotgrad = bufout[7];
 
+
+  // project tangent vector on tangent space
+
+  double buftan[3];
+  double tandots;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      tandots = tangent[i][0]*sp[i][0]+tangent[i][1]*sp[i][1]+
+	tangent[i][2]*sp[i][2];
+      buftan[0] = tangent[i][0]-tandots*sp[i][0];
+      buftan[1] = tangent[i][1]-tandots*sp[i][1];
+      buftan[2] = tangent[i][2]-tandots*sp[i][2];
+      tangent[i][0] = buftan[0];
+      tangent[i][1] = buftan[1];
+      tangent[i][2] = buftan[2];
+    }
+
+
   // normalize tangent vector
 
   if (tlen > 0.0) {
@@ -687,8 +683,6 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
       }
   }
 
-////////////////////////////////////////////////////////
-
   // first or last replica has no change to forces, just return
 
   if (ireplica > 0 && ireplica < nreplica-1)
@@ -700,93 +694,31 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
   if (ireplica < nreplica-1)
     dotgrad = dotgrad /(gradlen*gradnextlen);
 
+  // no Free End option active yet
+  // see fix_neb for example
   if (FreeEndIni && ireplica == 0) {
-    if (tlen > 0.0) {
-      double dotall;
-      MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
-      dot=dotall/tlen;
-
-      if (dot<0) prefactor = -dot - kspringIni*(veng-EIniIni);
-      else prefactor = -dot + kspringIni*(veng-EIniIni);
-
-      for (int i = 0; i < nlocal; i++)
-        if (mask[i] & groupbit) {
-          f[i][0] += prefactor *tangent[i][0];
-          f[i][1] += prefactor *tangent[i][1];
-          f[i][2] += prefactor *tangent[i][2];
-        }
-    }
+    error->all(FLERR,"NEB_spin Free End option not yet active");
   }
 
+  // no Free End option active yet
+  // see fix_neb for example
   if (FreeEndFinal && ireplica == nreplica -1) {
-    if (tlen > 0.0) {
-      double dotall;
-      MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
-      dot=dotall/tlen;
-
-      if (veng<EFinalIni) {
-        if (dot<0) prefactor = -dot - kspringFinal*(veng-EFinalIni);
-        else prefactor = -dot + kspringFinal*(veng-EFinalIni);
-      }
-      for (int i = 0; i < nlocal; i++)
-        if (mask[i] & groupbit) {
-          f[i][0] += prefactor *tangent[i][0];
-          f[i][1] += prefactor *tangent[i][1];
-          f[i][2] += prefactor *tangent[i][2];
-        }
-    }
+    error->all(FLERR,"NEB_spin Free End option not yet active");
   }
 
+  // no Free End option active yet
+  // see fix_neb for example
   if (FreeEndFinalWithRespToEIni&&ireplica == nreplica -1) {
-    if (tlen > 0.0) {
-      double dotall;
-      MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
-      dot=dotall/tlen;
-      if (veng<vIni) {
-        if (dot<0) prefactor = -dot - kspringFinal*(veng-vIni);
-        else prefactor = -dot + kspringFinal*(veng-vIni);
-      }
-      for (int i = 0; i < nlocal; i++)
-        if (mask[i] & groupbit) {
-          f[i][0] += prefactor *tangent[i][0];
-          f[i][1] += prefactor *tangent[i][1];
-          f[i][2] += prefactor *tangent[i][2];
-        }
-    }
+    error->all(FLERR,"NEB_spin Free End option not yet active");
   }
 
+  // no NEB_spin long range option 
+  // see fix_neb for example
   double lentot = 0;
   double meanDist,idealPos,lenuntilIm,lenuntilClimber;
   lenuntilClimber=0;
   if (NEBLongRange) {
-    if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) {
-      MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,uworld);
-    } else {
-      if (me == 0)
-        MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,rootworld);
-      MPI_Bcast(nlenall,nreplica,MPI_DOUBLE,0,world);
-    }
-
-    lenuntilIm = 0;
-    for (int i = 0; i < ireplica; i++)
-      lenuntilIm += nlenall[i];
-
-    for (int i = 0; i < nreplica; i++)
-      lentot += nlenall[i];
-
-    meanDist = lentot/(nreplica -1);
-
-    if (rclimber>0) {
-      for (int i = 0; i < rclimber; i++)
-        lenuntilClimber += nlenall[i];
-      double meanDistBeforeClimber = lenuntilClimber/rclimber;
-      double meanDistAfterClimber =
-        (lentot-lenuntilClimber)/(nreplica-rclimber-1);
-      if (ireplica<rclimber)
-        idealPos = ireplica * meanDistBeforeClimber;
-      else
-        idealPos = lenuntilClimber+ (ireplica-rclimber)*meanDistAfterClimber;
-    } else idealPos = ireplica * meanDist;
+    error->all(FLERR,"NEB_spin long range option not yet active");
   }
 
   if (ireplica == 0 || ireplica == nreplica-1) return ;
@@ -800,12 +732,19 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
-      dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] +
-        f[i][2]*tangent[i][2];
-      dotSpringTangent += springF[i][0]*tangent[i][0] +
+      dot += fm[i][0]*tangent[i][0] + fm[i][1]*tangent[i][1] +
+        fm[i][2]*tangent[i][2];
+      // springF defined for perp. spring option
+      // not defined here
+      //dotSpringTangent += springF[i][0]*tangent[i][0] +
         springF[i][1]*tangent[i][1] + springF[i][2]*tangent[i][2];}
+      //dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] +
+      //  f[i][2]*tangent[i][2];
+      //dotSpringTangent += springF[i][0]*tangent[i][0] +
+      //  springF[i][1]*tangent[i][1] + springF[i][2]*tangent[i][2];}
   }
 
+  // gather all dot and dotSpring for this replica (world)
   double dotSpringTangentall;
   MPI_Allreduce(&dotSpringTangent,&dotSpringTangentall,1,
                 MPI_DOUBLE,MPI_SUM,world);
@@ -814,20 +753,26 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
   MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
   dot=dotall;
 
-  if (ireplica == rclimber) prefactor = -2.0*dot;
-  else {
+
+  // implement climbing image here
+
+  if (ireplica == rclimber) {
+    error->all(FLERR,"NEB_spin climber option not yet active");
+    //prefactor = -2.0*dot;
+  } else {
     if (NEBLongRange) {
-      prefactor = -dot - kspring*(lenuntilIm-idealPos)/(2*meanDist);
+      error->all(FLERR,"NEB_spin climber option not yet active");
+      //prefactor = -dot - kspring*(lenuntilIm-idealPos)/(2*meanDist);
     } else if (StandardNEB) {
       prefactor = -dot + kspring*(nlen-plen);
     }
 
     if (FinalAndInterWithRespToEIni&& veng<vIni) {
       for (int i = 0; i < nlocal; i++)
-        if (mask[i] & groupbit) {
-          f[i][0] = 0;
-          f[i][1] = 0;
-          f[i][2] = 0;
+       if (mask[i] & groupbit) {
+          fm[i][0] = 0;
+          fm[i][1] = 0;
+          fm[i][2] = 0;
         }
       prefactor =  kspring*(nlen-plen);
       AngularContr=0;
@@ -836,20 +781,38 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
-      f[i][0] += prefactor*tangent[i][0] +
-        AngularContr*(springF[i][0] - dotSpringTangent*tangent[i][0]);
-      f[i][1] += prefactor*tangent[i][1] +
-        AngularContr*(springF[i][1] - dotSpringTangent*tangent[i][1]);
-      f[i][2] += prefactor*tangent[i][2] +
-        AngularContr*(springF[i][2] - dotSpringTangent*tangent[i][2]);
+      fm[i][0] += prefactor*tangent[i][0];
+      fm[i][1] += prefactor*tangent[i][1];
+      fm[i][2] += prefactor*tangent[i][2];
+      // springF and dotSpringTangent defined for the perp. spring
+      // option, not defined yet for spins
+      //fm[i][0] += prefactor*tangent[i][0] +
+      //  AngularContr*(springF[i][0] - dotSpringTangent*tangent[i][0]);
+      //fm[i][1] += prefactor*tangent[i][1] +
+      //  AngularContr*(springF[i][1] - dotSpringTangent*tangent[i][1]);
+      //fm[i][2] += prefactor*tangent[i][2] +
+      //  AngularContr*(springF[i][2] - dotSpringTangent*tangent[i][2]);
     }
+  
+  // project NEB force on tangent space
+
+  double fdots;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      fdots = fm[i][0]*sp[i][0]+fm[i][1]*sp[i][1]+
+	fm[i][2]*sp[i][2];
+      fm[i][0] -= fdots*sp[i][0];
+      fm[i][1] -= fdots*sp[i][1];
+      fm[i][2] -= fdots*sp[i][2];
+    }
+
 }
 
 /* ----------------------------------------------------------------------
    geodesic distance calculation (Vincenty's formula)
 ------------------------------------------------------------------------- */
 
-double geodesic_distance(double *spi, double *spj)
+double FixNEB_spin::geodesic_distance(double spi[3], double spj[3])
 {
   double dist;
   double crossx,crossy,crossz;
@@ -859,8 +822,7 @@ double geodesic_distance(double *spi, double *spj)
   crossx = spi[1]*spj[2]-spi[2]*spj[1];
   crossy = spi[2]*spj[0]-spi[0]*spj[2];
   crossz = spi[0]*spj[1]-spi[1]*spj[0];
-  crosslen = sqrt(crossspxp*crossspxp + crossspyp*crossspyp +
-      crossz*crossz);
+  crosslen = sqrt(crossx*crossx + crossy*crossy + crossz*crossz);
   dotx = spi[0]*spj[0];
   doty = spi[1]*spj[1];
   dotz = spi[2]*spj[2];
@@ -871,6 +833,149 @@ double geodesic_distance(double *spi, double *spj)
   return dist;
 }
 
+/* ----------------------------------------------------------------------
+   geometric damped advance os spins
+---------------------------------------------------------------------- */
+
+void FixNEB_spin::advance_spins(double dts)
+{
+  //int j=0;
+  //int *sametag = atom->sametag;
+  int nlocal = atom->nlocal;
+  int *mask = atom->mask;
+  double **sp = atom->sp;
+  double **fm = atom->fm;
+  double tdampx,tdampy,tdampz;
+  double msq,scale,fm2,energy,dts2;
+  double alpha;
+  double spi[3],fmi[3];
+  double cp[3],g[3];
+
+  //cp[0] = cp[1] = cp[2] = 0.0;
+  //g[0] = g[1] = g[2] = 0.0;
+  dts2 = dts*dts;		
+
+  // fictitious Gilbert damping of 1
+  alpha = 1.0;
+
+  // loop on all spins on proc.
+
+  if (ireplica != nreplica-1 && ireplica != 0)
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+
+	spi[0] = sp[i][0];
+	spi[1] = sp[i][1];
+	spi[2] = sp[i][2];
+
+	fmi[0] = fm[i][0];
+	fmi[1] = fm[i][1];
+	fmi[2] = fm[i][2];
+
+        // calc. damping torque
+      
+        tdampx = -alpha*(fmi[1]*spi[2] - fmi[2]*spi[1]);
+        tdampy = -alpha*(fmi[2]*spi[0] - fmi[0]*spi[2]);
+        tdampz = -alpha*(fmi[0]*spi[1] - fmi[1]*spi[0]);
+      
+        // apply advance algorithm (geometric, norm preserving)
+        
+        fm2 = (tdampx*tdampx+tdampy*tdampy+tdampz*tdampz);
+        energy = (sp[i][0]*tdampx)+(sp[i][1]*tdampy)+(sp[i][2]*tdampz);
+        
+        cp[0] = tdampy*sp[i][2]-tdampz*sp[i][1];
+        cp[1] = tdampz*sp[i][0]-tdampx*sp[i][2];
+        cp[2] = tdampx*sp[i][1]-tdampy*sp[i][0];
+      
+        g[0] = sp[i][0]+cp[0]*dts;
+        g[1] = sp[i][1]+cp[1]*dts;
+        g[2] = sp[i][2]+cp[2]*dts;
+      			
+        g[0] += (fm[i][0]*energy-0.5*sp[i][0]*fm2)*0.5*dts2;
+        g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts2;
+        g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts2;
+      			
+        g[0] /= (1+0.25*fm2*dts2);
+        g[1] /= (1+0.25*fm2*dts2);
+        g[2] /= (1+0.25*fm2*dts2);
+
+        sp[i][0] = g[0];
+        sp[i][1] = g[1];
+        sp[i][2] = g[2];			
+      
+        // 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;
+      
+        // comm. sp[i] to atoms with same tag (for serial algo)
+      
+        // no need for simplecticity
+        //if (sector_flag == 0) {
+        //  if (sametag[i] >= 0) {
+        //    j = sametag[i];
+        //    while (j >= 0) {
+        //      sp[j][0] = sp[i][0];
+        //      sp[j][1] = sp[i][1];
+        //      sp[j][2] = sp[i][2];
+        //      j = sametag[j];
+        //    }
+        //  }
+        //}
+        //
+
+      }
+}
+
+/* ----------------------------------------------------------------------
+   evaluate max timestep
+---------------------------------------------------------------------- */
+
+double FixNEB_spin::evaluate_dt()
+{
+  double dtmax;
+  double fmsq;
+  double fmaxsqone,fmaxsqloc,fmaxsqall;
+  int nlocal = atom->nlocal;
+  int *mask = atom->mask;
+  double **fm = atom->fm;
+
+  // finding max fm on this proc. 
+  
+  fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      fmsq = fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2];
+      fmaxsqone = MAX(fmaxsqone,fmsq);
+    }
+
+  // finding max fm on this replica
+ 
+  fmaxsqloc = fmaxsqone; 
+  MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world); 
+
+  // finding max fm over all replicas, if necessary
+  // this communicator would be invalid for multiprocess replicas
+
+  if (update->multireplica == 1) {
+    fmaxsqall = fmaxsqloc;
+    MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);
+  }
+
+  if (fmaxsqall < fmaxsqloc)
+    error->all(FLERR,"Incorrect fmaxall calc.");
+
+  // define max timestep
+  // dividing by 10 the inverse of max frequency
+
+  dtmax = MY_2PI/(10.0*sqrt(fmaxsqall));
+
+  return dtmax;
+}
+
 /* ----------------------------------------------------------------------
    send/recv NEB atoms to/from adjacent replicas
    received atoms matching my local atoms are stored in xprev,xnext
@@ -913,19 +1018,25 @@ void FixNEB_spin::inter_replica_comm()
   if (cmode == SINGLE_PROC_DIRECT) {
     if (ireplica > 0)
       MPI_Irecv(xprev[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld,&request);
+      MPI_Irecv(spprev[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld,&request);
     if (ireplica < nreplica-1)
       MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld);
+      MPI_Send(sp[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld);
     if (ireplica > 0) MPI_Wait(&request,MPI_STATUS_IGNORE);
     if (ireplica < nreplica-1)
       MPI_Irecv(xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request);
+      MPI_Irecv(spnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request);
     if (ireplica > 0)
       MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld);
+      MPI_Send(sp[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld);
     if (ireplica < nreplica-1) MPI_Wait(&request,MPI_STATUS_IGNORE);
 
     if (ireplica < nreplica-1)
       MPI_Irecv(fnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request);
+      MPI_Irecv(fmnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request);
     if (ireplica > 0)
       MPI_Send(f[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld);
+      MPI_Send(fm[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld);
     if (ireplica < nreplica-1) MPI_Wait(&request,MPI_STATUS_IGNORE);
 
     return;
@@ -1077,7 +1188,7 @@ void FixNEB_spin::inter_replica_comm()
                 fmsendall[0],counts,displacements,MPI_DOUBLE,0,world);
   } else {
     MPI_Gatherv(NULL,3*m,MPI_DOUBLE,
-                xsendall[0],counts,displacements,MPI_DOUBLE,0,world);
+                spsendall[0],counts,displacements,MPI_DOUBLE,0,world);
     MPI_Gatherv(NULL,3*m,MPI_DOUBLE,
                 fmsendall[0],counts,displacements,MPI_DOUBLE,0,world);
   }
diff --git a/src/REPLICA/fix_neb_spin.h b/src/REPLICA/fix_neb_spin.h
index e3bdd6889d..291341860e 100644
--- a/src/REPLICA/fix_neb_spin.h
+++ b/src/REPLICA/fix_neb_spin.h
@@ -13,8 +13,8 @@
 
 #ifdef FIX_CLASS
 
-FixStyle(neb,FixNEB)
-FixStyle(neb_spin,FixNEB_spin)
+//FixStyle(neb,FixNEB)
+FixStyle(neb/spin,FixNEB_spin)
 
 #else
 
@@ -36,6 +36,8 @@ class FixNEB_spin : public Fix {
   void init();
   void min_setup(int);
   void min_post_force(int);
+  void advance_spins(double);
+  double evaluate_dt();
 
  private:
   int me,nprocs,nprocs_universe;
@@ -79,6 +81,7 @@ class FixNEB_spin : public Fix {
 
   int *counts,*displacements;   // used for MPI_Gather
 
+  double geodesic_distance(double *, double *);
   void inter_replica_comm();
   void reallocate();
 };
diff --git a/src/REPLICA/neb_spin.cpp b/src/REPLICA/neb_spin.cpp
index 7860553532..60d44d1875 100644
--- a/src/REPLICA/neb_spin.cpp
+++ b/src/REPLICA/neb_spin.cpp
@@ -13,7 +13,7 @@
 
 // lmptype.h must be first b/c this file uses MAXBIGINT and includes mpi.h
 // due to OpenMPI bug which sets INT64_MAX via its mpi.h
-//   before lmptype.h can set flags to insure it is done correctly
+// before lmptype.h can set flags to insure it is done correctly
 
 #include "lmptype.h"
 #include <mpi.h>
@@ -21,7 +21,10 @@
 #include <cstdlib>
 #include <cstring>
 //#include "neb.h"
+// test spin
 #include "neb_spin.h"
+#include "compute.h"
+
 #include "universe.h"
 #include "atom.h"
 #include "update.h"
@@ -30,7 +33,10 @@
 #include "min.h"
 #include "modify.h"
 #include "fix.h"
-#include "fix_neb.h"
+//#include "fix_neb.h"
+// test spin
+#include "fix_neb_spin.h"
+
 #include "output.h"
 #include "thermo.h"
 #include "finish.h"
@@ -46,7 +52,7 @@ using namespace MathConst;
 #define MAXLINE 256
 #define CHUNK 1024
 //#define ATTRIBUTE_PERLINE 4
-// 8 attributes: tag number, coords, spin norm, spin dir.
+// 8 attributes: tag, spin norm, position (3), spin direction (3)
 #define ATTRIBUTE_PERLINE 8
 
 /* ---------------------------------------------------------------------- */
@@ -61,7 +67,8 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
          int n2steps_in, int nevery_in, double *buf_init, double *buf_final)
   : Pointers(lmp)
 {
-  double delx,dely,delz;
+  //double delx,dely,delz;
+  double delspx,delspy,delspz;
 
   etol = etol_in;
   ftol = ftol_in;
@@ -82,19 +89,40 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
   double fraction = ireplica/(nreplica-1.0);
 
   double **x = atom->x;
+  // spin quantitites
   double **sp = atom->sp;
   int nlocal = atom->nlocal;
 
-  // modif interp.
   int ii = 0;
+  double spinit[3],spfinal[3];
   for (int i = 0; i < nlocal; i++) {
-    delx = buf_final[ii] - buf_init[ii];
-    dely = buf_final[ii+1] - buf_init[ii+1];
-    delz = buf_final[ii+2] - buf_init[ii+2];
-    domain->minimum_image(delx,dely,delz);
-    x[i][0] = buf_init[ii] + fraction*delx;
-    x[i][1] = buf_init[ii+1] + fraction*dely;
-    x[i][2] = buf_init[ii+2] + fraction*delz;
+
+    spinit[0] = buf_init[ii];
+    spinit[1] = buf_init[ii+1];
+    spinit[2] = buf_init[ii+2];
+    spfinal[0] = buf_final[ii];
+    spfinal[1] = buf_final[ii+1];
+    spfinal[2] = buf_final[ii+2];
+
+    initial_rotation(spinit,spfinal,fraction);
+
+    sp[i][0] = spfinal[0];
+    sp[i][1] = spfinal[1];
+    sp[i][2] = spfinal[2];
+
+    //delx = buf_final[ii] - buf_init[ii];
+    //dely = buf_final[ii+1] - buf_init[ii+1];
+    //delz = buf_final[ii+2] - buf_init[ii+2];
+    
+    // adjust distance if pbc
+    // not implemented yet
+    //domain->minimum_image(delx,dely,delz);
+   
+    // need to define a procedure for circular initialization
+
+    //x[i][0] = buf_init[ii] + fraction*delx;
+    //x[i][1] = buf_init[ii+1] + fraction*dely;
+    //x[i][2] = buf_init[ii+2] + fraction*delz;
     ii += 3;
   }
 }
@@ -114,6 +142,15 @@ NEB_spin::~NEB_spin()
 
 void NEB_spin::command(int narg, char **arg)
 {
+
+  printf("test 1 \n");
+
+  // test 1
+  double **sp1 = atom->sp;
+  printf("test 1 atom: i=%d,%g,%g,%g \n",1,sp1[1][0],sp1[1][1],sp1[1][2]);
+  //error->all(FLERR,"end neb_spin test");
+
+
   if (domain->box_exist == 0)
     error->all(FLERR,"NEB_spin command before simulation box is defined");
 
@@ -141,6 +178,13 @@ void NEB_spin::command(int narg, char **arg)
   uworld = universe->uworld;
   MPI_Comm_rank(world,&me);
 
+  // check metal units and spin atom/style 
+
+  if (!atom->sp_flag)
+    error->all(FLERR,"neb/spin requires atom/spin style");
+  if (strcmp(update->unit_style,"metal") != 0)
+    error->all(FLERR,"neb/spin simulation requires metal unit style");
+
   // error checks
 
   if (nreplica == 1) error->all(FLERR,"Cannot use NEB_spin with a single replica");
@@ -149,6 +193,7 @@ void NEB_spin::command(int narg, char **arg)
 
   // process file-style setting to setup initial configs for all replicas
 
+  // check what options are available
   if (strcmp(arg[5],"final") == 0) {
     if (narg != 7 && narg !=8) error->universe_all(FLERR,"Illegal NEB_spin command");
     infile = arg[6];
@@ -163,6 +208,13 @@ void NEB_spin::command(int narg, char **arg)
 
   verbose=false;
   if (strcmp(arg[narg-1],"verbose") == 0) verbose=true;
+  
+
+  // test 1
+  double **sp = atom->sp;
+  printf("test 2 atom: i=%d,%g,%g,%g \n",1,sp[1][0],sp[1][1],sp[1][2]);
+  error->all(FLERR,"end neb_spin test");
+
   // run the NEB_spin calculation
 
   run();
@@ -181,17 +233,13 @@ void NEB_spin::run()
   else color = 1;
   MPI_Comm_split(uworld,color,0,&roots);
 
+  // search for neb_spin fix, allocate it
   int ineb;
   for (ineb = 0; ineb < modify->nfix; ineb++)
-    if (strcmp(modify->fix[ineb]->style,"neb_spin") == 0) break;
-  if (ineb == modify->nfix) error->all(FLERR,"NEB_spin requires use of fix neb_spin");
-  //int ineb;
-  //for (ineb = 0; ineb < modify->nfix; ineb++)
-  //  if (strcmp(modify->fix[ineb]->style,"neb") == 0) break;
-  //if (ineb == modify->nfix) error->all(FLERR,"NEB_spin requires use of fix neb");
+    if (strcmp(modify->fix[ineb]->style,"neb/spin") == 0) break;
+  if (ineb == modify->nfix) error->all(FLERR,"NEB_spin requires use of fix neb/spin");
 
-  //fneb = (FixNEB_spin *) modify->fix[ineb];
-  fneb_spin = (FixNEB_spin *) modify->fix[ineb];
+  fneb = (FixNEB_spin *) modify->fix[ineb];
   if (verbose) numall =7;
   else  numall = 4;
   memory->create(all,nreplica,numall,"neb:all");
@@ -206,8 +254,10 @@ void NEB_spin::run()
 
   lmp->init();
 
-  if (update->minimize->searchflag)
-    error->all(FLERR,"NEB_spin requires damped dynamics minimizer");
+  // put flag to check gilbert damping procedure is set
+  
+  //if (update->minimize->searchflag)
+  //  error->all(FLERR,"NEB_spin requires damped dynamics minimizer");
 
   // setup regular NEB_spin minimization
   FILE *uscreen = universe->uscreen;
@@ -264,8 +314,19 @@ void NEB_spin::run()
   timer->init();
   timer->barrier_start();
 
+  double dts;
   while (update->minimize->niter < n1steps) {
-    update->minimize->run(nevery);
+    //dts = evaluate_dt();
+    //advance_spins(dts);
+    dts = fneb->evaluate_dt();
+    fneb->advance_spins(dts);
+    
+    // no minimizer for spins
+    //update->minimize->run(nevery);
+    //
+    // evaluate dts
+    // loop on spins, damped advance
+    //
     print_status();
     if (update->minimize->stop_condition) break;
   }
@@ -309,8 +370,7 @@ void NEB_spin::run()
     error->all(FLERR,"Too many timesteps");
 
   update->minimize->init();
-  //fneb->rclimber = top;
-  fneb_spin->rclimber = top;
+  fneb->rclimber = top;
   update->minimize->setup();
 
   if (me_universe == 0) {
@@ -356,7 +416,11 @@ void NEB_spin::run()
   timer->barrier_start();
 
   while (update->minimize->niter < n2steps) {
-    update->minimize->run(nevery);
+    //dts = evaluate_dt();
+    //advance_spins(dts);
+    dts = fneb->evaluate_dt();
+    fneb->advance_spins(dts);
+    //update->minimize->run(nevery);
     print_status();
     if (update->minimize->stop_condition) break;
   }
@@ -373,6 +437,174 @@ void NEB_spin::run()
   update->beginstep = update->endstep = 0;
 }
 
+/* ----------------------------------------------------------------------
+   geodesic distance calculation (Vincenty's formula)
+------------------------------------------------------------------------- */
+
+//double NEB_spin::geodesic_distance2(double spi[3], double spj[3])
+//{
+//  double dist;
+//  double crossx,crossy,crossz;
+//  double dotx,doty,dotz;
+//  double crosslen,dots;
+//
+//  crossx = spi[1]*spj[2]-spi[2]*spj[1];
+//  crossy = spi[2]*spj[0]-spi[0]*spj[2];
+//  crossz = spi[0]*spj[1]-spi[1]*spj[0];
+//  crosslen = sqrt(crossx*crossx + crossy*crossy + crossz*crossz);
+//  dotx = spi[0]*spj[0];
+//  doty = spi[1]*spj[1];
+//  dotz = spi[2]*spj[2];
+//  dots = dotx+doty+dotz;
+//
+//  dist = atan2(crosslen,dots);
+//
+//  return dist;
+//}
+
+/* ----------------------------------------------------------------------
+   evaluate max timestep
+---------------------------------------------------------------------- */
+
+//double NEB_spin::evaluate_dt()
+//{
+//  double dtmax;
+//  double fmsq;
+//  double fmaxsqone,fmaxsqloc,fmaxsqall;
+//  int nlocal = atom->nlocal;
+//  int *mask = atom->mask;
+//  double **fm = atom->fm;
+//
+//  // finding max fm on this proc. 
+//  
+//  fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;
+//  for (int i = 0; i < nlocal; i++)
+//    if (mask[i] & groupbit) {
+//      fmsq = fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2];
+//      fmaxsqone = MAX(fmaxsqone,fmsq);
+//    }
+//
+//  // finding max fm on this replica
+// 
+//  fmaxsqloc = fmaxsqone; 
+//  MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world); 
+//
+//  // finding max fm over all replicas, if necessary
+//  // this communicator would be invalid for multiprocess replicas
+//
+//  if (update->multireplica == 1) {
+//    fmaxsqall = fmaxsqloc;
+//    MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);
+//  }
+//   
+//  if (fmaxsqall < fmaxsqloc)
+//    error->all(FLERR,"Incorrect fmaxall calc.");
+//
+//  // define max timestep
+//  // dividing by 10 the inverse of max frequency
+//
+//  dtmax = MY_2PI/(10.0*sqrt(fmaxsqall));
+//
+//  return dtmax;
+//}
+
+/* ----------------------------------------------------------------------
+   geometric damped advance os spins
+---------------------------------------------------------------------- */
+
+//void NEB_spin::advance_spins(double dts)
+//{
+//  //int j=0;
+//  //int *sametag = atom->sametag;
+//  int nlocal = atom->nlocal;
+//  int *mask = atom->mask;
+//  double **sp = atom->sp;
+//  double **fm = atom->fm;
+//  double tdampx,tdampy,tdampz;
+//  double msq,scale,fm2,energy,dts2;
+//  double alpha;
+//  double spi[3],fmi[3];
+//  double cp[3],g[3];
+//
+//  //cp[0] = cp[1] = cp[2] = 0.0;
+//  //g[0] = g[1] = g[2] = 0.0;
+//  dts2 = dts*dts;		
+//
+//  // fictitious Gilbert damping of 1
+//  alpha = 1.0;
+//
+//  // loop on all spins on proc.
+//
+//  if (ireplica != nreplica-1 && ireplica != 0)
+//    for (int i = 0; i < nlocal; i++)
+//      if (mask[i] & groupbit) {
+//
+//	spi[0] = sp[i][0];
+//	spi[1] = sp[i][1];
+//	spi[2] = sp[i][2];
+//
+//	fmi[0] = fm[i][0];
+//	fmi[1] = fm[i][1];
+//	fmi[2] = fm[i][2];
+//
+//        // calc. damping torque
+//      
+//        tdampx = -alpha*(fmi[1]*spi[2] - fmi[2]*spi[1]);
+//        tdampy = -alpha*(fmi[2]*spi[0] - fmi[0]*spi[2]);
+//        tdampz = -alpha*(fmi[0]*spi[1] - fmi[1]*spi[0]);
+//      
+//        // apply advance algorithm (geometric, norm preserving)
+//        
+//        fm2 = (tdampx*tdampx+tdampy*tdampy+tdampz*tdampz);
+//        energy = (sp[i][0]*tdampx)+(sp[i][1]*tdampy)+(sp[i][2]*tdampz);
+//        
+//        cp[0] = tdampy*sp[i][2]-tdampz*sp[i][1];
+//        cp[1] = tdampz*sp[i][0]-tdampx*sp[i][2];
+//        cp[2] = tdampx*sp[i][1]-tdampy*sp[i][0];
+//      
+//        g[0] = sp[i][0]+cp[0]*dts;
+//        g[1] = sp[i][1]+cp[1]*dts;
+//        g[2] = sp[i][2]+cp[2]*dts;
+//      			
+//        g[0] += (fm[i][0]*energy-0.5*sp[i][0]*fm2)*0.5*dts2;
+//        g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts2;
+//        g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts2;
+//      			
+//        g[0] /= (1+0.25*fm2*dts2);
+//        g[1] /= (1+0.25*fm2*dts2);
+//        g[2] /= (1+0.25*fm2*dts2);
+//      
+//        sp[i][0] = g[0];
+//        sp[i][1] = g[1];
+//        sp[i][2] = g[2];			
+//      
+//        // 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;
+//      
+//        // comm. sp[i] to atoms with same tag (for serial algo)
+//      
+//        // no need for simplecticity
+//        //if (sector_flag == 0) {
+//        //  if (sametag[i] >= 0) {
+//        //    j = sametag[i];
+//        //    while (j >= 0) {
+//        //      sp[j][0] = sp[i][0];
+//        //      sp[j][1] = sp[i][1];
+//        //      sp[j][2] = sp[i][2];
+//        //      j = sametag[j];
+//        //    }
+//        //  }
+//        //}
+//        //
+//
+//      }
+//}
+
 /* ----------------------------------------------------------------------
    read initial config atom coords from file
    flag = 0
@@ -395,8 +627,9 @@ void NEB_spin::readfile(char *file, int flag)
   char *eof,*start,*next,*buf;
   char line[MAXLINE];
   double xx,yy,zz,delx,dely,delz;
-  // creating new temp. sp
-  double spx,spy,spz,delspx,delspy,delspz;
+  // spin quantities
+  double musp,spx,spy,spz;
+    //,delx,dely,delz;
 
   if (me_universe == 0 && screen)
     fprintf(screen,"Reading NEB_spin coordinate file(s) ...\n");
@@ -440,10 +673,17 @@ void NEB_spin::readfile(char *file, int flag)
   double fraction = ireplica/(nreplica-1.0);
 
   double **x = atom->x;
-  // spin table
+  // spin quantities
   double **sp = atom->sp;
+  double spinit[3],spfinal[3];
   int nlocal = atom->nlocal;
 
+  // test 1.2
+  //double **sp = atom->sp;
+  printf("test 1.2 atom: i=%d,%g,%g,%g \n",1,sp[1][0],sp[1][1],sp[1][2]);
+  //error->all(FLERR,"end neb_spin test");
+
+
   // loop over chunks of lines read from file
   // two versions of read_lines_from_file() for world vs universe bcast
   // count # of atom coords changed so can check for invalid atom IDs in file
@@ -493,22 +733,59 @@ void NEB_spin::readfile(char *file, int flag)
       m = atom->map(tag);
       if (m >= 0 && m < nlocal) {
         ncount++;
-        xx = atof(values[1]);
-        yy = atof(values[2]);
-        zz = atof(values[3]);
+	musp = atof(values[1]);
+	xx = atof(values[2]);
+        yy = atof(values[3]);
+        zz = atof(values[4]);
+        spx = atof(values[5]);
+        spy = atof(values[6]);
+        spz = atof(values[7]);
+	//xx = atof(values[1]);
+        //yy = atof(values[2]);
+        //zz = atof(values[3]);
 
         if (flag == 0) {
-          delx = xx - x[m][0];
-          dely = yy - x[m][1];
-          delz = zz - x[m][2];
-          domain->minimum_image(delx,dely,delz);
-          x[m][0] += fraction*delx;
-          x[m][1] += fraction*dely;
-          x[m][2] += fraction*delz;
+          
+	  // here, function interp. spin states
+
+	  //spinit[0] = x[m][0];
+	  //spinit[1] = x[m][1];
+	  //spinit[2] = x[m][2];
+	  spinit[0] = sp[m][0];
+	  spinit[1] = sp[m][1];
+	  spinit[2] = sp[m][2];
+	  spfinal[0] = spx;
+	  spfinal[1] = spy;
+	  spfinal[2] = spz;
+          //domain->minimum_image(delx,dely,delz);
+
+	  // test 
+	  printf("spinit: %g %g %g \n",spinit[0],spinit[1],spinit[2]);
+	  printf("spfinal bef: %g %g %g \n",spfinal[0],spfinal[1],spfinal[2]);
+
+	  initial_rotation(spinit,spfinal,fraction);
+	  
+	  // test 
+	  printf("spfinal aft: %g %g %g \n",spfinal[0],spfinal[1],spfinal[2]);
+
+	  sp[m][0] = spfinal[0];
+	  sp[m][1] = spfinal[1];
+	  sp[m][2] = spfinal[2];
+	  sp[m][3] = musp;
+	  //delx = xx - x[m][0];
+          //dely = yy - x[m][1];
+          //delz = zz - x[m][2];
+          //x[m][0] += fraction*delx;
+          //x[m][1] += fraction*dely;
+          //x[m][2] += fraction*delz;
         } else {
-          x[m][0] = xx;
+          sp[m][3] = musp;
+	  x[m][0] = xx;
           x[m][1] = yy;
           x[m][2] = zz;
+	  sp[m][0] = spx;
+	  sp[m][1] = spy; 
+	  sp[m][2] = spz; 
         }
       }
 
@@ -518,6 +795,12 @@ void NEB_spin::readfile(char *file, int flag)
     nread += nchunk;
   }
 
+  // test 1.3
+  //double **sp = atom->sp;
+  printf("test 1.3 atom: i=%d,%g,%g,%g \n",1,sp[1][0],sp[1][1],sp[1][2]);
+  //error->all(FLERR,"end neb_spin test");
+
+
   // check that all atom IDs in file were found by a proc
 
   if (flag == 0) {
@@ -550,6 +833,75 @@ void NEB_spin::readfile(char *file, int flag)
   }
 }
 
+/* ----------------------------------------------------------------------
+   initial configuration of spin sploc using Rodrigues' formula
+   interpolates between initial (spi) and final (stored in sploc) 
+------------------------------------------------------------------------- */
+
+void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
+{
+  
+  
+  // implementing initial rotation using atan2
+  
+  //atan2(crosslen,dots);
+  
+  
+  
+  double theta,spdot;
+  double inormdot,ispinorm;
+  double kix,kiy,kiz;
+  double kinorm, ikinorm;
+  double crossx,crossy,crossz;
+
+  //printf("inside rot, spi %g, spf %g \n",spi[0],sploc[0]);
+
+  spdot = spi[0]*sploc[0]+spi[1]*sploc[1]+spi[2]*sploc[2];
+  theta = fraction*acos(spdot);
+ 
+  printf("inside rot, theta %g \n",theta);
+
+  kix = spi[1]*sploc[2]-spi[2]*sploc[1];
+  kiy = spi[2]*sploc[0]-spi[0]*sploc[2];
+  kiz = spi[0]*sploc[1]-spi[1]*sploc[0];
+  
+  //printf("inside rot1.1, ki %g %g %g \n",kix,kiy,kiz);
+
+  inormdot = 1.0/sqrt(spdot);
+  kinorm = kix*kix+kiy*kiy+kiz*kiz;
+  if (kinorm == 0.0) {
+    kix = 0.0; 
+    kiy = 0.0; 
+    kiz = 0.0;
+  } else {
+    ikinorm = 1.0/kinorm; 
+    kix *= ikinorm;
+    kiy *= ikinorm;
+    kiz *= ikinorm;
+  }
+ 
+  //printf("inside rot1.2, kin %g %g %g \n",kix,kiy,kiz);
+
+  crossx = kiy*spi[2]-kiz*spi[1];
+  crossy = kiz*spi[0]-kix*spi[2];
+  crossz = kix*spi[1]-kiy*spi[0];
+  
+  //printf("inside rot1.3, cross %g %g %g \n",crossx,crossy,crossz);
+
+  sploc[0] = spi[0]*cos(theta)+crossx*sin(theta);
+  sploc[1] = spi[1]*cos(theta)+crossy*sin(theta);
+  sploc[2] = spi[2]*cos(theta)+crossz*sin(theta);
+  
+  //printf("inside rot2, spf %g %g %g \n",sploc[0],sploc[1],sploc[2]);
+
+  ispinorm = 1.0/sqrt(sploc[0]*sploc[0]+sploc[1]*sploc[1]+sploc[2]*sploc[2]);
+
+  sploc[0] *= ispinorm;
+  sploc[1] *= ispinorm;
+  sploc[2] *= ispinorm;
+  printf("inside rot2, spf %g %g %g \n",sploc[0],sploc[1],sploc[2]);
+}
+
 /* ----------------------------------------------------------------------
    universe proc 0 opens NEB_spin data file
    test if gzipped
@@ -606,22 +958,15 @@ void NEB_spin::print_status()
   }
 
   double one[7];
-  //one[0] = fneb->veng;
-  //one[1] = fneb->plen;
-  //one[2] = fneb->nlen;
-  //one[3] = fneb->gradlen;
-  one[0] = fneb_neb->veng;
-  one[1] = fneb_neb->plen;
-  one[2] = fneb_neb->nlen;
-  one[3] = fneb_neb->gradlen;
+  one[0] = fneb->veng;
+  one[1] = fneb->plen;
+  one[2] = fneb->nlen;
+  one[3] = fneb->gradlen;
 
   if (verbose) {
-    //one[4] = fneb->dotpath;
-    //one[5] = fneb->dottangrad;
-    //one[6] = fneb->dotgrad;
-    one[4] = fneb_spin->dotpath;
-    one[5] = fneb_spin->dottangrad;
-    one[6] = fneb_spin->dotgrad;
+    one[4] = fneb->dotpath;
+    one[5] = fneb->dottangrad;
+    one[6] = fneb->dotgrad;
   }
 
   if (output->thermo->normflag) one[0] /= atom->natoms;
diff --git a/src/REPLICA/neb_spin.h b/src/REPLICA/neb_spin.h
index 3fa19460fc..6541658fd7 100644
--- a/src/REPLICA/neb_spin.h
+++ b/src/REPLICA/neb_spin.h
@@ -13,7 +13,8 @@
 
 #ifdef COMMAND_CLASS
 
-CommandStyle(neb_spin,NEB_SPIN)
+CommandStyle(neb/spin,NEB_spin)
+//CommandStyle(neb,NEB_spin)
 
 #else
 
@@ -56,7 +57,11 @@ class NEB_spin : protected Pointers {
   double *freplica;            // force on an image
   double *fmaxatomInRepl;      // force on an image
 
+  //double geodesic_distance2(double *, double *);
+  //double evaluate_dt();
+  //void advance_spins(double); 
   void readfile(char *, int);
+  void initial_rotation(double *, double *, double);
   void open(char *);
   void print_status();
 };
diff --git a/src/SPIN/pair_spin_exchange.cpp b/src/SPIN/pair_spin_exchange.cpp
index 8cd9d33abd..ec21fe8838 100644
--- a/src/SPIN/pair_spin_exchange.cpp
+++ b/src/SPIN/pair_spin_exchange.cpp
@@ -153,11 +153,12 @@ void PairSpinExchange::init_style()
   neighbor->requests[irequest]->half = 0;
   neighbor->requests[irequest]->full = 1;
 
-  // checking if nve/spin is a listed fix
+  // checking if nve/spin or neb/spin are a listed fix
 
   int ifix = 0;
   while (ifix < modify->nfix) {
     if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
+    if (strcmp(modify->fix[ifix]->style,"neb/spin") == 0) break;
     ifix++;
   }
   if (ifix == modify->nfix)