diff --git a/examples/SPIN/pppm_spin/Co_PurjaPun_2012.eam.alloy b/examples/SPIN/pppm_spin/Co_PurjaPun_2012.eam.alloy
new file mode 120000
index 0000000000..6a47c9eebe
--- /dev/null
+++ b/examples/SPIN/pppm_spin/Co_PurjaPun_2012.eam.alloy
@@ -0,0 +1 @@
+../cobalt_fcc/Co_PurjaPun_2012.eam.alloy
\ No newline at end of file
diff --git a/examples/SPIN/pppm_spin/exchange_fit_hcp_co/exchange_fit.py b/examples/SPIN/pppm_spin/exchange_fit_hcp_co/exchange_fit.py
new file mode 100644
index 0000000000..fa7dba417e
--- /dev/null
+++ b/examples/SPIN/pppm_spin/exchange_fit_hcp_co/exchange_fit.py
@@ -0,0 +1,32 @@
+#Program fitting the exchange interaction
+#Model curve: Bethe-Slater function
+import numpy as np, pylab, tkinter
+import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
+from decimal import *
+
+print("Loop begin")
+
+#Definition of the Bethe-Slater function
+def func(x,a,b,c):
+    return 4*a*((x/c)**2)*(1-b*(x/c)**2)*np.exp(-(x/c)**2)
+
+#Exchange coeff table (data to fit)
+rdata, Jdata = np.loadtxt('exchange_hcp_co.dat', usecols=(0,1), unpack=True)
+plt.plot(rdata, Jdata, 'b-', label='data')
+
+#Perform the fit
+popt, pcov = curve_fit(func, rdata, Jdata, bounds=(0, [500.,5.,5.]))
+plt.plot(rdata, func(rdata, *popt), 'r--', label='fit')
+
+#Print the fitted params
+print("Parameters: a={:.10} (in meV), b={:.10} (adim), c={:.10} (in Ang)".format(*popt))
+
+#Ploting the result
+plt.xlabel('r_ij')
+pylab.xlim([0,6.5])
+plt.ylabel('J_ij')
+plt.legend()
+plt.show()
+
+print("Loop end")
diff --git a/examples/SPIN/pppm_spin/exchange_fit_hcp_co/exchange_hcp_co.dat b/examples/SPIN/pppm_spin/exchange_fit_hcp_co/exchange_hcp_co.dat
new file mode 100644
index 0000000000..0968fa3edb
--- /dev/null
+++ b/examples/SPIN/pppm_spin/exchange_fit_hcp_co/exchange_hcp_co.dat
@@ -0,0 +1,9 @@
+2.25569176882662 73.37931034482759
+2.3817863397548162 47.99999999999999
+2.4518388791593697 34.39080459770115
+2.507880910683012 31.816091954022987
+2.5359019264448337 28.137931034482747
+2.5779334500875657 25.011494252873554
+2.6339754816112086 19.126436781609186
+2.760070052539404 13.241379310344826
+3.5446584938704033 6.068965517241367
diff --git a/examples/SPIN/pppm_spin/in.dipole.pppm_dipole b/examples/SPIN/pppm_spin/in.dipole.pppm_dipole
new file mode 100644
index 0000000000..804ddad1e2
--- /dev/null
+++ b/examples/SPIN/pppm_spin/in.dipole.pppm_dipole
@@ -0,0 +1,55 @@
+# 3d Lennard-Jones melt
+
+units		lj
+#atom_style	charge
+atom_style      hybrid sphere dipole
+processors      * 1 1
+
+lattice		fcc 0.8442
+#region		box block 0 10 0 10 0 10
+region		box block 0 5 0 5 0 5
+create_box	3 box
+create_atoms	1 box
+mass		* 1.0
+
+region		long block 3 6 0 10 0 10
+set             region long type 2
+set             group all dipole/random 98934 0.75
+#set type 1:2    charge 0.0
+
+velocity	all create 1.0 87287
+
+#pair_style	lj/long/coul/long long off 2.5
+#pair_coeff	* * 1.0 1.0 2.5
+#pair_coeff      * 2 1.0 1.0 5.0
+pair_style      lj/cut/dipole/long 3.0
+pair_coeff      * * 1.0 1.0
+
+#kspace_style    pppm/disp 1.0e-4
+kspace_style    pppm/dipole 1.0e-4
+kspace_modify   gewald/disp 0.1
+
+neighbor	0.3 bin
+neigh_modify	every 2 delay 4 check yes
+
+group		fast type 1
+group		slow type 2
+fix		0 all balance 20 1.0 shift x 5 1.0 &
+                weight group 2 fast 1.0 slow 2.0 weight time 0.66
+
+fix		1 all nve
+
+#dump		id all atom 50 dump.melt
+
+#dump		2 all image 25 image.*.jpg type type &
+#		axes yes 0.8 0.02 view 60 -30
+#dump_modify	2 pad 3
+
+#dump		3 all movie 25 movie.mpg type type &
+#		axes yes 0.8 0.02 view 60 -30
+#dump_modify	3 pad 3
+
+#thermo		50
+thermo		1
+#run		500
+run		5
diff --git a/examples/SPIN/pppm_spin/in.spin.pppm_spin b/examples/SPIN/pppm_spin/in.spin.pppm_spin
new file mode 100644
index 0000000000..87d18f4d16
--- /dev/null
+++ b/examples/SPIN/pppm_spin/in.spin.pppm_spin
@@ -0,0 +1,61 @@
+# hcp cobalt in a 3d periodic box
+
+clear 
+units	 	metal
+atom_style 	spin
+
+dimension 	3
+boundary 	p p p
+
+# necessary for the serial algorithm (sametag)
+atom_modify 	map array 
+
+lattice 	hcp 2.5071
+region 		box block 0.0 20.0 0.0 20.0 0.0 8.0
+create_box 	1 box
+create_atoms 	1 box
+
+# setting mass, mag. moments, and interactions for hcp cobalt
+
+mass		1 58.93
+
+#set 		group all spin/random 31 1.72
+set 		group all spin 1.72 0.0 0.0 1.0 
+velocity 	all create 100 4928459 rot yes dist gaussian
+
+pair_style 	hybrid/overlay eam/alloy spin/exchange 4.0 spin/long 8.0
+#pair_style 	hybrid/overlay eam/alloy spin/exchange 4.0
+pair_coeff 	* * eam/alloy ../examples/SPIN/pppm_spin/Co_PurjaPun_2012.eam.alloy Co
+pair_coeff 	* * spin/exchange exchange 4.0 0.3593 1.135028015e-05 1.064568567
+pair_coeff 	* * spin/long long 8.0
+
+neighbor 	0.1 bin
+neigh_modify 	every 10 check yes delay 20
+
+kspace_style pppm/dipole/spin 1.0e-4
+
+#fix 		1 all precession/spin zeeman 1.0 0.0 0.0 1.0
+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 		3 all nve/spin lattice yes
+
+timestep	0.0001
+
+
+compute 	out_mag    all compute/spin
+compute 	out_pe     all pe
+compute 	out_ke     all ke
+compute 	out_temp   all temp
+
+variable 	magz      equal c_out_mag[3]
+variable 	magnorm   equal c_out_mag[4]
+variable 	emag      equal c_out_mag[5]
+variable 	tmag      equal c_out_mag[6]
+
+thermo_style    custom step time v_magnorm v_emag temp etotal
+thermo          10
+
+compute 	outsp all property/atom spx spy spz sp fmx fmy fmz
+dump 		100 all custom 1 dump_cobalt_hcp.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3]
+
+run 		20000
diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp
index 132389b7d6..8eec8e2542 100644
--- a/src/KSPACE/pppm.cpp
+++ b/src/KSPACE/pppm.cpp
@@ -100,6 +100,9 @@ PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg),
   nyhi_in = nylo_in = nyhi_out = nylo_out = 0;
   nzhi_in = nzlo_in = nzhi_out = nzlo_out = 0;
 
