diff --git a/doc/src/Eqs/pair_spin_dmi_interaction.jpg b/doc/src/Eqs/pair_spin_dmi_interaction.jpg
index 74807b105d711c2724629968b417c043a852f18b..1d15b2199add75eefecb424f75ea2dbd363a1afe 100644
Binary files a/doc/src/Eqs/pair_spin_dmi_interaction.jpg and b/doc/src/Eqs/pair_spin_dmi_interaction.jpg differ
diff --git a/doc/src/Eqs/pair_spin_dmi_interaction.tex b/doc/src/Eqs/pair_spin_dmi_interaction.tex
index 1f26c0829586ff09c126fc5aecdac0cd9b87f68e..5a5a776f030088efe17fcbe5e131e4686e8de2a3 100644
--- a/doc/src/Eqs/pair_spin_dmi_interaction.tex
+++ b/doc/src/Eqs/pair_spin_dmi_interaction.tex
@@ -5,7 +5,12 @@
 \begin{document}
 \begin{varwidth}{50in}
   \begin{equation}
-    \bm{H}_{dm} = \sum_{{ i,j}=1,i\neq j}^{N} \vec{D}\left(r_{ij}\right) \cdot\left(\vec{s}_{i}\times \vec{s}_{j}\right), \nonumber
+    \bm{H}_{dm} = -\sum_{{ i,j}=1,i\neq j}^{N} 
+    \left( \vec{e}_{ij} \times \vec{D} \right)
+    \cdot\left(\vec{s}_{i}\times \vec{s}_{j}\right), 
+    \nonumber
   \end{equation}
 \end{varwidth}
 \end{document}
+    \vec{D}\left(r_{ij}\right) 
+    {\rm ~and~}  \vec{D}\left(r_{ij}\right) = \vec{e}_{ij} \times \vec{D}  
diff --git a/doc/src/Eqs/pair_spin_me_interaction.jpg b/doc/src/Eqs/pair_spin_me_interaction.jpg
index 7b43a993af82e9d7e98783e96a2d5b16edbeb02e..8ecddc959226c0e94fa533f36ad84f3afced6ba6 100644
Binary files a/doc/src/Eqs/pair_spin_me_interaction.jpg and b/doc/src/Eqs/pair_spin_me_interaction.jpg differ
diff --git a/doc/src/Eqs/pair_spin_neel_functions.jpg b/doc/src/Eqs/pair_spin_neel_functions.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..57023634bf9de762613994e44174702f83c630d0
Binary files /dev/null and b/doc/src/Eqs/pair_spin_neel_functions.jpg differ
diff --git a/doc/src/Eqs/pair_spin_neel_functions.tex b/doc/src/Eqs/pair_spin_neel_functions.tex
new file mode 100644
index 0000000000000000000000000000000000000000..ead5a2b87a2fc879c1b3739b627f74b90f05dddf
--- /dev/null
+++ b/doc/src/Eqs/pair_spin_neel_functions.tex
@@ -0,0 +1,13 @@
+\documentclass[preview]{standalone}
+\usepackage{varwidth}
+\usepackage[utf8x]{inputenc}
+\usepackage{amsmath,amssymb,amsthm,bm}
+\begin{document}
+\begin{varwidth}{50in}
+  \begin{eqnarray}
+    g_1(r_{ij}) &=& g(r_{ij}) + \frac{12}{35} q(r_{ij}) \nonumber \\
+    q_1(r_{ij}) &=& \frac{9}{5} q(r_{ij}) \nonumber \\
+    q_2(r_{ij}) &=& - \frac{2}{5} q(r_{ij}) \nonumber
+  \end{eqnarray}
+\end{varwidth}
+\end{document}
diff --git a/doc/src/Eqs/pair_spin_neel_interaction.jpg b/doc/src/Eqs/pair_spin_neel_interaction.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..c30600d840e56bd04d994a0d16d73c80ddf760dd
Binary files /dev/null and b/doc/src/Eqs/pair_spin_neel_interaction.jpg differ
diff --git a/doc/src/Eqs/pair_spin_neel_interaction.tex b/doc/src/Eqs/pair_spin_neel_interaction.tex
new file mode 100644
index 0000000000000000000000000000000000000000..8b4ac33e732cccf30e83e56269617dac780e7cac
--- /dev/null
+++ b/doc/src/Eqs/pair_spin_neel_interaction.tex
@@ -0,0 +1,16 @@
+\documentclass[preview]{standalone}
+\usepackage{varwidth}
+\usepackage[utf8x]{inputenc}
+\usepackage{amsmath,amssymb,amsthm,bm}
+\begin{document}
+\begin{varwidth}{50in}
+  \begin{equation}
+    \mathcal{H}_{N\acute{e}el}=-\sum_{{ i,j=1,i\neq j}}^N g_1(r_{ij})\left(({\bm e}_{ij}\cdot {\bm s}_{i})({\bm e}_{ij}
+    \cdot {\bm s}_{j})-\frac{{\bm s}_{i}\cdot{\bm s}_{j}}{3} \right)
+    +q_1(r_{ij})\left( ({\bm e}_{ij}\cdot {\bm s}_{i})^2 -\frac{{\bm s}_{i}\cdot{\bm s}_{j}}{3}\right)
+    \left( ({\bm e}_{ij}\cdot {\bm s}_{i})^2 -\frac{{\bm s}_{i}\cdot{\bm s}_{j}}{3} \right)
+    + q_2(r_{ij}) \Big( ({\bm e}_{ij}\cdot {\bm s}_{i}) ({\bm e}_{ij}\cdot {\bm s}_{j})^3 + ({\bm e}_{ij}\cdot
+    {\bm s}_{j}) ({\bm e}_{ij}\cdot {\bm s}_{i})^3\Big) \nonumber
+  \end{equation}
+\end{varwidth}
+\end{document}
diff --git a/doc/src/compute_spin.txt b/doc/src/compute_spin.txt
index df75480430faf8ef2e35f04f83af40a6414e7c35..fcc764c14cb0f7cc88869655bf8d2df1f745c777 100644
--- a/doc/src/compute_spin.txt
+++ b/doc/src/compute_spin.txt
@@ -21,7 +21,7 @@ compute out_mag all compute/spin :pre
 
 [Description:]
 