+  // test 
+  nlower = nupper = 0;
+
   density_brick = vdx_brick = vdy_brick = vdz_brick = NULL;
   density_fft = NULL;
   u_brick = NULL;
@@ -1428,12 +1431,13 @@ void PPPM::set_grid_local()
   double zprd = prd[2];
   double zprd_slab = zprd*slab_volfactor;
 
-  double dist[3];
+  double dist[3] = {0.0,0.0,0.0};
   double cuthalf = 0.5*neighbor->skin + qdist;
   if (triclinic == 0) dist[0] = dist[1] = dist[2] = cuthalf;
   else kspacebbox(cuthalf,&dist[0]);
 
   int nlo,nhi;
+  nlo = nhi = 0;
 
   nlo = static_cast<int> ((sublo[0]-dist[0]-boxlo[0]) *
                             nx_pppm/xprd + shift) - OFFSET;
diff --git a/src/KSPACE/pppm_dipole.cpp b/src/KSPACE/pppm_dipole.cpp
index a03f5b9980..4d2b594af8 100644
--- a/src/KSPACE/pppm_dipole.cpp
+++ b/src/KSPACE/pppm_dipole.cpp
@@ -116,9 +116,10 @@ void PPPMDipole::init()
   // error check
 
   dipoleflag = atom->mu?1:0;
-  qsum_qsq(0);
+  qsum_qsq(0); // q[i] might not be declared ?
+
   if (dipoleflag && q2)
-    error->all(FLERR,"Cannot (yet) uses charges with Kspace style PPPMDipole");
+    error->all(FLERR,"Cannot (yet) use charges with Kspace style PPPMDipole");
 
   triclinic_check();
 
@@ -168,7 +169,9 @@ void PPPMDipole::init()
   cutoff = *p_cutoff;
 
   // kspace TIP4P not yet supported
+  // qdist = offset only for TIP4P fictitious charge
 
+  qdist = 0.0;
   if (tip4pflag)
     error->all(FLERR,"Cannot yet use TIP4P with PPPMDipole");
 
@@ -206,7 +209,7 @@ void PPPMDipole::init()
                      "beyond nearest neighbor processor");
 
     compute_gf_denom();
-    set_grid_global(mu2);
+    set_grid_global();
     set_grid_local();
     if (overlap_allowed) break;
 
@@ -235,7 +238,7 @@ void PPPMDipole::init()
 
   // calculate the final accuracy
 
-  double estimated_accuracy = final_accuracy_dipole(mu2);
+  double estimated_accuracy = final_accuracy_dipole();
 
   // print stats
 
@@ -607,7 +610,7 @@ void PPPMDipole::allocate()
   // summation coeffs
 
   order_allocated = order;
-  memory->create(gf_b,order,"pppm_dipole:gf_b");
+  if (!gf_b) memory->create(gf_b,order,"pppm_dipole:gf_b");
   memory->create2d_offset(rho1d,3,-order/2,order/2,"pppm_dipole:rho1d");
   memory->create2d_offset(drho1d,3,-order/2,order/2,"pppm_dipole:drho1d");
   memory->create2d_offset(rho_coeff,order,(1-order)/2,order/2,"pppm_dipole:rho_coeff");
@@ -791,7 +794,8 @@ void PPPMDipole::deallocate_peratom()
    used for charge accumulation, FFTs, and electric field interpolation
 ------------------------------------------------------------------------- */
 