-Define a computation that calculates the magnetic quantities for a system 
+Define a computation that calculates magnetic quantities for a system 
 of atoms having spins.
 
 This compute calculates 6 magnetic quantities.
@@ -30,11 +30,11 @@ The three first quantities are the x,y and z coordinates of the total magnetizat
 
 The fourth quantity is the norm of the total magnetization.
 
-The fifth one is referred to as the spin temperature, according
+The fifth quantity is the magnetic energy. 
+
+The sixth one is referred to as the spin temperature, according
 to the work of "(Nurdin)"_#Nurdin1. 
   
-The sixth quantity is the magnetic energy. 
-
 The simplest way to output the results of the compute spin calculation
 is to define some of the quantities as variables, and to use the thermo and
 thermo_style commands, for example:
diff --git a/doc/src/fix_langevin_spin.txt b/doc/src/fix_langevin_spin.txt
index 20e98eba872a1a76fa1b4ca00c531c693b6d0522..246ee023312932cc3325d5c87f84d5ebe389fe0c 100644
--- a/doc/src/fix_langevin_spin.txt
+++ b/doc/src/fix_langevin_spin.txt
@@ -44,15 +44,16 @@ in metal units).
 
 More details about this implementation are reported in "(Tranchida)"_#Tranchida2.
 
-Note: due to the form of the sLLG equation, this fix has to be the last defined 
-magnetic fix before the integration/spin fix. As an example:
+Note: due to the form of the sLLG equation, this fix has to be defined just
+before the nve/spin fix (and after all other magnetic fixes). 
+As an example:
 
 fix 1 all force/spin zeeman 0.01 0.0 0.0 1.0
 fix 2 all langevin/spin 300.0 0.01 21 
 fix 3 all integration/spin lattice yes :pre
 
 is correct, but defining a force/spin command after the langevin/spin command
-would send an error message. 
+would give an error message. 
 
 Note: The random # {seed} must be a positive integer.  A Marsaglia random
 number generator is used.  Each processor uses the input seed to
@@ -79,7 +80,7 @@ The {langevin/spin} fix is part of the SPIN package.
 This style is only enabled if LAMMPS was built with this package.
 See the "Making LAMMPS"_Section_start.html#start_3 section for more info.
 
-The numerical integration has to be performed with {fix/integration/spin} 
+The numerical integration has to be performed with {fix/nve/spin} 
 when {langevin/spin} is enabled. 
 
 This fix has to be the last defined magnetic fix before the integration fix 
diff --git a/doc/src/fix_nve_spin.txt b/doc/src/fix_nve_spin.txt
index 4f25003c5f88103a06573d29d690e8a7f42dbc72..6ccebcebf66292459283bf6b4270868e99c9d16d 100644
--- a/doc/src/fix_nve_spin.txt
+++ b/doc/src/fix_nve_spin.txt
@@ -39,7 +39,7 @@ the equations of motion of the spin lattice system, following the scheme:
 
 according to the implementation reported in "(Omelyan)"_#Omelyan1.
 
-A sectoring enables this scheme for parallel calculations. 
+A sectoring method enables this scheme for parallel calculations. 
 The implementation of this sectoring algorithm is reported 
 in "(Tranchida)"_#Tranchida1.
 
diff --git a/doc/src/pair_spin_dmi.txt b/doc/src/pair_spin_dmi.txt
index 9ad74d567e356ef45b71a800fccfad218432140a..9fe53df18a5cf81b6c1a838b968567d8220f1d36 100644
--- a/doc/src/pair_spin_dmi.txt
+++ b/doc/src/pair_spin_dmi.txt
@@ -30,31 +30,15 @@ between pairs of magnetic spins:
 :c,image(Eqs/pair_spin_dmi_interaction.jpg)
 
 where si and sj are two neighboring magnetic spins of two particles, 