-void PPPMDipole::set_grid_global(double dipole2)
+//void PPPMDipole::set_grid_global(double dipole2)
+void PPPMDipole::set_grid_global()
 {
   // use xprd,yprd,zprd
   // adjust z dimension for 2d slab PPPMDipole
@@ -813,14 +817,11 @@ void PPPMDipole::set_grid_global(double dipole2)
   if (!gewaldflag) {
     if (accuracy <= 0.0)
       error->all(FLERR,"KSpace accuracy must be > 0");
-    //if (mu2 == 0.0)
-    if (dipole2 == 0.0)
+    if (mu2 == 0.0)
      error->all(FLERR,"Must use kspace_modify gewald for systems with no dipoles");
     g_ewald = (1.35 - 0.15*log(accuracy))/cutoff;
-    //Try Newton Solver
     double g_ewald_new =
-      find_gewald_dipole(g_ewald,cutoff,natoms,xprd*yprd*zprd,dipole2);
-      //find_gewald_dipole(g_ewald,cutoff,natoms,xprd*yprd*zprd,mu2);
+      find_gewald_dipole(g_ewald,cutoff,natoms,xprd*yprd*zprd,mu2);
     if (g_ewald_new > 0.0) g_ewald = g_ewald_new;
     else error->warning(FLERR,"PPPMDipole dipole Newton solver failed, "
                         "using old method to estimate g_ewald");
@@ -837,6 +838,7 @@ void PPPMDipole::set_grid_global(double dipole2)
     while (1) {
 
       // set grid dimension
+      
       nx_pppm = static_cast<int> (xprd/h_x);
       ny_pppm = static_cast<int> (yprd/h_y);
       nz_pppm = static_cast<int> (zprd_slab/h_z);
@@ -845,7 +847,8 @@ void PPPMDipole::set_grid_global(double dipole2)
       if (ny_pppm <= 1) ny_pppm = 2;
       if (nz_pppm <= 1) nz_pppm = 2;
 
-      //set local grid dimension
+      // set local grid dimension
+      
       int npey_fft,npez_fft;
       if (nz_pppm >= nprocs) {
         npey_fft = 1;
@@ -862,7 +865,7 @@ void PPPMDipole::set_grid_global(double dipole2)
       nzlo_fft = me_z*nz_pppm/npez_fft;
       nzhi_fft = (me_z+1)*nz_pppm/npez_fft - 1;
 
-      double df_kspace = compute_df_kspace_dipole(dipole2);
+      double df_kspace = compute_df_kspace_dipole();
 
       count++;
 
@@ -895,7 +898,7 @@ void PPPMDipole::set_grid_global(double dipole2)
    compute estimated kspace force error for dipoles
 ------------------------------------------------------------------------- */
 
-double PPPMDipole::compute_df_kspace_dipole(double dipole2)
+double PPPMDipole::compute_df_kspace_dipole()
 {
   double xprd = domain->xprd;
   double yprd = domain->yprd;
@@ -903,8 +906,7 @@ double PPPMDipole::compute_df_kspace_dipole(double dipole2)
   double zprd_slab = zprd*slab_volfactor;
   bigint natoms = atom->natoms;
   double qopt = compute_qopt_dipole();
-  //double df_kspace = sqrt(qopt/natoms)*mu2/(3.0*xprd*yprd*zprd_slab);
-  double df_kspace = sqrt(qopt/natoms)*dipole2/(3.0*xprd*yprd*zprd_slab);
+  double df_kspace = sqrt(qopt/natoms)*mu2/(3.0*xprd*yprd*zprd_slab);
   return df_kspace;
 }
 
@@ -1101,27 +1103,27 @@ void PPPMDipole::compute_gf_dipole()
    calculate f(x) for use in Newton-Raphson solver
 ------------------------------------------------------------------------- */
 
-//double PPPMDipole::newton_raphson_f()
-//{
-//  double xprd = domain->xprd;
-//  double yprd = domain->yprd;
-//  double zprd = domain->zprd;
-//  bigint natoms = atom->natoms;
-//
-//  double df_rspace,df_kspace;
-//  double vol = xprd*yprd*zprd;
-//  double a = cutoff*g_ewald;
-//  double rg2 = a*a;
-//  double rg4 = rg2*rg2;
-//  double rg6 = rg4*rg2;
-//  double Cc = 4.0*rg4 + 6.0*rg2 + 3.0;
-//  double Dc = 8.0*rg6 + 20.0*rg4 + 30.0*rg2 + 15.0;
-//  df_rspace = (mu2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) *
-//      sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) * exp(-rg2));
-//  df_kspace = compute_df_kspace_dipole();
-//
-//  return df_rspace - df_kspace;
-//}
+double PPPMDipole::newton_raphson_f()
+{
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+  bigint natoms = atom->natoms;
+
+  double df_rspace,df_kspace;
+  double vol = xprd*yprd*zprd;
+  double a = cutoff*g_ewald;
+  double rg2 = a*a;
+  double rg4 = rg2*rg2;
+  double rg6 = rg4*rg2;
+  double Cc = 4.0*rg4 + 6.0*rg2 + 3.0;
+  double Dc = 8.0*rg6 + 20.0*rg4 + 30.0*rg2 + 15.0;
+  df_rspace = (mu2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) *
+      sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) * exp(-rg2));
+  df_kspace = compute_df_kspace_dipole();
+
+  return df_rspace - df_kspace;
+}
 
 /* ----------------------------------------------------------------------
    find g_ewald parameter for dipoles based on desired accuracy
@@ -1184,7 +1186,7 @@ double PPPMDipole::derivf_dipole(double x, double Rc,
    calculate the final estimate of the accuracy
 ------------------------------------------------------------------------- */
 
-double PPPMDipole::final_accuracy_dipole(double dipole2)
+double PPPMDipole::final_accuracy_dipole()
 {
   double xprd = domain->xprd;
   double yprd = domain->yprd;
@@ -1193,7 +1195,7 @@ double PPPMDipole::final_accuracy_dipole(double dipole2)
   bigint natoms = atom->natoms;
   if (natoms == 0) natoms = 1; // avoid division by zero
 
-  double df_kspace = compute_df_kspace_dipole(mu2);
+  double df_kspace = compute_df_kspace_dipole();
 
   double a = cutoff*g_ewald;
   double rg2 = a*a;
@@ -1201,10 +1203,7 @@ double PPPMDipole::final_accuracy_dipole(double dipole2)
   double rg6 = rg4*rg2;
   double Cc = 4.0*rg4 + 6.0*rg2 + 3.0;
   double Dc = 8.0*rg6 + 20.0*rg4 + 30.0*rg2 + 15.0;
-  //double df_rspace = (mu2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) *
-  //  sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) *
-  //  exp(-rg2));
-  double df_rspace = (dipole2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) *
+  double df_rspace = (mu2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) *
     sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) *
     exp(-rg2));
 
@@ -2521,11 +2520,11 @@ double PPPMDipole::memory_usage()
   double bytes = nmax*3 * sizeof(double);
   int nbrick = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) *
     (nzhi_out-nzlo_out+1);
-  bytes += 6 * nfft_both * sizeof(double); // vg
-  bytes += nfft_both * sizeof(double); // greensfn
+  bytes += 6 * nfft_both * sizeof(double);   // vg
+  bytes += nfft_both * sizeof(double); 	     // greensfn
   bytes += nfft_both*5 * sizeof(FFT_SCALAR); // work*2*2
-  bytes += 9 * nbrick * sizeof(FFT_SCALAR); // ubrick*3 + vdbrick*6
-  bytes += nfft_both*7 * sizeof(FFT_SCALAR); //density_ffx*3 + work*2*2
+  bytes += 9 * nbrick * sizeof(FFT_SCALAR);  // ubrick*3 + vdbrick*6
+  bytes += nfft_both*7 * sizeof(FFT_SCALAR); // density_ffx*3 + work*2*2
 
   if (peratom_allocate_flag)
     bytes += 21 * nbrick * sizeof(FFT_SCALAR);
diff --git a/src/KSPACE/pppm_dipole.h b/src/KSPACE/pppm_dipole.h
index 8db28b540a..52bd2e5a9d 100644
--- a/src/KSPACE/pppm_dipole.h
+++ b/src/KSPACE/pppm_dipole.h
@@ -37,8 +37,8 @@ class PPPMDipole : public PPPM {
   double memory_usage();
 
  protected:
-  void set_grid_global(double);
-  //double newton_raphson_f();
+  void set_grid_global();
+  double newton_raphson_f(); 
 
   void allocate();
   void allocate_peratom();
@@ -76,7 +76,7 @@ class PPPMDipole : public PPPM {
   double find_gewald_dipole(double, double, bigint, double, double);
   double newton_raphson_f_dipole(double, double, bigint, double, double);
   double derivf_dipole(double, double, bigint, double, double);
-  double compute_df_kspace_dipole(double);
+  double compute_df_kspace_dipole();
   double compute_qopt_dipole();
   void compute_gf_dipole();
   void make_rho_dipole();
@@ -85,7 +85,7 @@ class PPPMDipole : public PPPM {
   void poisson_peratom_dipole();
   void fieldforce_ik_dipole();
   void fieldforce_peratom_dipole();
-  double final_accuracy_dipole(double dipole2);
+  double final_accuracy_dipole();
   void musum_musq();
 
 };
diff --git a/src/KSPACE/pppm_dipole_spin.cpp b/src/KSPACE/pppm_dipole_spin.cpp
index 4fde7ba101..a5aee7150c 100644
--- a/src/KSPACE/pppm_dipole_spin.cpp
+++ b/src/KSPACE/pppm_dipole_spin.cpp
@@ -50,8 +50,8 @@ using namespace MathSpecial;
 #define SMALL 0.00001
 #define EPS_HOC 1.0e-7
 
-enum{REVERSE_SP};
-enum{FORWARD_SP,FORWARD_SP_PERATOM};
+enum{REVERSE_MU};
+enum{FORWARD_MU,FORWARD_MU_PERATOM};
 
 #ifdef FFT_SINGLE
 #define ZEROF 0.0f
@@ -68,10 +68,12 @@ PPPMDipoleSpin::PPPMDipoleSpin(LAMMPS *lmp, int narg, char **arg) :
 {
   dipoleflag = 0;
   spinflag = 1;
-  group_group_enable = 0;
-
-  cg_dipole = NULL;
-  cg_peratom_dipole = NULL;
+  
+  hbar = force->hplanck/MY_2PI;         // eV/(rad.THz)
+  mub = 5.78901e-5;                     // in eV/T
+  mu_0 = 1.2566370614e-6;               // in T.m/A
+  mub2mu0 = mub * mub * mu_0;           // in eV
+  mub2mu0hbinv = mub2mu0 / hbar;        // in rad.THz
 }
 
 /* ----------------------------------------------------------------------
@@ -104,7 +106,10 @@ void PPPMDipoleSpin::init()
   // error check
 
   spinflag = atom->sp?1:0;
-
+  //qsum_qsq(0); // q[i] is probably not declared ?
+  //if (spinflag && q2)
+  //  error->all(FLERR,"Cannot use charges with Kspace style PPPMDipoleSpin");
+  
   triclinic_check();
 
   if (triclinic != domain->triclinic)
@@ -118,8 +123,8 @@ void PPPMDipoleSpin::init()
 
   if (!atom->sp) error->all(FLERR,"Kspace style requires atom attribute sp");
 
-  if (atom->sp && differentiation_flag == 1) error->all(FLERR,"Cannot (yet) use kspace_modify diff"
-       " ad with spins");
+  if (atom->sp && differentiation_flag == 1) error->all(FLERR,"Cannot (yet) use"
+     " kspace_modify diff ad with spins");
 
   if (spinflag && strcmp(update->unit_style,"metal") != 0)
     error->all(FLERR,"'metal' units have to be used with spins");
@@ -148,21 +153,19 @@ void PPPMDipoleSpin::init()
 
   int itmp = 0;
   double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
+  // probably not the correct extract here
   if (p_cutoff == NULL)
     error->all(FLERR,"KSpace style is incompatible with Pair style");
   cutoff = *p_cutoff;
 
   // kspace TIP4P not yet supported
-
+  // qdist = offset only for TIP4P fictitious charge
+  
+  qdist = 0.0; 
   if (tip4pflag)
     error->all(FLERR,"Cannot yet use TIP4P with PPPMDipoleSpin");
 
   scale = 1.0;
-  hbar = force->hplanck/MY_2PI;         // eV/(rad.THz)
-  mub = 5.78901e-5;                     // in eV/T
-  mu_0 = 1.2566370614e-6;               // in T.m/A
-  mub2mu0 = mub * mub * mu_0;           // in eV
-  mub2mu0hbinv = mub2mu0 / hbar;        // in rad.THz
   spsum_spsq();
   natoms_original = atom->natoms;
 
@@ -188,14 +191,14 @@ void PPPMDipoleSpin::init()
 
   GridComm *cgtmp = NULL;
   int iteration = 0;
-
+    
   while (order >= minorder) {
     if (iteration && me == 0)
       error->warning(FLERR,"Reducing PPPMDipoleSpin order b/c stencil extends "
                      "beyond nearest neighbor processor");
 
     compute_gf_denom();
-    set_grid_global(sp2);
+    set_grid_global();
     set_grid_local();
     if (overlap_allowed) break;
 
@@ -224,7 +227,7 @@ void PPPMDipoleSpin::init()
 
   // calculate the final accuracy
 
-  double estimated_accuracy = final_accuracy_dipole(sp2);
+  double estimated_accuracy = final_accuracy_dipole();
 
   // print stats
 
@@ -310,7 +313,7 @@ void PPPMDipoleSpin::compute(int eflag, int vflag)
 
   // return if there are no spins
 
-  if (spsqsum == 0.0) return;
+  if (musqsum == 0.0) return;
 
   // convert atoms from box to lamda coords
 
@@ -334,7 +337,7 @@ void PPPMDipoleSpin::compute(int eflag, int vflag)
   //   to fully sum contribution in their 3d bricks
   // remap from 3d decomposition to FFT decomposition
 
-  cg_dipole->reverse_comm(this,REVERSE_SP);
+  cg_dipole->reverse_comm(this,REVERSE_MU);
   brick2fft_dipole();
 
   // compute potential gradient on my FFT grid and
@@ -347,12 +350,12 @@ void PPPMDipoleSpin::compute(int eflag, int vflag)
   // all procs communicate E-field values
   // to fill ghost cells surrounding their 3d bricks
 
-  cg_dipole->forward_comm(this,FORWARD_SP);
+  cg_dipole->forward_comm(this,FORWARD_MU);
 
   // extra per-atom energy/virial communication
 
   if (evflag_atom) {
-    cg_peratom_dipole->forward_comm(this,FORWARD_SP_PERATOM);
+    cg_peratom_dipole->forward_comm(this,FORWARD_MU_PERATOM);
   }
 
   // calculate the force on my particles
@@ -374,7 +377,7 @@ void PPPMDipoleSpin::compute(int eflag, int vflag)
     energy = energy_all;
 
     energy *= 0.5*volume;
-    energy -= spsqsum*2.0*g3/3.0/MY_PIS;
+    energy -= musqsum*2.0*g3/3.0/MY_PIS;
     energy *= spscale;
   }
 
@@ -510,7 +513,7 @@ void PPPMDipoleSpin::fieldforce_ik_spin()
   double spx,spy,spz;
   double **x = atom->x;
   double **f = atom->f;
-  double **fm = atom->fm;
+  double **fm_long = atom->fm_long;
 
   int nlocal = atom->nlocal;
 
@@ -548,7 +551,7 @@ void PPPMDipoleSpin::fieldforce_ik_spin()
       }
     }
 
-    // convert M-field to mech. and mag. forces
+    // convert M-field and store mech. forces
 
     const double spfactor = mub2mu0 * scale;
     spx = sp[i][0]*sp[i][3];
@@ -558,13 +561,12 @@ void PPPMDipoleSpin::fieldforce_ik_spin()
     f[i][1] += spfactor*(vxy*spx + vyy*spy + vyz*spz);
     f[i][2] += spfactor*(vxz*spx + vyz*spy + vzz*spz);
 
+    // store long-range mag. precessions
+    
     const double spfactorh = mub2mu0hbinv * scale;
-    fm[i][0] += spfactorh*ex;
-    fm[i][1] += spfactorh*ey;
-    fm[i][2] += spfactorh*ez;
-
-    // create a new vector (in atom_spin style ?) to store long-range fm tables
-
+    fm_long[i][0] += spfactorh*ex;
+    fm_long[i][1] += spfactorh*ey;
+    fm_long[i][2] += spfactorh*ez;
   }
 }
 
@@ -708,9 +710,11 @@ void PPPMDipoleSpin::slabcorr()
   // add on mag. force corrections
 
   double ffact = spscale * (-4.0*MY_PI/volume);
-  double **fm = atom->fm;
+  //double **fm = atom->fm;
+  double **fm_long = atom->fm_long;
   for (int i = 0; i < nlocal; i++) {
-    fm[i][2] += ffact * spin_all;
+    //fm[i][2] += ffact * spin_all;
+    fm_long[i][2] += ffact * spin_all;
   }
 }
 
@@ -723,13 +727,13 @@ void PPPMDipoleSpin::spsum_spsq()
 {
   const int nlocal = atom->nlocal;
 
-  spsum = spsqsum = sp2 = 0.0;
+  musum = musqsum = mu2 = 0.0;
   if (atom->sp_flag) {
     double **sp = atom->sp;
     double spx, spy, spz;
     double spsum_local(0.0), spsqsum_local(0.0);
 
-    // not exactly the good loop: need to add norm of spins
+    // sum (direction x norm) of all spins
 
     for (int i = 0; i < nlocal; i++) {
       spx = sp[i][0]*sp[i][3];
@@ -739,12 +743,14 @@ void PPPMDipoleSpin::spsum_spsq()
       spsqsum_local += spx*spx + spy*spy + spz*spz;
     }
 
-    MPI_Allreduce(&spsum_local,&spsum,1,MPI_DOUBLE,MPI_SUM,world);
-    MPI_Allreduce(&spsqsum_local,&spsqsum,1,MPI_DOUBLE,MPI_SUM,world);
+    // store results into pppm_dipole quantities 
 
-    sp2 = spsqsum * mub2mu0;
+    MPI_Allreduce(&spsum_local,&musum,1,MPI_DOUBLE,MPI_SUM,world);
+    MPI_Allreduce(&spsqsum_local,&musqsum,1,MPI_DOUBLE,MPI_SUM,world);
+
+    mu2 = musqsum * mub2mu0;
   }
 
-  if (sp2 == 0 && comm->me == 0)
+  if (mu2 == 0 && comm->me == 0)
     error->all(FLERR,"Using kspace solver PPPMDipoleSpin on system with no spins");
 }
diff --git a/src/KSPACE/pppm_dipole_spin.h b/src/KSPACE/pppm_dipole_spin.h
index 347006e586..8d6c5d4eb2 100644
--- a/src/KSPACE/pppm_dipole_spin.h
+++ b/src/KSPACE/pppm_dipole_spin.h
@@ -42,7 +42,6 @@ class PPPMDipoleSpin : public PPPMDipole {
 
   // spin
 
-  double spsum,spsqsum,sp2;
   void make_rho_spin();
   void fieldforce_ik_spin();
   void fieldforce_peratom_spin();
diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp
index fb2b6dd797..477da613d2 100644
--- a/src/SPIN/atom_vec_spin.cpp
+++ b/src/SPIN/atom_vec_spin.cpp
@@ -106,7 +106,6 @@ void AtomVecSpin::grow_reset()
   sp = atom->sp; fm = atom->fm; fm_long = atom->fm_long;
 }
 
-
 /* ----------------------------------------------------------------------
    copy atom I info to atom J
 ------------------------------------------------------------------------- */
diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp
index b75f03212a..996bd3c2da 100644
--- a/src/SPIN/fix_nve_spin.cpp
+++ b/src/SPIN/fix_nve_spin.cpp
@@ -126,6 +126,7 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) :
   // initialize the magnetic interaction flags
 
   pair_spin_flag = 0;
+  long_spin_flag = 0;
   precession_spin_flag = 0;
   maglangevin_flag = 0;
   tdamp_flag = temp_flag = 0;
@@ -209,8 +210,16 @@ void FixNVESpin::init()
   if (count != npairspin)
     error->all(FLERR,"Incorrect number of spin pairs");
 
+  // set pair/spin and long/spin flags
+
   if (npairspin >= 1) pair_spin_flag = 1;
 
+  for (int i = 0; i<npairs; i++) {
+    if (force->pair_match("spin/long",0,i)) {
+      long_spin_flag = 1;
+    }
+  }
+
   // ptrs FixPrecessionSpin classes
 
   int iforce;
@@ -425,6 +434,7 @@ void FixNVESpin::ComputeInteractionsSpin(int i)
 
   double **sp = atom->sp;
   double **fm = atom->fm;
+  double **fm_long = atom->fm_long;
 
   // force computation for spin i
 
@@ -442,6 +452,14 @@ void FixNVESpin::ComputeInteractionsSpin(int i)
     }
   }
 