-rij = ri - rj is the inter-atomic distance between the two particles,
-and D(rij) is the DM vector defining the intensity and the 
-sign of the exchange interaction.
+eij = (ri - rj)/|ri-rj| is the normalized separation vector between the 
+two particles, and D is the DM vector defining the intensity and the 
+sign of the interaction.
 
-This function is defined as:
-
-:c,image(Eqs/pair_spin_exchange_function.jpg)
-
-where a, b and d are the three constant coefficients defined in the associated 
-"pair_coeff" command.   
-
-The coefficients a, b, and d need to be fitted so that the function above matches with 
-the value of the DM interaction for the N neighbor shells taken 
-into account.
-
-Examples and more explanations about this function and its parametrization are reported 
-in "(Tranchida)"_#Tranchida5.
+Examples and more explanations about this interaction and its parametrization are 
+reported in "(Tranchida)"_#Tranchida5.
 
 From this DM interaction, each spin i will be submitted to a magnetic torque
-omega and its associated atom to a force F (for spin-lattice calculations only), 
-such as: 
-
-:c,image(Eqs/pair_spin_soc_dmi_forces.jpg)
-
-with h the Planck constant (in metal units).
+omega and its associated atom to a force F (for spin-lattice calculations only). 
 
 More details about the derivation of these torques/forces are reported in
 "(Tranchida)"_#Tranchida5.
diff --git a/doc/src/pair_spin_me.txt b/doc/src/pair_spin_me.txt
index 54c6d433003f72b14b9f1638889a5a886305f81c..8186775b7795f62110588c03d4ad9112ae609ef3 100644
--- a/doc/src/pair_spin_me.txt
+++ b/doc/src/pair_spin_me.txt
@@ -30,7 +30,7 @@ pairs of magnetic spins. According to the derivation reported in
 :c,image(Eqs/pair_spin_me_interaction.jpg)
 
 where si and sj are neighboring magnetic spins of two particles, 
-rij = ri - rj is the normalized separation vector between the 
+eij = (ri - rj)/|ri-rj| is the normalized separation vector between the 
 two particles, and E is an electric polarization vector.
 The norm and direction of E are giving the intensity and the 
 direction of a screened dielectric atomic polarization (in eV).
diff --git a/doc/src/pair_spin_neel.txt b/doc/src/pair_spin_neel.txt
index 5c82b6ef6f9c3481227aff56e6d21d1ea04c6b99..f7c9608a9398638378b0f9b1d039f56b07b199c8 100644
--- a/doc/src/pair_spin_neel.txt
+++ b/doc/src/pair_spin_neel.txt
@@ -27,14 +27,18 @@ pair_coeff 1 2 neel 4.0 0.0048 0.234 1.168 0.0 0.0 1.0 :pre
 Style {spin/neel} computes the Neel pair anisotropy model 
 between pairs of magnetic spins: 
 
-:c,image(Eqs/pair_spin_dmi_interaction.jpg)
+:c,image(Eqs/pair_spin_neel_interaction.jpg)
 
 where si and sj are two neighboring magnetic spins of two particles, 
 rij = ri - rj is the inter-atomic distance between the two particles,
-and D(rij) is the DM vector defining the intensity and the 
-sign of the exchange interaction.
+eij = (ri - rj)/|ri-rj| is their normalized separation vector 
+and g1, q1 and q2 are three functions defining the intensity of the
+dipolar and quadrupolar contributions, with:
 
-This function is defined as:
+:c,image(Eqs/pair_spin_neel_functions.jpg)
+
+With the functions g(rij) and q(rij) defined and fitted according to the same
+Bethe-Slater function used to fit the exchange interaction: 
 
 :c,image(Eqs/pair_spin_exchange_function.jpg)
 
@@ -42,19 +46,14 @@ where a, b and d are the three constant coefficients defined in the associated
 "pair_coeff" command.   
 
 The coefficients a, b, and d need to be fitted so that the function above matches with 
-the value of the DM interaction for the N neighbor shells taken 
-into account.
+the values of the magneto-elastic constant of the materials at stake. 
 
 Examples and more explanations about this function and its parametrization are reported 
-in "(Tranchida)"_#Tranchida6.
+in "(Tranchida)"_#Tranchida6. More examples of parametrization will be provided in 
+future work.
 
 From this DM interaction, each spin i will be submitted to a magnetic torque
-omega and its associated atom to a force F (for spin-lattice calculations only), 
-such as: 
-
-:c,image(Eqs/pair_spin_soc_dmi_forces.jpg)
-
-with h the Planck constant (in metal units).
+omega and its associated atom to a force F (for spin-lattice calculations only). 
 
 More details about the derivation of these torques/forces are reported in
 "(Tranchida)"_#Tranchida6.