+  // update magnetic long-range components
+
+  if (long_spin_flag) {
+    fmi[0] += fm_long[i][0];
+    fmi[1] += fm_long[i][1];
+    fmi[2] += fm_long[i][2];
+  }
+
   // update magnetic precession interactions
 
   if (precession_spin_flag) {
diff --git a/src/SPIN/fix_nve_spin.h b/src/SPIN/fix_nve_spin.h
index afc1db14d6..565de13e92 100644
--- a/src/SPIN/fix_nve_spin.h
+++ b/src/SPIN/fix_nve_spin.h
@@ -58,7 +58,8 @@ friend class PairSpin;
   int nlocal_max;			// max value of nlocal (for lists size)
 
   int pair_spin_flag;			// magnetic pair flags
-  int precession_spin_flag;			// magnetic precession flags
+  int long_spin_flag;			// magnetic long-range flag
+  int precession_spin_flag;		// magnetic precession flags
   int maglangevin_flag;			// magnetic langevin flags
   int tdamp_flag, temp_flag;
 
diff --git a/src/SPIN/pair_spin_long.cpp b/src/SPIN/pair_spin_long.cpp
index ca4f3d5e24..efedea3247 100644
--- a/src/SPIN/pair_spin_long.cpp
+++ b/src/SPIN/pair_spin_long.cpp
@@ -81,6 +81,7 @@ PairSpinLong::~PairSpinLong()
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cut_spin_long);
+    memory->destroy(cutsq);
   }
 }
 
@@ -97,7 +98,6 @@ void PairSpinLong::settings(int narg, char **arg)
     error->all(FLERR,"Spin simulations require metal unit style");
 
   cut_spin_long_global = force->numeric(FLERR,arg[0]);
-  //cut_spin = force->numeric(FLERR,arg[0]);
   
   // reset cutoffs that have been explicitly set
   
@@ -126,7 +126,7 @@ void PairSpinLong::coeff(int narg, char **arg)
 
   if (strcmp(arg[2],"long") != 0)
     error->all(FLERR,"Incorrect args in pair_style command");
-  if (narg != 3) 
+  if (narg < 1 || narg > 4)
     error->all(FLERR,"Incorrect args in pair_style command");
   
   int ilo,ihi,jlo,jhi;
@@ -197,9 +197,9 @@ void PairSpinLong::init_style()
 double PairSpinLong::init_one(int i, int j)
 {
   if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
- 
+  
   cut_spin_long[j][i] = cut_spin_long[i][j];
-
+  
   return cut_spin_long_global;
 }
 
@@ -505,7 +505,8 @@ void PairSpinLong::allocate()
     for (int j = i; j <= n; j++)
       setflag[i][j] = 0;
 
-  memory->create(cut_spin_long,n+1,n+1,"pair:cut_spin_long");
+  memory->create(cut_spin_long,n+1,n+1,"pair/spin/long:cut_spin_long");
+  memory->create(cutsq,n+1,n+1,"pair/spin/long:cutsq");
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/atom.cpp b/src/atom.cpp
index cf4d20a71e..58b00712f7 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -98,7 +98,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
 
   // SPIN package
 
-  sp = fm = NULL;
+  sp = fm = fm_long = NULL;
   
   // USER-DPD
 
@@ -277,6 +277,7 @@ Atom::~Atom()
 
   memory->destroy(sp);
   memory->destroy(fm);
+  memory->destroy(fm_long);
 
   memory->destroy(vfrac);
   memory->destroy(s0);
diff --git a/src/atom.h b/src/atom.h
index 7e003dff5e..95c4a9c30f 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -65,6 +65,7 @@ class Atom : protected Pointers {
 
   double **sp;
   double **fm;
+  double **fm_long;
 
   // PERI package
 
diff --git a/src/pair.h b/src/pair.h
index f830b7c035..0b46d58782 100644
--- a/src/pair.h
+++ b/src/pair.h
@@ -61,7 +61,7 @@ class Pair : protected Pointers {
   int dispersionflag;            // 1 if compatible with LJ/dispersion solver
   int tip4pflag;                 // 1 if compatible with TIP4P solver
   int dipoleflag;                // 1 if compatible with dipole solver
-  int spinflag;                // 1 if compatible with spin long solver
+  int spinflag;			 // 1 if compatible with spin solver
   int reinitflag;                // 1 if compatible with fix adapt and alike
 
   int tail_flag;                 // pair_modify flag for LJ tail correction