Merge lp:~albertog/siesta/merge-OSSO into lp:siesta

[Alberto Garcia]

Merge the 'offsite spin-orbit' implementation by Ramon Cuadrado.

Reference:

R. Cuadrado and J. I. Cerda, J. Phys.: Condens. Matter 24, 086005,
(2012) (DOI:10.1088/0953-8984/24/8/086005)).

The new 'offsite' implementation is now the default when

   Spin { SO, SOC, S+O }
   SpinOrbit T (deprecated)

is specified in the input fdf file. To request the 'onsite' approximation,
use

   Spin { SO+onsite, SOC+onsite, S+O+onsite}

The only basic change from the existing 'onsite' implementation is in
the construction of the SO hamiltonian, which involves the full set of
lj KB projectors. These are constructed by new code in 'atom', and
processed in the new routine 'nlefsm_SO_off', which has roughly the
same structure as 'nlefsm', but constructs at the same time the 'ion'
(nl) and 'SO' pieces from the relativistic projectors.

This routine calls 'calc_Vj_offsiteSO', where VSO and Vion and the
corresponding forces are computed using the Clebsch–Gordan
coefficients needed to change from the basis |l,m,sigma> to |j,mj>.

The conventions for structure and signs in H and the DM are the same
as in the existing 'onsite' implementation, so there are no changes in
the diagonalization routines, or in the analysis routines and tools.

Eventually, the 'offsite' qualifier might be removed, as this is a
full spin-orbit implementation.

Note that this merge focuses on the core electronic-structure
functionality of 'full' spin-orbit coupling, and does not provide any
spin-orbit enhancements to the analysis tools.

Some polishing, shorter-running tests for SOC, etc.

This should be mature enough to go to the trunk, as an implementation of the 'offsite' SOC.
Further enhancements to analysis tools, etc, should be worked on in the other series,

[Nick Papior (nickpapior)]

There are changes in parts of the code which are not relevant.

Specifically the cold smearing method has been removed! Possibly due to a faulty merge.

Also there were changes to TD-DFT.

I would like a clean-up of comments which is scattered all-over.

Lastly, the use of the spin type should be used throughout (we might as well do it now).

My personal opinion would be to:
spin%SO_offsite change to spin%SO_off

spin%SO_offsite should *only* be true for SO calculations AND off site version. Otherwise false.
This allows easier if-statements throughout the code.

review: Needs Fixing

[Nick Papior (nickpapior)]

Great work Alberto!

I have amended comments again.

If you want, I can have a go on a PR for this branch?

review: Needs Fixing

[Alberto Garcia (albertog)]

Ok. Dfscf and rhoofd are going to be changed later.

On Fri, Apr 20, 2018, 11:59 Nick Papior <email address hidden> wrote:

> Review: Needs Fixing
>
> Great work Alberto!
>
> I have amended comments again.
>
> If you want, I can have a go on a PR for this branch?
>
> Diff comments:
>
> > === modified file 'Docs/siesta.tex'
> > --- Docs/siesta.tex 2018-04-17 13:07:32 +0000
> > +++ Docs/siesta.tex 2018-04-20 09:43:43 +0000
> > @@ -3721,10 +3728,12 @@
> >
> > \option[spin-orbit]%
> > \fdfindex*{Spin:spin-orbit}%
> > - Perform a calculation with spin-orbit coupling. This requires the
> > + Performs calculations including the spin--orbit coupling. By
> default the
> > + off--site SO option is set to \fdftrue. To perform an on--site SO
> calculations
> > + this option has to be {\bf spin-orbit+onsite}. This requires the
>
> No \bf, \textbf{...} ;)
>
> > pseudopotentials to be relativistic.
> >
> > - See Sect.~\ref{sec:spin-orbit}.
> > + See Sect.~\ref{sec:spin-orbit} for further specific spin--orbit
> options.
> >
> > \end{fdfoptions}
> >
> > @@ -3775,96 +3784,148 @@
> > \siesta\ includes the posibility to perform fully relativistic
> > calculations by means of the inclusion in the total Hamiltonian not
> > only the Darwin and velocity correction terms~(Scalar--Relativistic
> > -calculations), but also the spin-orbit~(SO) contribution. The
> > -implementation is based on the on-site SO approximation, where only
> > -the intra-SO contribution of each atom is taken into account. See
> > -\fdf{Spin} on how to turn on the spin-orbit coupling.
> > +calculations), but also the spin--orbit~(SO) contribution. There are
> > +two approaches regarding the SO formalism: on--site and off--site.
> > +Within the on--site approximation only the intra--atomic SO
> contribution is taken
> > +into account. In the off--site scheme additional neighboring
> > +interactions are also included in the SO term. By default, the
> off--site SO
> > +formalism is switched on, being necessary to change the \fdf{Spin} flag
> > +in the input file if the on--site approximation wants to be used. See
> > +\fdf{Spin} on how to handle the spin--orbit coupling.
> >
> > -The current implementation in \siesta\ has been implemented by
> > -Dr. Ram\'on Cuadrado based on the original on-site SO formalism and
> > -implementation developed by Prof. Jaime Ferrer, \textit{et al}~(L
> > +The on--site spin-orbit scheme in this version of \siesta\ has been
> implemented by
> > +Dr. Ram\'on Cuadrado based on the original on--site SO formalism and
> > +implementation developed by Prof. Jaime Ferrer and his collaborators
> \textit{et al}~(L
> > Fern\'andez--Seivane, M Oliveira, S Sanvito, and J Ferrer, Journal of
> > -Physics: Condensed Matter, 2006 vol. 18 pp. 7999; L
> > -Fern\'andez--Seivane and Jaime Ferrer, Phys. Rev. Lett. 99, 2007,
> > -183401).
> > -
> > -The inclusion of the SO term in the Hamiltonian~(and in the Density
> > -Matrix) will involve the increase of non-zero elements in their
> > -off-diagonal parts, i.e., for some $\mu\nu$ orbitals,
> > -H$^{\sigma\sigma'}_{\mu\nu}$(DM$^{\sigma\sigma'}_{\mu\nu}$)
> > +Physics: Condensed Matter, {\bf 18...

[Alberto Garcia (albertog)]

Please re-check physical results in view of the fix of the neglect of the highest-l SO component in the pseudopotential.

[Nick Papior (nickpapior)]

Very good Alberto.

I have re-assessed the comments and I only have two more points:

1) The WriteOrbMoms seems like a "bad" option name.

2) I can make a PR for clarifying H_so in sparse_matrices. I think having the classes and names of the variables to coincide with the type to be good (H_so_on_2D, H_so_off_2D) for instance?

[Ramon Cuadrado (ramon-cuadrado)]

Hi Alberto, all,

    Just finished some FePt and Pt2 tests (not Z-matrix) and it seems that for trunk-688--merge-OSSO-698 every test gives right results! I mean, anisotropy energies, MMs, forces and rotational tests are ok. Maybe Roberto wants to perform an additional test.

    The code does not finish properly. Just stops after write "Target enthalpy (eV/cell)" or Mulliken pop which I guess that will be fixed in the final merge...

    One more thing, I am not sure if the makefile within "offsite_SOC_" works fine. Just for info.

Best,

Ramón

review: Approve

[Alberto Garcia (albertog)]

> Very good Alberto.
>
> I have re-assessed the comments and I only have two more points:
>
> 1) The WriteOrbMoms seems like a "bad" option name.

I agree. We can change that in 4.1 and propagate.

>
> 2) I can make a PR for clarifying H_so in sparse_matrices. I think having the
> classes and names of the variables to coincide with the type to be good
> (H_so_on_2D, H_so_off_2D) for instance?

Please go ahead. At some point I would like the offSO H to be a real array with the same
conventions as the rest, to clarify the code throughout.

[Alberto Garcia (albertog)]

> Hi Alberto, all,
>
> Just finished some FePt and Pt2 tests (not Z-matrix) and it seems that for
> trunk-688--merge-OSSO-698 every test gives right results! I mean, anisotropy
> energies, MMs, forces and rotational tests are ok. Maybe Roberto wants to
> perform an additional test.

So you find no changes stemming from the fixing of the wrong-highest-l pseudo SO contribution?

>
> The code does not finish properly. Just stops after write "Target enthalpy
> (eV/cell)" or Mulliken pop which I guess that will be fixed in the final
> merge...

We need to look at this. Please provide fdf file and psf's.

>

[Ramon Cuadrado (ramon-cuadrado)]

> > Hi Alberto, all,
> >
> > Just finished some FePt and Pt2 tests (not Z-matrix) and it seems that
> for
> > trunk-688--merge-OSSO-698 every test gives right results! I mean, anisotropy
> > energies, MMs, forces and rotational tests are ok. Maybe Roberto wants to
> > perform an additional test.
>
> So you find no changes stemming from the fixing of the wrong-highest-l pseudo
> SO contribution?

  No changes.

>
> >
> > The code does not finish properly. Just stops after write "Target
> enthalpy
> > (eV/cell)" or Mulliken pop which I guess that will be fixed in the final
> > merge...
>
> We need to look at this. Please provide fdf file and psf's.

The PPs and the fdf’s are thosein the Test directory for FePt and Pt2...

>
> >

[Roberto Robles (roberto-robles)]

El 23/04/18 a las 11:35, Ramon Cuadrado escribió:
>>> Hi Alberto, all,
>>>
>>> Just finished some FePt and Pt2 tests (not Z-matrix) and it seems that
>> for
>>> trunk-688--merge-OSSO-698 every test gives right results! I mean, anisotropy
>>> energies, MMs, forces and rotational tests are ok. Maybe Roberto wants to
>>> perform an additional test.
>> So you find no changes stemming from the fixing of the wrong-highest-l pseudo
>> SO contribution?
> No changes.
>
>
>>> The code does not finish properly. Just stops after write "Target
>> enthalpy
>>> (eV/cell)" or Mulliken pop which I guess that will be fixed in the final
>>> merge...
>> We need to look at this. Please provide fdf file and psf's.
> The PPs and the fdf’s are thosein the Test directory for FePt and Pt2...
Is this the bug with state_analysis.F compiled with an intel compiler?
This is a very annoying problem, up to three people has contacted me
with that problem. If we can not fix it it should be stated in the
manual in the instructions for compiling, since usually people recycle
their arch.make instead of modifying the supplied intel.make

Roberto

>

[Ramon Cuadrado (ramon-cuadrado)]

Ok. Go ahead!

review: Approve

[Roberto Robles (roberto-robles)]

review: Approve

[Nick Papior (nickpapior)]

review: Approve

[lp:~albertog/siesta/merge-OSSO]

713. By Alberto Garcia on 2018-04-30

Update SOC_offsite notes in Docs

Also, remove comment in siesta.tex regarding structural optimization
with the on-site approximation.

(Thanks to Roberto Robles for input)

=== modified file 'Docs/Contributors.txt'
--- Docs/Contributors.txt	2016-09-28 07:49:18 +0000
+++ Docs/Contributors.txt	2018-04-30 20:53:35 +0000
@@ -30,7 +30,8 @@
 Thomas Archer
 Luis C. Balbas
 Xavier Blase
-Ramon Cuadrado
+Jorge I. Cerda,
+Ramon Cuadrado,
 Michele Ceriotti
 Raul de la Cruz
 Gabriel Fabricius

=== added file 'Docs/SOC_offsite_notes.txt'
--- Docs/SOC_offsite_notes.txt	1970-01-01 00:00:00 +0000
+++ Docs/SOC_offsite_notes.txt	2018-04-30 20:53:35 +0000
@@ -0,0 +1,33 @@
+Notes to the new 'offsite spin-orbit' implementation by Ramon Cuadrado.
+
+Reference:
+
+R. Cuadrado and J. I. Cerda,
+"Fully relativistic pseudopotential formalism under an atomic orbital
+basis: spin–orbit splittings and magnetic anisotropies",
+J. Phys.: Condens. Matter 24, 086005, (2012)
+(DOI:10.1088/0953-8984/24/8/086005)
+
+In this 'offsite' implementation the introduction of a fully
+relativistic Hamiltonian is done by using fully non-local operators
+for the pseudopotentials. In this way it is possible to avoid the
+'onsite' approximation needed to reduce the computational effort
+required when explicitly computing the LS term.
+
+The construction of a fully relativistic Hamiltonian involves the use
+of a full set of lj KB projectors. These are constructed by new code
+in 'atom', and processed in the new routine 'nlefsm_SO_off', which has
+roughly the same structure as 'nlefsm', but constructs at the same
+time the 'ion' and 'SO' pieces from the relativistic projectors.
+
+This routine calls 'calc_Vj_offsiteSO', where VSO and Vion and the
+corresponding forces are computed using the Clebsch–Gordan
+coefficients needed to change from the basis |l,m,sigma> to |j,mj>.
+
+The conventions for structure and signs in H and the DM are the same
+as in the existing 'onsite' implementation, so there are no changes in
+the diagonalization routines, or in the analysis routines and tools.
+
+Eventually, the 'offsite' qualifier might be removed, as this is a
+full spin-orbit implementation which involves a similar computational
+effort using less drastic approximations.

=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex	2018-04-24 10:24:02 +0000
+++ Docs/siesta.tex	2018-04-30 20:53:35 +0000
@@ -214,7 +214,8 @@
 Thomas Archer,
 Luis C. Balbas,
 Xavier Blase,
-Ramon Cuadrado,
+Jorge I. Cerd\'a,
+Ram\'on Cuadrado,
 Michele Ceriotti,
 Raul de la Cruz,
 Gabriel Fabricius,
@@ -2241,6 +2242,12 @@
 The default is \emph{one} KB projector from each angular momentum,
 constructed from the nodeless eigenfunction.
 
+For full spin-orbit calculations, the program generates $lj$
+projectors using the $l+1/2$ and $l-1/2$ components of the
+(relativistic) pseudopotentials. In this case the specification of the
+reference energies for projectors is not changed: only $l$ is
+relevant.
+
 \end{fdfentry}  
 
 
@@ -3251,7 +3258,7 @@
     coordinates rigidly to have them positive, by using
     \fdf{AtomicCoordinatesOrigin}.  See the
     \program{Sies2arc}\index{Sies2arc@\textsc{Sies2arc}} utility in the
-    \program{Util/} directory for generating \sysfile*{.arc} files for CERIUS animation.
+    \program{Util/} directory for generating \sysfile*{arc} files for CERIUS animation.
 
   \end{fdflogicalF}
 
@@ -3481,8 +3488,9 @@
 symmetry is valid in the absence of external magnetic fields or
 non-collinear/spin-orbit interaction.
 
-This flag should be used with care, as the code will produce wrong
-results if there is no support for the appropriate symmetrization.
+This flag is only honored for spinless or collinear-spin calculations,
+as the code will produce wrong results if there is no support for the
+appropriate symmetrization.
 
 The default value is \fdftrue unless: a) the option \fdf{Spin!Spiral} 
 is used. In this case time-reversal-symmetry is broken explicitly. b)
@@ -3721,10 +3729,12 @@
 
     \option[spin-orbit]%
     \fdfindex*{Spin:spin-orbit}%
-    Perform a calculation with spin-orbit coupling.  This requires the
+    Performs calculations including the spin-orbit coupling. By default the 
+    off-site SO option is set to \fdftrue. To perform an on-site SO calculations 
+    this option has to be \fdf*{spin-orbit+onsite}. This requires the
     pseudopotentials to be relativistic.
 
-    See Sect.~\ref{sec:spin-orbit}.
+    See Sect.~\ref{sec:spin-orbit} for further specific spin-orbit options.
 
   \end{fdfoptions}
 
@@ -3769,96 +3779,127 @@
 \end{fdflogicalF}
 
 
-\subsection{Spin--Orbit coupling}
+\subsection{Spin-Orbit coupling}
 \label{sec:spin-orbit}
 
 \siesta\ includes the posibility to perform fully relativistic
 calculations by means of the inclusion in the total Hamiltonian not
 only the Darwin and velocity correction terms~(Scalar--Relativistic
-calculations), but also the spin-orbit~(SO) contribution. The
-implementation is based on the on-site SO approximation, where only
-the intra-SO contribution of each atom is taken into account. See
-\fdf{Spin} on how to turn on the spin-orbit coupling.
+calculations), but also the spin-orbit~(SO) contribution. There are
+two approaches regarding the SO formalism: on-site and off-site.
+Within the on-site approximation only the intra-atomic SO
+contribution is taken into account. In the off-site scheme additional
+neighboring interactions are also included in the SO term. By default,
+the off-site SO formalism is switched on, being necessary to change
+the \fdf{Spin} flag in the input file if the on-site approximation
+wants to be used. See \fdf{Spin} on how to handle the spin-orbit
+coupling.
 
-The current implementation in \siesta\ has been implemented by
+The on-site spin-orbit scheme in this version of \siesta\ has been implemented by
 Dr. Ram\'on Cuadrado based on the original on-site SO formalism and
-implementation developed by Prof. Jaime Ferrer, \textit{et al}~(L
+implementation developed by Prof. Jaime Ferrer and his collaborators \textit{et al}~(L
 Fern\'andez--Seivane, M Oliveira, S Sanvito, and J Ferrer, Journal of
-Physics: Condensed Matter, 2006 vol. 18 pp. 7999; L
-Fern\'andez--Seivane and Jaime Ferrer, Phys. Rev. Lett. 99, 2007,
-183401).
-
-The inclusion of the SO term in the Hamiltonian~(and in the Density
-Matrix) will involve the increase of non-zero elements in their
-off-diagonal parts, i.e., for some $\mu\nu$ orbitals,
-H$^{\sigma\sigma'}_{\mu\nu}$(DM$^{\sigma\sigma'}_{\mu\nu}$)
-[$\sigma,\sigma'$=$\uparrow,\downarrow$] will be $\neq$0. This is
+Physics: Condensed Matter, \textbf{18}, 7999 (2006); L Fern\'andez--Seivane 
+and Jaime Ferrer, Phys. Rev. Lett. \textbf{99}, 183401 (2007)).
+
+The off-site scheme has been implemented by
+Dr. Ram\'on Cuadrado and Dr. Jorge I. Cerd\'a based on their initial 
+work~(R. Cuadrado and J. I. Cerd\'a ``Fully relativistic pseudopotential 
+formalism under an atomic orbital basis: spin-orbit splittings and 
+magnetic anisotropies'', J. Phys.: Condens. Matter \textbf{24}, 086005 (2012); 
+``In-plane/out-of-plane disorder influence on the magnetic anisotropy of 
+Fe$_{1-y}$Mn$_y$Pt-L1(0) bulk alloy'', R. Cuadrado, Kai Liu, Timothy 
+J. Klemmer and R. W. Chantrell, Applied Physics Letters, \textbf{108}, 
+123102 (2016)).  
+
+The inclusion of the SO term in the Hamiltonian (and in the Density
+Matrix) causes an increase in the number of non-zero elements in their
+off-diagonal parts, i.e., for some $(\mu,\nu)$ pair of basis
+orbitals, $\mathbf H^{\sigma\sigma'}_{\mu\nu}$ ($\mathbf{DM}^{\sigma\sigma'}_{\mu\nu}$)
+[$\sigma,\sigma'=\uparrow,\downarrow$] will be $\neq0$. This is
 mainly due to the fact that the $\mathbf L\cdot\mathbf S$ operator
-will promote the mixing between different spin-up/down components. The
-terms responsible of this matrices expansion are the
-exchange-correlation potential and the SO. The remaining terms such as
-the kinetic energy or Hartree contribution do not depend of the spin
-orientations and hence will be only added to the total
-Hamiltonian~(and DM) to their diagonal parts.
-
-The current SO formalism enables the possibility of several types of calculations:
+will promote the mixing between different spin-up/down components.
+In addition, these $\mathbf H^{\sigma\sigma'}_{\mu\nu}$ (and
+$\mathbf{DM}^{\sigma\sigma'}_{\mu\nu}$) elements will be complex, in contrast
+with typical polarized/non-polarized calculations where these
+matrices are purely real. Since the spin-up and spin-down manifolds
+are essentially mixed, the solver has to deal with matrices whose
+dimensions are twice as large as for the collinear (unmixed) spin
+problem. Due to this, we advise to take special
+attention to the memory needed to perform a spin-orbit calculation.
+
+
+Unless explicitly advised the following type of calculation can be carried out 
+regardless of whether on-site or off-site approximation is employed: 
 \begin{itemize}
   % 
   \item Selfconsistent calculations for gamma point as well as for
-  bulks~(Not yet implemented for optimizations).
-  % 
-  \item Magnetic Anisotropy Energy~(MAE) can be easily calculated. From
-  first principles calculations, MAE is obtained after subtract the
-  total selfconsistent energy in two different orientations, usually the
-  total energy associated with easy axis from the hard axis. In \siesta\
-  it is possible to perform several self-consistent calculations for
-  different magnetization orientations using the specific block
-  \fdf{DM.InitSpin} in the fdf file. In doing so one will be able to
-  include the initial orientation angles of the magnetization for each
-  atom, as well as an initial value of their net magnetic moments.
-  % 
-  \item By means of Mulliken analysis, after the self-consistent
+  bulks.
+  %
+  \item Structure optimizations 
+  %
+  %%% *** Incompatible... \item LDA+U calculations~(See Sect.\ref{sec:lda+u} for further info).
+  % 
+  \item Magnetic Anisotropy Energy~(MAE) can be easily
+    calculated. From first principles it is obtained after subtracting
+    the total selfconsistent energy calculated for two different
+    magnetic orientations. In \siesta\ it is possible to perform
+    calculations with different initial magnetic orderings
+    by means of the use of the block \fdf{DM.InitSpin} in the fdf
+    file. In doing so one will be able to include the initial
+    orientation angles of the magnetization for each atom, as well as
+    an initial value of its net magnetic moments.
+  % 
+  \item By means of Mulliken analysis, after the selfconsistent
   procedure, local spin and orbital moments can be calculated by means
   of the flags \fdf{WriteMullikenPop} and \fdf{WriteOrbMom}.
+  % 
 \end{itemize}
 
-Note: Due to the small SO energy value contribution to the total
-energy, the level of precision requiered to perform a proper fully
-relativistic calculation during the selfconsistent process is quite
-demanding. The following values must be carefully converged and
-checked for each specific system to assure that the results are
-accurate enough: \fdf{SCF.H!Tolerance} during the
-selfconsistency~(typically <10$^{-5}$eV), \fdf{ElectronicTemperature},
-\textbf{k}-point sampling and high values of
-\fdf{MeshCutoff}~(specifically for extended solids). In general, one
-can say that a good calculation will have high number of k--points,
-low \fdf{ElectronicTemperature}, extremely small \fdf{SCF.H!Tolerance}
-and high values of \fdf{MeshCutoff}.  We encourage the user to test
-carefully these options for each system. An additional point to take
-into account when the spin--orbit contribution is included is the
-mixing scheme employed. You are encouraged to use \fdf{SCF.Mix}
-\fdf*{hamiltonian} instead of the density matrix, due to the fact that
-the convergence speed increases considerably for the first case. In
-addition, the pseudopotentials have to be well generated and tested
-for each specific system and they have to be generated in their fully
-relativistic form and use the non-linear core corrections.
- 
+Note: Due to the small SO contribution to the total energy, the level
+of precision required to perform a proper fully relativistic
+calculation during the selfconsistent process is quite demanding. The
+following values must be carefully converged and checked for each
+specific system to assure that the results are accurate enough:
+\fdf{SCF.H!Tolerance} during the selfconsistency (typically between
+$10^{-3}\,\mathrm{eV}$ -- $10^{-4}\,\mathrm{eV}$),
+\fdf{ElectronicTemperature}, \textbf{k}--point sampling and high
+values of \fdf{MeshCutoff}~(specifically for extended solids). In
+general, one can say that a good calculation will have high number of
+\textbf{k}--points, low \fdf{ElectronicTemperature}, extremely small
+\fdf{SCF.H!Tolerance} and high values of \fdf{MeshCutoff}.  We
+encourage the user to test carefully these options for each system. An
+additional point to take into account is the mixing scheme
+employed. You are encouraged to use \fdf{SCF.Mix:hamiltonian}
+(currently is set up by default) instead of density matrix mixing,
+since it speeds up the convergence.  The pseudopotentials have to be
+properly generated and tested for each specific system and they have
+to be in their fully relativistic form, together with the non-linear
+core corrections. Finally it is worth to mention that the
+selfconsistent convergence for some non-highly symmetric
+magnetizations directions with respect to the physical symmetry axis
+could still be difficult.
+
 \begin{fdfentry}{Spin!OrbitStrength}[real]<1.0>
 
-  It allows to vary the strength of the spin-orbit interaction from
-  zero to any positive value, including the physical value. This flag
-  is only active when \fdf{Spin} is set to \fdf*{spin-orbit}.
-
+  It allows to vary the strength of the 
+  spin-orbit interaction from zero to any positive value. It can be
+  used for both the on-site and off-site SOC flavors, but only for
+  debugging and testing purposes, as the only physical value is 1.0.
+  Note that this feature is implemented by modifying the SO parts of the
+  semilocal potentials read from a \code{.psf} file. Care must be
+  taken when re-using any \code{.ion} files produced.
+  
 \end{fdfentry}
 
 \begin{fdflogicalF}{WriteOrbMom}
 
-  If \fdftrue, a table is provided in the main output file, which
+  If \fdftrue, a table is provided in the output file that
   includes an estimation of the vector orbital magnetic
-  moments, in units of the Bohr magneton, projected onto each orbital
-  and also onto each atom. The estimation for the orbital moments is
-  based on a two-center approximation, and makes use of the Mulliken
-  population analysis.
+  moments, in units of the Bohr magneton, projected 
+  onto each orbital and also onto each atom. The estimation for the 
+  orbital moments is based on a two-center approximation, and makes use 
+  of the Mulliken population analysis.
 
   If \fdf{MullikenInScf} is \fdftrue, this information is printed at
   every scf step.
@@ -3866,6 +3907,7 @@
 \end{fdflogicalF}
 
 
+
 \subsection{The self-consistent-field loop}
 
 \textbf{IMPORTANT NOTE: Convergence of the Kohn-Sham energy and forces}
@@ -4988,7 +5030,7 @@
   \index{reading saved data!density matrix}
   
   Instructs to read the density matrix stored in file
-  \sysfile{.DM} by a previous run.
+  \sysfile{DM} by a previous run.
   
   \siesta\ will continue even if \sysfile*{DM} is not found.
 
@@ -9014,8 +9056,8 @@
 Enable an experimental timer which is based on wall time on the master
 node and is aware of the tree-structure of the timed sections. At the
 end of the program, a report is generated in the output file, and a
-{\tt time.json} file in JSON format is also written. \index{JSON
-  timing report@{\bf JSON timing report}} This file can be used by
+\file{time.json} file in JSON format is also written. \index{JSON
+  timing report@\textbf{JSON timing report}} This file can be used by
 third-party scripts to process timing data.
 
   \note, if used with the PEXSI solver (see Sec.~\ref{SolverPEXSI})
@@ -10567,9 +10609,9 @@
 
 \begin{fdflogicalF}{TDED.Saverho}
 
-If \fdftrue\ the instantaneous time-dependent density is saved to
-\texttt{ <istep>}.\texttt{TDRho}\index{TDRho@{\bf TDRho}} after every
-\fdf{TDED.Nsaverho} number of steps.
+  If \fdftrue\ the instantaneous time-dependent density is saved to
+  \file{<istep>.TDRho} after every \fdf{TDED.Nsaverho} number of
+  steps.
 
 \end{fdflogicalF}
 

=== modified file 'Src/Makefile'
--- Src/Makefile	2018-04-24 10:24:02 +0000
+++ Src/Makefile	2018-04-30 20:53:35 +0000
@@ -635,10 +635,11 @@
 compute_dm.o: normalize_dm.o ordern.o parallel.o precision.o siesta_geom.o
 compute_dm.o: siesta_options.o sparse_matrices.o sys.o units.o
 compute_ebs_shift.o: m_mpi_utils.o parallel.o precision.o
-compute_energies.o: atomlist.o class_SpData1D.o class_SpData2D.o dhscf.o
-compute_energies.o: files.o m_dipol.o m_energies.o m_mpi_utils.o m_ntm.o
-compute_energies.o: m_rhog.o m_spin.o precision.o siesta_geom.o
-compute_energies.o: siesta_options.o sparse_matrices.o
+compute_energies.o: atomlist.o class_SpData1D.o class_SpData2D.o
+compute_energies.o: class_SpData2D.o dhscf.o files.o m_dipol.o m_energies.o
+compute_energies.o: m_mpi_utils.o m_ntm.o m_rhog.o m_spin.o parallel.o
+compute_energies.o: precision.o siesta_geom.o siesta_options.o
+compute_energies.o: sparse_matrices.o
 compute_max_diff.o: m_mpi_utils.o precision.o
 compute_norm.o: m_mpi_utils.o m_spin.o precision.o sparse_matrices.o
 compute_pw_matrix.o: alloc.o m_planewavematrix.o m_planewavematrixvar.o
@@ -725,15 +726,15 @@
 fft.o: alloc.o fft1d.o m_timer.o mesh.o parallel.o parallelsubs.o precision.o
 fft.o: sys.o
 fft1d.o: parallel.o precision.o sys.o
-final_H_f_stress.o: alloc.o atomlist.o class_SpData2D.o compute_max_diff.o
-final_H_f_stress.o: dnaefs.o files.o grdsam.o kinefsm.o ldau.o ldau_specs.o
-final_H_f_stress.o: m_dipol.o m_energies.o m_forces.o m_gamma.o m_hsx.o
-final_H_f_stress.o: m_mpi_utils.o m_ncdf_siesta.o m_ntm.o m_spin.o m_steps.o
-final_H_f_stress.o: m_stress.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-final_H_f_stress.o: m_ts_options.o metaforce.o molecularmechanics.o naefs.o
-final_H_f_stress.o: nlefsm.o overfsm.o parallel.o siesta_geom.o
-final_H_f_stress.o: siesta_options.o sparse_matrices.o spinorbit.o sys.o
-final_H_f_stress.o: units.o
+final_H_f_stress.o: alloc.o atomlist.o class_SpData2D.o class_SpData2D.o
+final_H_f_stress.o: compute_max_diff.o dnaefs.o files.o grdsam.o kinefsm.o
+final_H_f_stress.o: ldau.o ldau_specs.o m_dipol.o m_energies.o m_forces.o
+final_H_f_stress.o: m_gamma.o m_hsx.o m_mpi_utils.o m_ncdf_siesta.o m_ntm.o
+final_H_f_stress.o: m_spin.o m_steps.o m_stress.o m_ts_global_vars.o m_ts_io.o
+final_H_f_stress.o: m_ts_kpoints.o m_ts_options.o metaforce.o
+final_H_f_stress.o: molecularmechanics.o naefs.o nlefsm.o overfsm.o parallel.o
+final_H_f_stress.o: siesta_geom.o siesta_options.o sparse_matrices.o
+final_H_f_stress.o: spinorbit.o sys.o units.o
 find_kgrid.o: alloc.o minvec.o parallel.o precision.o units.o
 fire_optim.o: alloc.o m_fire.o m_mpi_utils.o parallel.o precision.o
 fire_optim.o: siesta_options.o units.o
@@ -932,7 +933,7 @@
 m_sparsity_handling.o: class_SpData2D.o class_Sparsity.o geom_helper.o
 m_sparsity_handling.o: intrinsic_missing.o m_interpolate.o m_region.o
 m_sparsity_handling.o: precision.o
-m_spin.o: alloc.o parallel.o precision.o sys.o units.o
+m_spin.o: alloc.o files.o m_cite.o parallel.o precision.o sys.o units.o
 m_stress.o: precision.o
 m_supercell.o: atom_graph.o class_OrbitalDistribution.o class_SpData2D.o
 m_supercell.o: intrinsic_missing.o parallel.o parallelsubs.o precision.o
@@ -1135,7 +1136,7 @@
 new_matel.o: alloc.o errorf.o interpolation.o matel_registry.o parallel.o
 new_matel.o: precision.o radfft.o spher_harm.o sys.o
 nlefsm.o: alloc.o atm_types.o atmfuncs.o atomlist.o chemical.o mneighb.o
-nlefsm.o: new_matel.o parallel.o parallelsubs.o precision.o
+nlefsm.o: new_matel.o parallel.o parallelsubs.o precision.o sparse_matrices.o
 normalize_dm.o: atomlist.o m_mpi_utils.o m_spin.o parallel.o precision.o
 normalize_dm.o: siesta_options.o sparse_matrices.o sys.o
 obc.o: alloc.o precision.o
@@ -1227,15 +1228,17 @@
 setatomnodes.o: alloc.o parallel.o precision.o spatial.o sys.o
 setspatial.o: alloc.o parallel.o precision.o spatial.o
 setup_H0.o: alloc.o atmfuncs.o atomlist.o class_SpData1D.o class_SpData2D.o
-setup_H0.o: dhscf.o dnaefs.o kinefsm.o m_energies.o m_mpi_utils.o m_ntm.o
-setup_H0.o: m_spin.o metaforce.o molecularmechanics.o naefs.o nlefsm.o
-setup_H0.o: siesta_geom.o siesta_options.o sparse_matrices.o spinorbit.o
+setup_H0.o: class_SpData2D.o dhscf.o dnaefs.o kinefsm.o m_energies.o
+setup_H0.o: m_mpi_utils.o m_ntm.o m_spin.o metaforce.o molecularmechanics.o
+setup_H0.o: naefs.o nlefsm.o siesta_geom.o siesta_options.o sparse_matrices.o
+setup_H0.o: spinorbit.o
 setup_hamiltonian.o: alloc.o atmfuncs.o atomlist.o class_SpData1D.o
-setup_hamiltonian.o: class_SpData2D.o dhscf.o files.o ldau.o ldau_specs.o
-setup_hamiltonian.o: m_dipol.o m_energies.o m_gamma.o m_hsx.o m_mpi_utils.o
-setup_hamiltonian.o: m_ntm.o m_partial_charges.o m_rhog.o m_spin.o m_steps.o
-setup_hamiltonian.o: m_stress.o metaforce.o molecularmechanics.o parallel.o
-setup_hamiltonian.o: siesta_geom.o siesta_options.o sparse_matrices.o sys.o
+setup_hamiltonian.o: class_SpData2D.o class_SpData2D.o dhscf.o files.o ldau.o
+setup_hamiltonian.o: ldau_specs.o m_dipol.o m_energies.o m_gamma.o m_hsx.o
+setup_hamiltonian.o: m_mpi_utils.o m_ntm.o m_partial_charges.o m_rhog.o
+setup_hamiltonian.o: m_spin.o m_steps.o m_stress.o metaforce.o
+setup_hamiltonian.o: molecularmechanics.o parallel.o siesta_geom.o
+setup_hamiltonian.o: siesta_options.o sparse_matrices.o sys.o
 setup_ordern_indexes.o: alloc.o domain_decom.o parallel.o spatial.o
 shaper.o: atmfuncs.o mneighb.o precision.o
 show_distribution.o: atomlist.o parallel.o parallelsubs.o siesta_geom.o sys.o
@@ -1305,26 +1308,27 @@
 siesta_tddft.o: wavefunctions.o
 sparse_matrices.o: alloc.o class_Fstack_Pair_Geometry_SpData2D.o
 sparse_matrices.o: class_OrbitalDistribution.o class_SpData1D.o
-sparse_matrices.o: class_SpData2D.o class_Sparsity.o precision.o
+sparse_matrices.o: class_SpData2D.o class_SpData2D.o class_Sparsity.o
+sparse_matrices.o: precision.o
 spatial.o: precision.o
 spher_harm.o: alloc.o precision.o sys.o
 spinorbit.o: atm_types.o atmfuncs.o atmparams.o basis_types.o m_mpi_utils.o
-spinorbit.o: parallel.o parallelsubs.o precision.o pseudopotential.o sys.o
+spinorbit.o: parallel.o parallelsubs.o precision.o pseudopotential.o
 state_analysis.o: atomlist.o born_charge.o flook_siesta.o m_energies.o
 state_analysis.o: m_fixed.o m_forces.o m_ntm.o m_spin.o m_stress.o
 state_analysis.o: m_wallclock.o parallel.o remove_intramol_pressure.o
 state_analysis.o: siesta_cml.o siesta_geom.o siesta_options.o sparse_matrices.o
 state_analysis.o: units.o write_subs.o zmatrix.o
 state_init.o: alloc.o atomlist.o class_Data2D.o class_SpData1D.o
-state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
-state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
-state_init.o: iotdxv.o ioxv.o kpoint_grid.o ldau_specs.o m_chess.o m_energies.o
-state_init.o: m_eo.o m_gamma.o m_mixing.o m_mixing_scf.o m_mpi_utils.o
-state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
-state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
-state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
-state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
-state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: class_SpData2D.o class_SpData2D.o class_Sparsity.o
+state_init.o: create_Sparsity_SC.o domain_decom.o files.o hsparse.o
+state_init.o: iodm_netcdf.o iodmhs_netcdf.o iotdxv.o ioxv.o kpoint_grid.o
+state_init.o: ldau_specs.o m_chess.o m_energies.o m_eo.o m_gamma.o m_mixing.o
+state_init.o: m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o m_pivot_methods.o
+state_init.o: m_rmaxh.o m_sparse.o m_sparsity_handling.o m_spin.o m_steps.o
+state_init.o: m_supercell.o m_test_io.o m_ts_charge.o m_ts_electype.o
+state_init.o: m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o m_ts_options.o
+state_init.o: m_ts_sparse.o m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
 state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
 state_init.o: siesta_options.o sparse_matrices.o sys.o units.o write_subs.o
 state_init.o: zmatrix.o

=== modified file 'Src/atm_types.f'
--- Src/atm_types.f	2018-04-07 19:24:04 +0000
+++ Src/atm_types.f	2018-04-30 20:53:35 +0000
@@ -23,12 +23,12 @@
       integer, parameter, public  :: maxnorbs = 100
 !       Maximum number of nlm orbitals
 !
-      integer, parameter, public  :: maxn_pjnl = 10
+      integer, parameter, public  :: maxn_pjnl = 20
 !       Maximum number of projectors (not counting different "m" copies)
       integer, parameter, public  :: maxn_orbnl = 200
 !       Maximum number of nl orbitals (not counting different "m" copies)
 !       Now very large to accommodate filteret basis sets
-      integer, parameter, public  :: maxnprojs = 50
+      integer, parameter, public  :: maxnprojs = 100
 !       Maximum number of nlm projectors
 !
 
@@ -66,9 +66,11 @@
 !             1 to the total number of projectors at that l.
 !             
 !
+         logical                         ::  lj_projs = .false.
          integer                         ::  n_pjnl=0   ! num of "nl" projs
          integer                         ::  lmax_projs=0 ! l cutoff for projs
          integer, dimension(maxn_pjnl)   ::  pjnl_l     ! l of each nl proj
+         real(dp), dimension(maxn_pjnl)  ::  pjnl_j     ! j of each nl proj
          integer, dimension(maxn_pjnl)   ::  pjnl_n     ! n of each nl proj
          real(dp), dimension(maxn_pjnl)
      $                                   ::  pjnl_ekb   ! energy of
@@ -92,9 +94,10 @@
          integer, dimension(maxnprojs)   ::  pj_index
          integer, dimension(maxnprojs)   ::  pj_n
          integer, dimension(maxnprojs)   ::  pj_l
+         real(dp), dimension(maxnprojs)  ::  pj_j
          integer, dimension(maxnprojs)   ::  pj_m
          integer, dimension(maxnprojs)   ::  pj_gindex
-!
+!----------------------------
 !        LDA+U Projectors
 !        Here we follow the scheme used for the KB projectors
 !        

=== modified file 'Src/atmfuncs.f'
--- Src/atmfuncs.f	2016-04-29 07:31:13 +0000
+++ Src/atmfuncs.f	2018-04-30 20:53:35 +0000
@@ -20,6 +20,8 @@
       use atm_types
       use radial, only: rad_get, rad_func
       use spher_harm, only: rlylm
+ 
+
 
       implicit none 
 !

=== modified file 'Src/atmparams.f'
--- Src/atmparams.f	2016-01-25 16:00:16 +0000
+++ Src/atmparams.f	2018-04-30 20:53:35 +0000
@@ -21,7 +21,11 @@
 C INTEGER  NKBMX    : Maximum number of Kleinman-Bylander projectors
 C                     for each angular momentum
 
-         integer, parameter, public  :: nkbmx  =    2  
+C         For the off-site SO calculation plus semicore states
+C         there will be at least 4 KBs for each l angular momentum
+C         (for each l shell we have J = l +/- 1/2 )
+         integer, parameter, public  :: nkbmx  =    4
+
 
 C INTEGER  NSMX    : Maximum number of semicore shells for each angular
 C                    momentum present in the atom ( for normal atom nsmx=0)

=== modified file 'Src/atom.F'
--- Src/atom.F	2018-04-26 13:16:13 +0000
+++ Src/atom.F	2018-04-30 20:53:35 +0000
@@ -97,7 +97,7 @@
      .                       rco,lambda_in,atm_label,npolorb,semic,
      .                       nsemic,cnfigmx,charge_in,smass,basistype, 
      .                       ISIN,RINN,VCTE,qcoe,qyuk,qwid,
-     .                       split_norm,filtercut_in,basp, spp)
+     .                       split_norm,filtercut_in,basp,spp,lj_projs)
 
         use atm_types, only: species_info
         
@@ -105,6 +105,7 @@
 
         type(basis_def_t), pointer   :: basp
         type(species_info), intent(inout)  :: spp
+        logical, intent(in)                :: lj_projs
 
         integer, intent(in) :: izin  
 !       Atomic number (if IZ is negative it will be regarded as a set of 
@@ -242,6 +243,7 @@
         real(dp)
      .    rofi(nrmax), drdi(nrmax), s(nrmax),
      .    vps(nrmax,0:lmaxd), rphi(nrmax,0:lmaxd,nsemx),
+     .    vps_u(nrmax,0:lmaxd),
      .    vlocal(nrmax), vxc(nrmax), ve(nrmax),
      .    rho(nrmax), chcore(nrmax), rho_PAO(nrmax), auxrho(nrmax),
      .    vePAO(nrmax), qPAO(0:lmaxd,nsemx), chlocal(nrmax),
@@ -329,7 +331,7 @@
 !
 ! Reading pseudopotentials
 ! 
-          call read_vps(lmxo, lmxkb, nrval,a,b,rofi,drdi,s,vps,
+          call read_vps(lmxo, lmxkb, nrval,a,b,rofi,drdi,s,vps,vps_u,
      .        rho, chcore, zval, chgvps, nicore, irel, icorr,basp) 
 
           do ir=1,nrval
@@ -569,8 +571,8 @@
 ! Calculation of the Kleinman-Bylander projector functions
 ! 
             call KBgen(is, a,b,rofi,drdi,s,
-     .                 vps, vlocal, ve, nrval, Zval, lmxkb, 
-     .                 nkbl, erefkb, nkb, spp) 
+     .                 vps, vps_u, vlocal, ve, nrval, Zval, lmxkb, 
+     .                 nkbl, erefkb, nkb, spp, lj_projs) 
 
           elseif(flting.lt.0.0d0) then 
 !
@@ -1856,7 +1858,7 @@
 
 
       subroutine KBproj(ikb,rofi,drdi,vps,vlocal,nrwf,l,rphi,
-     .    dkbcos,ekb,proj,nrc)  
+     .    dkbcos,ekb,proj,nrc,rphi2,vii)  
 C
 C  This routine calculates the Kleinman-Bylander projector
 C  with angular momentum l.
@@ -1872,6 +1874,7 @@
      .                           rphi(*), vlocal(*)
         real(dp), intent(out) :: proj(*), dkbcos, ekb
         integer, intent(out)  :: nrc
+        real(dp), intent(inout) :: rphi2(:,:), vii(:)
 C
 C Internal variables
 C
@@ -1880,8 +1883,6 @@
      .    dnrm, vl, vphi, avgv, r, phi, dknrm,
      .    dincv, rc, sum, vij(nkbmx)
 
-        real(dp), save :: rphi2(nrmax,nkbmx), vii(nkbmx)
-
         integer   ir, jkb
 
         real(dp), parameter  :: eps=1.0d-6
@@ -2702,7 +2703,7 @@
 
 !
         subroutine read_vps(lmxo, lmxkb,
-     .             nrval,a,b,rofi,drdi,s,vps,
+     .             nrval,a,b,rofi,drdi,s,vps, vps_u,
      .             rho, chcore, zval, chgvps,
      .             nicore, irel, icorr,basp)
 
@@ -2718,14 +2719,14 @@
       
         real(dp)
      .    rofi(nrmax), drdi(nrmax), s(nrmax), vps(nrmax,0:lmaxd),
+     .    vps_u(nrmax,0:lmaxd),
      .    rho(nrmax), chcore(nrmax) 
 
         real(dp)
      .    a, b, zval
-           
         integer  nrval, lmxo, lmxkb 
-
         character nicore*4, irel*3, icorr*2
+        integer :: nup_end
 C
 C Internal variables
 C
@@ -2736,12 +2737,14 @@
      .    ea, rpb, chgvps
 
         integer  
-     .    nr, nodd, lmax, linput, npotd, npotu, ndown, l, ir, i
+     .    nr, nodd, lmax, linput, npotd, npotu, ndown, l, ir, i, nup
 
         character method(6)*10,text*70
 
         type(pseudopotential_t), pointer :: vp
 
+        real(dp) :: so_strength
+        
         vp => basp%pseudopotential
 
         icorr = vp%icorr
@@ -2792,12 +2795,8 @@
         write(6,'(7a)')   
      .    'read_vps: ',method(1),(method(i),i=3,6)
           
-C We are going to find the charge configuration
-C used for the pseudopotential generation using the information given in
-C the 'text' variable.
-
+        !Total charge density used for the pseudopotential generation
         chgvps = vp%gen_zval
-
         write(6,'(a,f10.5)') 'Total valence charge: ', chgvps
 
         if (nicore.ne.'nc  ') then
@@ -2842,9 +2841,9 @@
         enddo 
 
 
-!       Ionic pseudopotentials   (Only 'down' used)
+!       Ionic pseudopotentials (down and up)
 
-        do 20 ndown=1,lmax+1
+        do ndown=1,lmax+1
           l = vp%ldown(ndown) 
           if (l.ne.ndown-1) then
             write(6,'(a)')
@@ -2858,40 +2857,51 @@
           enddo
           vps(1,l) = vps(2,l)     ! AG
 
-  20    continue
-
+       enddo
+       
+       vps_u = 0.0_dp
+       ! For backward compatibility with original OSSO branch,
+       ! put 'keep-npotu-bug T' in fdf file
+       ! To be removed altogether for releases
+       if (fdf_get('keep-npotu-bug',.false.))  then
+          nup_end =  min(lmax,npotu-1)
+          write(6,'(a)') 'atom: *** Buggy code in highest-l ' //
+     $                   'SO component setup enabled ***'
+       else
+          nup_end = npotu
+       endif
+
+       ! Allow an overall factor for the SO part (for debugging only!)
+       ! Note that (notably for the lj-based SO implementation) this
+       ! trick is based on the availability of a semilocal pseudo.  Care
+       ! should also be taken with the onward use of any .ion files
+       ! produced.
+       
+       so_strength = fdf_get('spin-orbit-strength', 1.0_dp)
+       if (so_strength /= 1.0_dp)  then
+          write(6,'(a,f10.6)') 'atom: *** SO enhancement factor: ',
+     $         so_strength
+       endif
+
+       do nup=1,nup_end
+          l = vp%lup(nup)
+          vp%vup(nup,1:nrval) = so_strength * vp%vup(nup,1:nrval)
+          vps_u(1:nrval,l) = vp%vup(nup,1:nrval)
+          do ir=2,nrval
+            vps_u(ir,l)=vps_u(ir,l)/rofi(ir)
+          enddo
+          vps_u(1,l) = vps_u(2,l)    
+       enddo
 
 !       Core and valence charge density
 
         chcore(1:nrval) = vp%chcore(1:nrval)
         rho(1:nrval) = vp%chval(1:nrval)
-C
-C Obtain an ionic-pseudopotential if core correction for Hartree
-C potential
-!
-! AG: OBSOLETE, as the program will stop in this case.
-
-        if ((nicore.eq.'pche').or.(nicore.eq.'fche')) then
-          call vhrtre(chcore,ve,rofi,drdi,s,nrval,a)
-          do l=0,lmax
-            do ir=2,nrval
-              vps(ir,l)=vps(ir,l)+ve(ir)
-            enddo
-            vps(1,l) = vps(2,l)    ! AG
-          enddo
-        endif
-
-        return
-
- 5000   continue
-        write(6,*)
-     .      'ERROR: You are using an old pseudopotential file.',
-     .      ' Siesta needs a newer version.'
-        call die
-          
+
         end subroutine read_vps
 !
-        subroutine comKB(is,a,b,rofi,proj,l,ikb,rc,ekb,nrc,spp)
+        subroutine comKB(is,a,b,rofi,proj,l,lj_projs,jk,
+     $                   ikb,rc,ekb,nrc,spp)
 C
 C  Creates the common block with all the information about the 
 C  Kleinman-Bylander projectors.
@@ -2906,6 +2916,8 @@
         integer, intent(in)  :: l, nrc,is, ikb 
         real(dp), intent(in) :: rc, ekb, proj(nrmax), a, b,
      .       rofi(nrmax)
+        logical, intent(in)  :: lj_projs
+        integer, intent(in)  :: jk
         type(species_info), intent(inout)  :: spp
         
         character(len=40) filename
@@ -2932,10 +2944,17 @@
         
         spp%pjnl_n(n) = ikb
         spp%pjnl_l(n) = l
+        if (lj_projs) then
+           spp%pjnl_j(n) = l + (2*jk-3)*0.5_dp ! l -/+ 1/2
+           if (l == 0 ) spp%pjnl_j(n) = 0.5_dp ! special case: only jk=1
+        else
+           spp%pjnl_j(n) = 0.0_dp
+        endif
         do m = -l, l
            ntot = ntot + 1
            spp%pj_index(ntot) = n
            spp%pj_l(ntot) = l
+           spp%pj_j(ntot) = spp%pjnl_j(n)
            spp%pj_m(ntot) = m
         enddo
         spp%nprojs = ntot
@@ -2983,8 +3002,8 @@
        end subroutine comkb
 !
         subroutine KBgen(is, a,b,rofi,drdi,s, 
-     .         vps, vlocal, ve, nrval, Zval, lmxkb, 
-     .         nkbl, erefkb, nkb, spp)
+     .         vps, vps_u, vlocal, ve, nrval, Zval, lmxkb, 
+     .         nkbl, erefkb, nkb, spp, lj_projs)
 
         use basis_specs, only: restricted_grid
         use atom_options, only: new_kb_reference_orbitals
@@ -3002,28 +3021,34 @@
 
         implicit none
 
-        real(dp) 
-     .    a, b, rofi(nrmax), vps(nrmax,0:lmaxd),
+        real(dp), intent(in) ::
+     .    a, b, rofi(nrmax), vps(nrmax,0:lmaxd), vps_u(nrmax,0:lmaxd),
      .    drdi(nrmax), s(nrmax), ve(nrmax),vlocal(nrmax),
      .    Zval, erefkb(nkbmx,0:lmaxd)
 
-        integer
+        integer, intent(in) ::
      .    nrval, lmxkb, nkb, is, nkbl(0:lmaxd)
 
         type(species_info), intent(inout)  :: spp
+        real(dp), allocatable :: rphi2(:,:,:), vii(:,:)
+        real(dp), allocatable :: rphi(:,:,:)
+        real(dp), allocatable :: eref(:,:)
+
+        logical, intent(in) :: lj_projs
 C
 C Internal variables
 C
         integer 
      .    l,nprin, nnodes, ighost, nrwf, ikb, ir,
-     .    nrc, nrlimit, n_pjnl, nkbs_tot
-        real(dp)
-     .    rc(nkbmx,0:lmaxd), dkbcos(nkbmx,0:lmaxd),
-     .    ekb(nkbmx,0:lmaxd)
+     .    nrc, nrlimit, n_pjnl, nkbs_tot, jk, nj
+        real(dp) :: rc, dkbcos, ekb
                
         real(dp)
-     .    rphi(nrmax,nkbmx), rmax, dnrm, 
-     .    proj(nrmax)
+     .    rmax, dnrm, 
+     .    proj(nrmax), vp(nrmax)
+
+        character(len=2) :: jstr
+        logical :: multiple_projectors
                  
 C  The atomic wavefunctions and/or its energy derivatives are
 C  calculated only inside a sphere of radius Rmax. To define the
@@ -3056,13 +3081,20 @@
 
         spp%lmax_projs = lmxkb
         
-        ! Generalize this for lj projectors
         nkbs_tot = 0
         n_pjnl = 0
         do l = 0, lmxkb
-           n_pjnl = n_pjnl + nkbl(l)
-           nkbs_tot = nkbs_tot + nkbl(l)*(2*l+1)
+           do ikb = 1, nkbl(l)
+              n_pjnl = n_pjnl + 1
+              nkbs_tot = nkbs_tot + (2*l+1)
+              if (lj_projs .and. l>0) then  ! j=l+1/2 and j=l-1/2 shells
+                 n_pjnl = n_pjnl + 1
+                 nkbs_tot = nkbs_tot + (2*l+1)
+              endif
+           enddo
         enddo
+
+        spp%lj_projs = lj_projs
         
         ! Eventually will make the arrays allocatable
         if (n_pjnl > maxn_pjnl) then
@@ -3076,88 +3108,128 @@
         ! zero out for build-up in comkb
         spp%n_pjnl = 0   
         spp%nprojs = 0
+
+        write(6,'(/,a)')'KBgen: Kleinman-Bylander projectors: '
+
+        multiple_projectors = .false.
+
+        if (lj_projs) then
+           nj = 2
+        else
+           nj = 1
+        endif
         
+        ! Some output related to ghosts might be changed
+        ! Note ordering of projectors
+        ! Independent projector stacks for different j's
+        allocate (rphi2(nrmax,nkbmx,nj), vii(nkbmx,nj))
+        ! For the derivative option, we need a previous projector 
+        ! and the previous reference energy
+        allocate (rphi(nrmax,nkbmx,nj))
+        allocate (eref(nkbmx,nj))
         do l=0,lmxkb
-          do ir=1,nrmax
-            do ikb=1,nkbmx
-              rphi(ir,ikb)=0.0d0
-            enddo 
-            proj(ir)=0.0d0
-          enddo
+          rphi(:,:,:) = 0.0_dp
           do ikb= 1, nkbl(l)
-C
-C Atomic wavefunctions and eigenvalues
-C for the construction of the KB projectors
-C
-C If the reference energies have not been specifed, the eigenstates
-C with the condition of being zero at r(nrval) will be used.
-C       
+            do jk = 1, nj      ! j-, j+
+                
+             if (jk == 2 .and. l==0) CYCLE  ! only one proj for l=0
+             if (lj_projs .and. l>0 ) then
+              ! Set up the appropriate potential and label
+              if (jk == 1) then
+                 ! V(l-j/2) = Vdn - (l+1)/2 Vup = Vion - (l+1)/2 Vso
+                 vp(:) = vps(:,l) - 0.5_dp*(l+1)*vps_u(:,l)
+                 jstr = 'j-'
+              else
+                 ! V(l+1/2) = Vdn + l/2 Vup = Vion + l/2 Vso
+                 vp(:) = vps(:,l) + 0.5_dp*l*vps_u(:,l)
+                 jstr = 'j+'
+              endif
+             else
+              vp(:) = vps(:,l)
+              jstr = '  '
+             endif
+             proj(:) = 0.0_dp
+           
+! Atomic wavefunctions and eigenvalues
+! for the construction of the KB projectors
+
+           ! We use the same user-entered reference energies for both j's
+
            if (erefkb(ikb,l).ge.1.0d3) then             
+             ! If the reference energies have not been specifed, the
+             ! eigenstates with the condition of being zero at r(nrval)
+             ! will be used.
              nnodes=ikb
              nprin=l+1
-
-             call schro_eq(Zval,rofi,vps(1,l),ve,s,drdi,
+             ! use eref to avoid overwriting erefkb(,) for next j
+             call schro_eq(Zval,rofi,vp,ve,s,drdi,
      .                     nrlimit,l,a,b,nnodes,nprin,
-     .                     erefkb(ikb,l),rphi(1,ikb)) 
+     .                     eref(ikb,jk),rphi(1,ikb,jk)) 
 C Normalization of the eigenstates inside a sphere of radius Rmax
              dnrm=0.0d0
              do ir=1,nrwf
-               dnrm=dnrm+drdi(ir)*rphi(ir,ikb)**2
+               dnrm=dnrm+drdi(ir)*rphi(ir,ikb,jk)**2
              enddo 
              dnrm=sqrt(dnrm)
              do ir=1,nrwf
-               rphi(ir,ikb)=rphi(ir,ikb)/dnrm
+               rphi(ir,ikb,jk)=rphi(ir,ikb,jk)/dnrm
              enddo 
 C
            elseif((erefkb(ikb,l).le.-1.0d3).and.
      .       (ikb.gt.1) ) then 
-C If the energy is specified to be 1000 Ry, the energy derivative
-C of the previous wavefunction will be used
-C
-             call energ_deriv(a,rofi,rphi(1,ikb-1),vps(1,l),
-     .                        ve,drdi,nrwf,l,erefkb(ikb-1,l),
-     .                        rphi(1,ikb),nrval)
-             erefkb(ikb,l)=0.0d0
-C
+             ! If the energy is specified to be 1000 Ry, the energy
+             ! derivative of the previous wavefunction will be used
+             call energ_deriv(a,rofi,rphi(1,ikb-1,jk),vp(1),
+     .                        ve,drdi,nrwf,l,eref(ikb-1,jk),
+     .                        rphi(1,ikb,jk),nrval)
+             eref(ikb,jk)=0.0_dp
+
            else 
-C If the reference energies have been specified, we just use them
-C
-             call rphi_vs_e(a,b,rofi,vps(1,l),
+             ! If the reference energies have been specified, we just
+             ! use them
+             call rphi_vs_e(a,b,rofi,vp(1),
      .                      ve,nrval,l,erefkb(ikb,l),
-     .                      rphi(1,ikb),Rmax)
-C
+     .                      rphi(1,ikb,jk),Rmax)
+             eref(ikb,jk) = erefkb(ikb,l)
            endif 
-C
-C Ghost analysis
-C 
+
+           ! Ghost analysis
+
            if (nkbl(l).eq.1) then
-             call ghost(Zval,rofi,vps(:,l),vlocal,
+             call ghost(Zval,rofi,vp(:),vlocal,
      .                  ve,s,drdi,nrlimit,l,a,b,nrwf,
-     .                  erefkb(ikb,l),rphi(:,ikb),ighost)
+     .                  eref(ikb,jk),rphi(:,ikb,jk),ighost)
            else 
-             if (ikb.eq.1)
-     .         write(6,'(a,i3,/a)') 
-     .         'KBgen: More than one KB projector for l=',l,
-     .         'KBgen: ghost states analysis will be not performed'
+              multiple_projectors = .true.
            endif
-C
-C KB Projectors
-C
-           call KBproj(ikb,rofi,drdi,vps(1,l),vlocal,nrwf,l,
-     .                rphi(1,ikb),dkbcos(ikb,l),ekb(ikb,l),proj,nrc)
-
-C
-           rc(ikb,l)=rofi(nrc)
-C
-C Common block with the information about the  KB projectors
-C
-          call comKB(is,a,b,rofi,proj,
-     .                 l,ikb,rc(ikb,l),ekb(ikb,l),nrc, spp)
-C
-
-          enddo 
-         enddo   
+
+           ! The projector stack and vii are passed as arguments
+           ! to enable the needed loop ordering
+           call KBproj(ikb,rofi,drdi,vp,vlocal,nrwf,l,
+     .                rphi(1,ikb,jk),dkbcos,ekb,proj,nrc,
+     .                rphi2(:,:,jk),vii(:,jk))
+
+           rc = rofi(nrc)
+
+          ! Store KB projectors in data structure
+
+          call comKB(is,a,b,rofi,proj,l,lj_projs,jk,
+     .               ikb,rc,ekb,nrc, spp)
+
+          write(6,'(a2,1x,a,i2,4(3x,a,f10.6))') jstr,
+     .        'l=',l, 'rc=',rc, 'el=',eref(ikb,jk), 
+     .        'Ekb=',ekb,'kbcos=',dkbcos
+
+          enddo ! jk
+         enddo  ! ikb:1,nkbl
+        enddo   ! l
+        deallocate(rphi,rphi2,vii,eref)
         
+        if (multiple_projectors) then
+         write(6,'(a,/a)') 
+     .  'KBgen: More than one KB projector for some l shell(s)',
+     .  'KBgen: ghost-state analysis will not be performed'
+        endif
          if (ighost.eq.1) then
             write(6,"(2a)")'KBgen: WARNING: ',
      .            'Ghost states have been detected'
@@ -3169,34 +3241,14 @@
                write(6,"(a)") " KBgen: *** Warning Ignored by User"
                ighost = 0
             else
-               call die
+               call die("Ghost states detected")
             endif
          endif
 
-         write(6,'(/,a)')'KBgen: Kleinman-Bylander projectors: '
-         do l=0,lmxkb
-           do ikb=1, nkbl(l)
-              write(6,'(3x,a,i2,4(3x,a,f10.6))')
-     .        'l=',l, 'rc=',rc(ikb,l), 'el=',erefkb(ikb,l), 
-     .        'Ekb=',ekb(ikb,l),'kbcos=',dkbcos(ikb,l)
-           enddo
-         enddo
-C
-C Total number of Kleinman-Bylander projectors
-C
-           nkb=0
-           do l=0,lmxkb
-              do ikb=1,nkbl(l)
-                 nkb=nkb+(2*l+1) 
-              enddo 
-           enddo 
-           write(6,'(/,a, i4)')
-     .'KBgen: Total number of  Kleinman-Bylander projectors: ', nkb
-           if (nkb /= spp%nprojs) then
-              print *, "spp%nprojs: ", spp%nprojs
-              call die("Mismatch nkb spp%nprojs")
-           endif
-C     
+         write(6,'(/,a, i4)')
+     .'KBgen: Total number of  Kleinman-Bylander projectors: ',
+     $          spp%nprojs
+
         end subroutine KBgen
 !
         subroutine Basis_gen(Zval,is, iz, a,b,rofi,drdi,s,

=== modified file 'Src/bands.F'
--- Src/bands.F	2018-04-25 11:33:32 +0000
+++ Src/bands.F	2018-04-30 20:53:35 +0000
@@ -366,7 +366,7 @@
       use atmfuncs,     only : symfio, cnfigfio, labelfis, nofis
       use writewave,    only : wfs_filename
 
-      use m_spin,       only: NoMagn, SPpol, NonCol, SpOrb
+      use m_spin,       only: spin 
 
       use m_diag_option, only: ParallelOverK, Serial
       use m_diag, only: diag_init
@@ -459,7 +459,7 @@
       endif
 
 C Find the band energies
-      if (NoMagn .or. SPpol) then
+      if ( spin%none .or. spin%Col ) then
 C fixspin and qs are not used in diagk, since getD=.false. ...
         qs(1) = 0.0_dp
         qs(2) = 0.0_dp
@@ -488,7 +488,7 @@
         ParallelOverK = SaveParallelOverK
         if ( ParallelOverK ) Serial = .true.
 
-      elseif ( NonCol ) then
+      elseif ( spin%NCol ) then
          if (getPSI) then
             if (node==0) then
                write(6,*) "No WFS for non-colinear spin, yet..."
@@ -503,7 +503,7 @@
      .              Haux, Saux, psi, Haux, Saux, aux,
      .              no_u, occtol, 1, no_u )
 
-      elseif ( SpOrb ) then
+      elseif ( spin%SO ) then
         if (getPSI) then
            if (node==0) then
               write(6,*) "No WFS for spin-orbit, yet..."

=== modified file 'Src/basis_io.F'
--- Src/basis_io.F	2017-12-22 10:26:20 +0000
+++ Src/basis_io.F	2018-04-30 20:53:35 +0000
@@ -401,8 +401,14 @@
 ! KBs
          read(lun,*)
          do i=1,spp%n_pjnl
-            read(lun,*)
+            if (spp%lj_projs) then
+               read(lun,*)
+     $           spp%pjnl_l(i), spp%pjnl_j(i), spp%pjnl_n(i),
+     $           spp%pjnl_ekb(i)
+            else
+               read(lun,*)
      $           spp%pjnl_l(i), spp%pjnl_n(i), spp%pjnl_ekb(i)
+            endif
             call radial_read_ascii(spp%pjnl(i),lun,
      $                             yp1=0.0_dp,ypn=huge(1.0_dp))
          enddo
@@ -417,6 +423,7 @@
                nk = nk+1
                spp%pj_n(nk) = spp%pjnl_n(i)
                spp%pj_l(nk) = spp%pjnl_l(i)
+               spp%pj_j(nk) = spp%pjnl_j(i)
                spp%pj_m(nk) = m
                spp%pj_index(nk) = i
             enddo
@@ -450,6 +457,7 @@
       integer, intent(in)         :: unit
       
       character(len=78) line
+      integer :: iostat
 
       read(unit,'(a)') line
       if (trim(line) .eq. '<preamble>') then
@@ -465,7 +473,12 @@
       read(unit,*) p%mass
       read(unit,*) p%self_energy
       read(unit,*) p%lmax_basis, p%n_orbnl
-      read(unit,*) p%lmax_projs, p%n_pjnl
+      read(unit,fmt=*,iostat=iostat) p%lmax_projs, p%n_pjnl, p%lj_projs
+      if (iostat /=0 ) then
+         backspace(unit)
+         read(unit,*) p%lmax_projs, p%n_pjnl
+         p%lj_projs = .false.
+      endif
 
       allocate(p%orbnl(p%n_orbnl))
       allocate(p%pjnl(p%n_pjnl))
@@ -817,9 +830,16 @@
          do i=1,spp%n_pjnl
             zeta = spp%pjnl_n(i)
             l = spp%pjnl_l(i)
-            write(lun,'(2i3,f22.16,2x,a)')
+            if (spp%lj_projs) then
+               write(lun,'(i3,f4.1,i3,f22.16,2x,a)')
+     $           spp%pjnl_l(i), spp%pjnl_j(i), spp%pjnl_n(i),
+     $           spp%pjnl_ekb(i),
+     $           " #kb l, j, n (sequence number), Reference energy"
+            else
+               write(lun,'(2i3,f22.16,2x,a)')
      $           spp%pjnl_l(i), spp%pjnl_n(i), spp%pjnl_ekb(i),
      $           " #kb l, n (sequence number), Reference energy"
+            endif
             call radial_dump_ascii(spp%pjnl(i),lun)
 
            if (write_ion_plot_files) then
@@ -949,8 +969,14 @@
          write(unit,'(g22.12,4x,a)') p%self_energy, "# Self energy "
          write(unit,'(2i4,22x,a)') p%lmax_basis, p%n_orbnl,
      $        "# Lmax for basis, no. of nl orbitals "
-         write(unit,'(2i4,22x,a)') p%lmax_projs, p%n_pjnl,
-     $        "# Lmax for projectors, no. of nl KB projectors "
+         if (p%lj_projs) then
+            write(unit,'(2i4,1x,l1,20x,a)') p%lmax_projs, p%n_pjnl,
+     $       p%lj_projs,
+     $       "# Lmax for projectors, no. of nl KB projectors, LJ projs?"
+         else
+            write(unit,'(2i4,22x,a)') p%lmax_projs, p%n_pjnl,
+     $           "# Lmax for projectors, no. of nl KB projectors"
+         endif
 
          end subroutine write_header
 

=== modified file 'Src/broadcast_basis.F'
--- Src/broadcast_basis.F	2016-06-28 05:57:03 +0000
+++ Src/broadcast_basis.F	2018-04-30 20:53:35 +0000
@@ -104,10 +104,14 @@
      $                 0,MPI_Comm_World,MPIerror)
         call MPI_Bcast(spp%lmax_projs,1,MPI_integer,
      $                 0,MPI_Comm_World,MPIerror)
+        call MPI_Bcast(spp%lj_projs,1,MPI_logical,
+     $                 0,MPI_Comm_World,MPIerror)
         call MPI_Bcast(spp%pjnl_l,maxn_pjnl,MPI_integer,
      $                 0,MPI_Comm_World,MPIerror)
         call MPI_Bcast(spp%pjnl_n,maxn_pjnl,MPI_integer,
      $                 0,MPI_Comm_World,MPIerror)
+        call MPI_Bcast(spp%pjnl_j,maxn_pjnl,MPI_double_precision,
+     $                 0,MPI_Comm_World,MPIerror)
         call MPI_Bcast(spp%pjnl_ekb,maxn_pjnl,MPI_double_precision,
      $                 0,MPI_Comm_World,MPIerror)
 
@@ -134,6 +138,8 @@
      $                 0,MPI_Comm_World,MPIerror)
         call MPI_Bcast(spp%pj_l,maxnprojs,MPI_integer,
      $                 0,MPI_Comm_World,MPIerror)
+        call MPI_Bcast(spp%pj_j,maxnprojs,MPI_double_precision,
+     $                 0,MPI_Comm_World,MPIerror)
         call MPI_Bcast(spp%pj_m,maxnprojs,MPI_integer,
      $                 0,MPI_Comm_World,MPIerror)
 

=== modified file 'Src/compute_energies.F90'
--- Src/compute_energies.F90	2017-12-21 15:49:49 +0000
+++ Src/compute_energies.F90	2018-04-30 20:53:35 +0000
@@ -47,8 +47,10 @@
                                  lastkb, no_s, rmaxv, indxua, iphorb, lasto, &
                                  rmaxo, no_l
       use m_ntm,           only: ntm
-      use m_spin,          only: NoMagn, SPpol, NonCol, SpOrb 
-      use m_spin,          only: nspin, h_spin_dim, spinor_dim
+
+      use m_spin,          only: spin
+      use parallel,        only: IONode 
+
       use m_dipol,         only: dipol
       use siesta_geom,     only: na_u, na_s, xa, isa
       use m_rhog,          only: rhog
@@ -58,22 +60,20 @@
 
       integer, intent(in)   :: iscf
 
-      integer               :: ihmat, ifa, istr, ispin, io
+      integer               :: ispin, io
 #ifdef MPI
       real(dp) :: buffer1
 #endif
 
-      real(dp) :: const
-      real(dp) :: dummy_stress(3,3), dummy_fa(1,1)
-      real(dp) :: dummy_E, g2max, dummy_H(1,1)
       logical  :: mixDM
 
+
       mixDM = (.not. (mixH .or. mix_charge))
 
 !     Compute the band-structure energy
 
       call compute_EBS()
-    
+
       ! These energies were calculated in the latest call to
       ! setup_hamiltonian, using as ingredient D_in 
 
@@ -126,36 +126,89 @@
 
    subroutine compute_EBS()
 
+     real(dp)    :: Ebs_SO(4)
+     complex(dp) :: Ebs_Daux(2,2), Ebs_Haux(2,2)
+
+      integer :: i, j, ind
+
       Ebs = 0.0_dp
-!       const factor takes into account that there are two nondiagonal
-!       elements in non-collinear spin density matrix, stored as
-!       ispin=1 => D11; ispin=2 => D22, ispin=3 => Real(D12);
-!       ispin=4 => Imag(D12)
-
-      if ( SpOrb ) then
-        do io = 1,maxnh
-          Ebs = Ebs      + H(io,1) * ( Dscf(io,1) )   &
-                         + H(io,2) * ( Dscf(io,2) )   &
-                         + H(io,3) * ( Dscf(io,7) )   &
-                         + H(io,4) * ( Dscf(io,8) )   &
-                         - H(io,5) * ( Dscf(io,5) )   &
-                         - H(io,6) * ( Dscf(io,6) )   &
-                         + H(io,7) * ( Dscf(io,3) )   &
-                         + H(io,8) * ( Dscf(io,4) )
+
+! Modifed for Off-Site Spin-orbit coupling by R. Cuadrado, Feb. 2018
+!
+!*****************************************************************************
+!  Note about Ebs and E_Harris calculation for the Spin-Orbit:
+!*****************************************************************************
+!
+!  E_bs and E_Harris are calculated by means of the following 
+! complex matrix multiplication: 
+!
+!                 E_bs = Re { Tr[ H * DM ] } 
+!             E_Harris = Re { Tr[ H * (DM-DM_old) ] }  
+!
+!  In the following: DM/H(1,1) --> up/up     <--> uu
+!                    DM/H(2,2) --> down/down <--> dd
+!                    DM/H(1,2) --> up/down   <--> ud
+!                    DM/H(2,1) --> down/up   <--> du
+!
+!  Using DM/H components, E_bs would be sum(E_bs(1:4)), where
+!
+!   E_bs(1)=Re{sum_ij(H_ij(1,1)*D_ji(1,1))}=Re{sum_ij(H_ij^uu*(DM_ij^uu)^*)}
+!   E_bs(2)=Re{sum_ij(H_ij(2,2)*D_ji(2,2))}=Re{sum_ij(H_ij^dd*(DM_ij^dd)^*)}
+!   E_bs(3)=Re{sum_ij(H_ij(1,2)*D_ji(2,1))}=Re{sum_ij(H_ij^ud*(DM_ij^ud)^*)}
+!   E_bs(4)=Re{sum_ij(H_ij(2,1)*D_ji(1,2))}=Re{sum_ij(H_ij^du*(DM_ij^du)^*)}
+!
+!         since, due to overall hermiticity,  DM_ij^ab = (DM_ji^ba)^*
+!
+!   The trace operation is then an extended dot product over the "ij"
+!   sparse index, which can also be conveniently done in parallel, as
+!   each processor handles the same indexes in H and the DM. Only a
+!   global reduction is needed at the end.
+      
+!  Same comments are valid for the E_Harris calculation.
+!
+!*****************************************************************************
+!
+      if ( spin%SO ) then
+
+        Ebs_SO = 0.0_dp  
+        Ebs_Daux = dcmplx(0.0_dp, 0.0_dp)
+        Ebs_Haux = dcmplx(0.0_dp, 0.0_dp)
+
+        do io = 1, maxnh
+
+          Ebs_Haux(1,1) = dcmplx(H(io,1),H(io,5))  
+          Ebs_Haux(2,2) = dcmplx(H(io,2),H(io,6))  
+          Ebs_Haux(1,2) = dcmplx(H(io,3),-H(io,4)) 
+          Ebs_Haux(2,1) = dcmplx(H(io,7),H(io,8)) 
+
+          Ebs_Daux(1,1) = dcmplx(Dscf(io,1),Dscf(io,5)) 
+          Ebs_Daux(2,2) = dcmplx(Dscf(io,2),Dscf(io,6)) 
+          Ebs_Daux(1,2) = dcmplx(Dscf(io,3),-Dscf(io,4))
+          Ebs_Daux(2,1) = dcmplx(Dscf(io,7),Dscf(io,8)) 
+
+
+          Ebs_SO(1) = Ebs_SO(1) + real( Ebs_Haux(1,1)*dconjg(Ebs_Daux(1,1)) )
+          Ebs_SO(2) = Ebs_SO(2) + real( Ebs_Haux(2,2)*dconjg(Ebs_Daux(2,2)) )
+          Ebs_SO(3) = Ebs_SO(3) + real( Ebs_Haux(1,2)*dconjg(Ebs_Daux(1,2)) )
+          Ebs_SO(4) = Ebs_SO(4) + real( Ebs_Haux(2,1)*dconjg(Ebs_Daux(2,1)) )
+
         enddo
-      else if ( NonCol ) then
+
+        Ebs = sum ( Ebs_SO )
+
+      else if ( spin%NCol ) then
         do io = 1,maxnh
           Ebs    = Ebs    + H(io,1) * ( Dscf(io,1)  )   &
                           + H(io,2) * ( Dscf(io,2)  )   &
                  + 2.0_dp * H(io,3) * ( Dscf(io,3)  )   &
                  + 2.0_dp * H(io,4) * ( Dscf(io,4)  )
         enddo
-      else if ( SPpol )  then
+      else if ( spin%Col )  then
         do io = 1,maxnh
           Ebs    = Ebs    + H(io,1) * Dscf(io,1)  &
                           + H(io,2) * Dscf(io,2)
         enddo
-      else if ( NoMagn ) then
+      else if ( spin%none ) then
         do io = 1,maxnh
           Ebs    = Ebs + H(io,1) *  Dscf(io,1)
         enddo
@@ -170,36 +223,65 @@
 
     subroutine compute_DEharr()
 
+      real(dp)    :: DEharr_SO(4)
+      complex(dp) :: DEharr_Daux(2,2), DEharr_Haux(2,2), DEharr_Daux_old(2,2)
+
       DEharr = 0.0_dp
-             !       const factor takes into account that there are two nondiagonal
-             !       elements in non-collinear spin density matrix, stored as
-             !       ispin=1 => D11; ispin=2 => D22, ispin=3 => Real(D12);
-             !       ispin=4 => Imag(D12)
-
-      if ( SpOrb ) then
-        do io = 1,maxnh
-          DEharr = DEharr + H(io,1) * ( Dscf(io,1) - Dold(io,1) )  &
-                          + H(io,2) * ( Dscf(io,2) - Dold(io,2) )  &
-                          + H(io,3) * ( Dscf(io,7) - Dold(io,7) )  &
-                          + H(io,4) * ( Dscf(io,8) - Dold(io,8) )  &
-                          - H(io,5) * ( Dscf(io,5) - Dold(io,5) )  &
-                          - H(io,6) * ( Dscf(io,6) - Dold(io,6) )  &
-                          + H(io,7) * ( Dscf(io,3) - Dold(io,3) )  &
-                          + H(io,8) * ( Dscf(io,4) - Dold(io,4) )
-        enddo
-      else if ( NonCol ) then
+
+      if ( spin%SO ) then
+
+        DEharr_SO = 0.0_dp 
+        DEharr_Daux     = dcmplx(0.0_dp, 0.0_dp)
+        DEharr_Daux_old = dcmplx(0.0_dp, 0.0_dp)
+        DEharr_Haux     = dcmplx(0.0_dp, 0.0_dp)
+
+        do io = 1, maxnh
+
+          DEharr_Haux(1,1) = dcmplx( H(io,1),H(io,5) )
+          DEharr_Haux(2,2) = dcmplx( H(io,2),H(io,6) )
+          DEharr_Haux(1,2) = dcmplx( H(io,3),-H(io,4) )
+          DEharr_Haux(2,1) = dcmplx( H(io,7),H(io,8) )
+
+          DEharr_Daux(1,1) = dcmplx( Dscf(io,1),Dscf(io,5) )
+          DEharr_Daux(2,2) = dcmplx( Dscf(io,2),Dscf(io,6) )
+          DEharr_Daux(1,2) = dcmplx( Dscf(io,3),-Dscf(io,4) )
+          DEharr_Daux(2,1) = dcmplx( Dscf(io,7),Dscf(io,8) )
+
+          DEharr_Daux_old(1,1) = dcmplx( Dold(io,1),Dold(io,5) )
+          DEharr_Daux_old(2,2) = dcmplx( Dold(io,2),Dold(io,6) ) 
+          DEharr_Daux_old(1,2) = dcmplx( Dold(io,3),-Dold(io,4) ) 
+          DEharr_Daux_old(2,1) = dcmplx( Dold(io,7),Dold(io,8) )
+
+
+          DEharr_SO(1) = DEharr_SO(1) &
+             + real( DEharr_Haux(1,1)*dconjg(DEharr_Daux(1,1)-DEharr_Daux_old(1,1)) ) 
+
+          DEharr_SO(2) = DEharr_SO(2) &
+             + real( DEharr_Haux(2,2)*dconjg(DEharr_Daux(2,2)-DEharr_Daux_old(2,2)) )
+
+          DEharr_SO(3) = DEharr_SO(3) &
+             + real( DEharr_Haux(1,2)*dconjg(DEharr_Daux(1,2)-DEharr_Daux_old(1,2)) )
+
+          DEharr_SO(4) = DEharr_SO(4) &
+             + real( DEharr_Haux(2,1)*dconjg(DEharr_Daux(2,1)-DEharr_Daux_old(2,1)) )
+
+         enddo
+
+         DEharr = sum ( DEharr_SO )
+
+      else if ( spin%NCol ) then
         do io = 1,maxnh
           DEharr = DEharr + H(io,1) * ( Dscf(io,1) - Dold(io,1) )  &
                           + H(io,2) * ( Dscf(io,2) - Dold(io,2) )  &
                  + 2.0_dp * H(io,3) * ( Dscf(io,3) - Dold(io,3) )  &
                  + 2.0_dp * H(io,4) * ( Dscf(io,4) - Dold(io,4) )
         enddo
-      elseif (SPpol)  then
+      elseif ( spin%Col )  then
         do io = 1,maxnh
           DEharr = DEharr + H(io,1) * ( Dscf(io,1) - Dold(io,1) )  &
                           + H(io,2) * ( Dscf(io,2) - Dold(io,2) )
         enddo
-      elseif (NoMagn) then
+      elseif ( spin%none ) then
         do io = 1,maxnh
           DEharr = DEharr + H(io,1) * ( Dscf(io,1) - Dold(io,1) )
         enddo
@@ -217,11 +299,21 @@
       use files, only : filesOut_t    ! derived type for output file names
       use class_dSpData1D, only : val
       use class_dSpData2D, only : val
-      use sparse_matrices, only: H_kin_1D, H_vkb_1D, H_so_2D
+      use class_zSpData2D, only : val
+      use sparse_matrices, only: H_kin_1D, H_vkb_1D
+      use sparse_matrices, only: H_so_on_2D, H_so_off_2D
+
 
       type(filesOut_t)  :: filesOut  ! blank output file names
-      real(dp), pointer :: H_vkb(:), H_kin(:), H_so(:,:)
-
+      real(dp), pointer :: H_vkb(:), H_kin(:), H_so_on(:,:)
+      complex(dp), pointer :: H_so_off(:,:)
+
+
+      complex(dp) :: Hc, Dc
+      real(dp)    :: dummy_stress(3,3), dummy_fa(1,1)
+      real(dp)    :: dummy_E, g2max, dummy_H(1,1)
+      integer     :: ihmat, ifa, istr
+      
       ! Compute E_KS(DM_out)
 
       g2max = g2cut
@@ -233,7 +325,7 @@
 
       ! Remove unwanted arguments...
 
-      call dhscf( nspin, no_s, iaorb, iphorb, no_l,      &
+      call dhscf( spin%Grid, no_s, iaorb, iphorb, no_l,      &
                   no_u, na_u, na_s, isa, xa, indxua,                        &
                   ntm, ifa, istr, ihmat, filesOut,                          &
                   maxnh, numh, listhptr, listh, Dscf, Datm,                 &
@@ -250,7 +342,7 @@
       
       Ekin = 0.0_dp
       Enl  = 0.0_dp
-      do ispin = 1, spinor_dim
+      do ispin = 1, spin%spinor
          do io = 1,maxnh
             Ekin = Ekin + H_kin(io) * Dscf(io,ispin)
             Enl  = Enl  + H_vkb(io) * Dscf(io,ispin)
@@ -266,21 +358,48 @@
 #endif
 
       Eso = 0._dp
-      if ( SpOrb ) then
+
+      if ( spin%SO_offsite ) then
+         H_so_off => val(H_so_off_2D)
+
+         ! The computation of the trace is different here, as H_so_off has
+         ! a different structure from H and the DM.
+        do io = 1, maxnh
+
+!-------- Eso(u,u)
+          Dc = cmplx(Dscf(io,1),Dscf(io,5), dp)
+          Eso = Eso + real( H_so_off(io,1)*Dc, dp)
+!-------- Eso(d,d)
+          Dc = cmplx(Dscf(io,2),Dscf(io,6), dp)
+          Eso = Eso + real( H_so_off(io,2)*Dc, dp)
+!-------- Eso(u,d)
+          Dc = cmplx(Dscf(io,3),Dscf(io,4), dp)
+          Eso = Eso + real( H_so_off(io,4)*Dc, dp)
+!-------- Eso(d,u)
+          Dc = cmplx(Dscf(io,7),-Dscf(io,8), dp)
+          Eso = Eso + real( H_so_off(io,3)*Dc, dp)
+
+        end do
+
+      else if ( spin%SO_onsite ) then
          ! Sadly some compilers (g95), does
          ! not allow bounds for pointer assignments :(
-         H_so => val(H_so_2D)
+         H_so_on => val(H_so_on_2D)
          do io = 1,maxnh
-            Eso = Eso + H_so(io,1)*Dscf(io,7) + H_so(io,2)*Dscf(io,8) &
-                 + H_so(io,5)*Dscf(io,3) + H_so(io,6)*Dscf(io,4) &
-                 - H_so(io,3)*Dscf(io,5) - H_so(io,4)*Dscf(io,6)
-         end do
+            Eso = Eso + H_so_on(io,1)*Dscf(io,7) + H_so_on(io,2)*Dscf(io,8) &
+                 + H_so_on(io,5)*Dscf(io,3) + H_so_on(io,6)*Dscf(io,4) &
+                 - H_so_on(io,3)*Dscf(io,5) - H_so_on(io,4)*Dscf(io,6)
+          end do
+          
+      end if
+
 #ifdef MPI
+      if (spin%SO) then
          ! Global reduction of Eso
          call globalize_sum( Eso, buffer1 )
          Eso = buffer1
+      endif
 #endif
-      end if
       
       ! E0 = Ena + Ekin + Enl + Eso - Eions
 

=== modified file 'Src/dfscf.f'
--- Src/dfscf.f	2018-04-22 00:13:43 +0000
+++ Src/dfscf.f	2018-04-30 20:53:35 +0000
@@ -20,6 +20,7 @@
 C Adds the SCF contribution to atomic forces and stress.
 C Written by P.Ordejon, J.M.Soler, and J.Gale.
 C Last modification by J.M.Soler, October 2000.
+C Modifed for Off-Site Spin-orbit coupling by R. Cuadrado, Feb. 2018
 C *********************** INPUT **************************************
 C integer ifa             : Are forces required? (1=yes,0=no)
 C integer istr            : Is stress required? (1=yes,0=no)
@@ -334,7 +335,7 @@
 
 C     Factor two for nondiagonal elements for non-collinear spin (and SO)
         if ( nspin == 4 ) then
-           V(1:nsp,3:4) = 2.0_dp * V(1:nsp,3:4)
+         V(1:nsp,3:4) = 2.0_dp * V(1:nsp,3:4)
         end if
 
 C     Loop on first orbital of mesh point

=== modified file 'Src/final_H_f_stress.F'
--- Src/final_H_f_stress.F	2017-11-23 14:17:27 +0000
+++ Src/final_H_f_stress.F	2018-04-30 20:53:35 +0000
@@ -21,8 +21,9 @@
       use siesta_options, only: recompute_H_after_scf
       use sparse_matrices, only: numh, listh, listhptr
       use sparse_matrices, only: H, S, Dscf, Escf, maxnh, xijo
-      use sparse_matrices, only: H_ldau_2D
+      use sparse_matrices, only: H_ldau_2D, H_so_off_2D
       use class_dSpData2D, only: val
+      use class_zSpData2D, only: val
 
       use siesta_geom
       use atomlist, only: no_u, iaorb, iphkb, qtot, indxuo, datm, 
@@ -30,7 +31,7 @@
      &                    rmaxo, no_l, iza
       use metaforce, only: lMetaForce, meta
       use molecularmechanics, only : twobody
-      use m_nlefsm,     only: nlefsm
+      use m_nlefsm,     only: nlefsm, nlefsm_SO_off
       use m_overfsm,    only: overfsm
       use m_kinefsm,    only: kinefsm
       use m_naefs,      only: naefs
@@ -46,8 +47,8 @@
       use m_energies
       use m_steps, only: istp
       use m_ntm
-      use m_spin,         only: spin
-      use spinorbit,      only: spinorb
+      use m_spin,          only: spin
+      use spinorbit,       only: spinorb
       use m_dipol
       use m_forces,          only: fa
       use alloc, only: re_alloc, de_alloc
@@ -91,6 +92,8 @@
       real(dp), pointer     :: H_tmp(:,:) => null()
       real(dp), pointer     :: S_tmp(:) => null()  
       real(dp), pointer     :: H_ldau(:,:)
+      complex(dp), pointer  :: H_so_off(:,:)
+
 #ifdef FINAL_CHECK_HS
       real(dp) :: diff_H, diff_S
 #endif
@@ -98,6 +101,7 @@
       real(dp)              :: stresstmp(3,3)
       real(dp), pointer     :: fatmp(:,:)
       real(dp)              :: buffer1
+      integer               :: MPIerror 
 #endif
       character(len=label_length+13) :: fname
       integer :: io_istep, io_ia1
@@ -106,6 +110,8 @@
       type(dict) :: dic_save
 #endif
 #endif
+
+
 !------------------------------------------------------------------------- BEGIN
 
 !     Initialize Hamiltonian ........................................
@@ -158,28 +164,42 @@
 ! ..................
 
 ! Non-local-pseudop: energy, forces, stress and matrix elements .......
-      call nlefsm( scell, na_u, na_s, isa, xa, indxua, 
-     &             maxnh, maxnh, lasto, lastkb, iphorb, iphKB, 
-     &             numh, listhptr, listh, numh, listhptr, listh, 
-     &             spin%spinor, Dscf, Enl, fal, stressl, H_tmp,
-     &             matrix_elements_only=.false. ) 
-
-#ifdef MPI
-! Global reduction of energy terms
-      call globalize_sum(Enl,buffer1)
-      Enl = buffer1
-#endif
-
-! If in the future the spin-orbit routine is able to compute
-! forces and stresses, then "last" will be needed. If we are not
-! computing forces and stresses, calling it in the first iteration
-! should be enough
-!
-      if ( spin%SO ) then
+
+      Eso = 0.0_dp
+      
+      if ( .not.spin%SO_offsite ) then
+       call nlefsm( scell, na_u, na_s, isa, xa, indxua, 
+     &              maxnh, maxnh, lasto, lastkb, iphorb, iphKB, 
+     &              numh, listhptr, listh, numh, listhptr, listh, 
+     &              spin%spinor, Dscf, Enl, fal, stressl, H_tmp,
+     &              matrix_elements_only=.false. ) 
+
+#ifdef MPI
+! Global reduction of energy terms
+       call globalize_sum(Enl,buffer1)
+       Enl = buffer1
+#endif
+      else
+        H_so_off => val(H_so_off_2D)
+       call nlefsm_SO_off(scell, na_u, na_s, isa, xa, indxua,
+     &                maxnh, maxnh, lasto, lastkb, iphorb, iphKB,
+     &                numh, listhptr, listh, numh, listhptr, listh,
+     &                spin%Grid,
+     &                Enl, Eso, fal,
+     &                stressl, H_tmp, H_so_off,
+     &                matrix_elements_only=.false.)
+#ifdef MPI
+! Global reduction of energy terms
+       call globalize_sum(Enl,buffer1)
+       Enl = buffer1
+       call globalize_sum(Eso,buffer1)
+       Eso = buffer1
+#endif
+      endif
+
+      if ( spin%SO_onsite ) then
          call spinorb(no_u,no_l,iaorb,iphorb,isa,indxuo,
      &        maxnh,numh,listhptr,listh,Dscf,H_tmp(:,3:),Eso)
-      else
-         Eso = 0._dp
       endif
       
 !     Non-SCF part of total energy
@@ -276,7 +296,8 @@
      &        diff_S
       end if
 #endif
-      
+
+      ! TODO I am not sure this works with SO_offsite?
       if (recompute_H_after_scf) then
          if (ionode) then
             write(6,"(a)") ":!: Updating H after scf cycle" //

=== modified file 'Src/initatom.f'
--- Src/initatom.f	2018-04-07 19:24:04 +0000
+++ Src/initatom.f	2018-04-30 20:53:35 +0000
@@ -50,10 +50,10 @@
       use ldau_specs, only: ldau_proj_gen
       use ldau_specs, only: populate_species_info_ldau
       use pseudopotential, only: pseudo_read
-    
+      
       use chemical
 
-      use m_spin, only: SpOrb
+      use m_spin, only: spin
 
       use atm_types, only: species, species_info, nspecies
       
@@ -63,6 +63,7 @@
       integer                      :: is
       logical                      :: user_basis, user_basis_netcdf
       logical :: req_init_setup
+      logical :: lj_projs
 
       type(basis_def_t),   pointer :: basp
       type(species_info),  pointer :: spp
@@ -106,12 +107,13 @@
         call elec_corr_setup()
       else if (user_basis) then
 
-       if ( SpOrb ) then  
+       if ( spin%SO_onsite ) then  
           ! We still need to read the pseudopotential information
-          write(6,'(a)') ' initatom: Spin configuration = spin-orbit'
+          ! because the .ion files do not contain V_so information
+          write(6,'(a)') ' initatom: spin-orbit-onsite with user-basis'
+          write(6,'(a)') ' initatom: Still need to read the psf files.'
           do is = 1 , nsp
              basp => basis_parameters(is)
-             
              basp%label = species_label(is)
              call pseudo_read(basp%label,basp%pseudopotential)
           end do
@@ -132,6 +134,8 @@
         nspecies = nsp              ! For atm_types module
         call setup_atom_tables(nsp)
 
+        lj_projs = (spin%SO_offsite)
+
         allocate(species(nspecies))
         do is = 1,nsp
           call write_basis_specs(6,is)
@@ -151,7 +155,8 @@
      &                    qyuk(0:lmaxd,1:nsemx,is),
      &                    qwid(0:lmaxd,1:nsemx,is),
      &                    split_norm(0:lmaxd,1:nsemx,is), 
-     &                    filtercut(0:lmaxd,1:nsemx,is), basp, spp)
+     &                    filtercut(0:lmaxd,1:nsemx,is), basp, spp,
+     $                    lj_projs)
 !         Generate the projectors for the LDA+U simulations (if requested)
           call ldau_proj_gen(is)
         enddo 
@@ -165,11 +170,12 @@
         call elec_corr_setup()
         ns = nsp               ! Set number of species for main program
 
+        call dump_basis_ascii()
+        call dump_basis_netcdf()
+        call dump_basis_xml()
+        
       endif
 
-      call dump_basis_ascii()
-      call dump_basis_netcdf()
-      call dump_basis_xml()
 
       if (.not. user_basis .and. .not. user_basis_netcdf) then
         call deallocate_spec_arrays()

=== modified file 'Src/m_cite.F90'
--- Src/m_cite.F90	2016-12-29 10:36:10 +0000
+++ Src/m_cite.F90	2018-04-30 20:53:35 +0000
@@ -74,7 +74,7 @@
   ! Increment this after having added a new
   ! citation!
   ! OTHERWISE YOU WILL EXPERIENCE A SEG-FAULT.
-  integer, parameter :: N_citations = 7
+  integer, parameter :: N_citations = 9
 
   private
 
@@ -242,7 +242,29 @@
        cit%issue = "30"
        cit%cite_key = "Lin2014"
        ID = 7
+
+    case ( "10.1088/0953-8984/24/8/086005" )
+       ! Off-Site Spin-Orbit 
+       cit%comment = "Off-Site Spin-Orbit Implementation"
+       cit%doi = "10.1088/0953-8984/24/8/086005"
+       cit%journal = "Journal of Physics: Condensed Matter"
+       cit%year = 2012
+       cit%volume = "24"
+       cit%issue = "8"
+       cit%cite_key = "Cuadrado2012"
+       ID = 8
        
+    case ( "10.1088/0953-8984/18/34/012")
+       ! On-Site Spin-Orbit 
+       cit%comment = "On-Site Spin-Orbit Implementation"
+       cit%doi = "10.1088/0953-8984/18/34/012"
+       cit%journal = "Journal of Physics: Condensed Matter"
+       cit%year = 2006
+       cit%volume = "18"
+       cit%issue = "0"
+       cit%cite_key = "FernandezSeivane2006"
+       ID = 9
+
     end select
 
     ! probably we should error out...

=== modified file 'Src/m_pulay.F90'
--- Src/m_pulay.F90	2016-06-04 20:06:11 +0000
+++ Src/m_pulay.F90	2018-04-30 20:53:35 +0000
@@ -8,7 +8,7 @@
 module m_pulay
 
   use precision, only: dp
-  use m_spin,    only: h_spin_dim, SpOrb
+  use m_spin,    only: spin 
   implicit none
 
   private
@@ -48,7 +48,7 @@
     if (maxsav .le. 0) then
        ! No need for auxiliary arrays
     else
-       nauxpul = sum(numh(1:no_l)) * h_spin_dim * maxsav
+       nauxpul = sum(numh(1:no_l)) * spin%H * maxsav
        !
        call re_alloc(auxpul,1,nauxpul,1,2,name="auxpul",        &
             routine="pulay")
@@ -370,8 +370,8 @@
        ! B(i,i) = dot_product(Res(i)*Res(i))
        b(i,i) = 0.0_dp
        ssum=0.0_dp
-! CC RC Added for on-site SO
-       if ( .not. SpOrb ) then
+
+       if ( .not. spin%SO ) then
         do is=1,nspin
          do ii=1,no_l
           do jj=1,numd(ii)
@@ -380,7 +380,7 @@
           enddo
          enddo
         enddo
-       elseif ( SpOrb ) then
+       elseif ( spin%SO ) then
         do ii=1,no_l
          do jj=1,numd(ii)
           ind = listdptr(ii) + jj
@@ -393,7 +393,6 @@
          enddo
         enddo
        endif
-! CC RC
 
        b(i,i)=ssum
        !
@@ -414,7 +413,7 @@
 
           b(i,j)=0.0_dp
           ssum=0.0_dp
-          if ( .not. SpOrb ) then
+          if ( .not. spin%SO ) then
            do is=1,nspin
             do ii=1,no_l
              do jj=1,numd(ii)
@@ -423,7 +422,7 @@
              enddo
             enddo
            enddo
-          elseif ( SpOrb ) then
+          elseif ( spin%SO ) then
            do ii=1,no_l
             do jj=1,numd(ii)
              ind = listdptr(ii) + jj

=== modified file 'Src/m_spin.F90'
--- Src/m_spin.F90	2018-02-27 21:05:45 +0000
+++ Src/m_spin.F90	2018-04-30 20:53:35 +0000
@@ -7,13 +7,16 @@
 ! ---
 module t_spin
 
+  use precision, only: dp
+  
   implicit none
-
+  private
+  
   !> Type containing a simulations spin configuration
   !!
   !! Thus this type contains _all_ relevant information
   !! regarding the spin-configuration.
-  type tSpin
+  type, public :: tSpin
 
      !> Dimensionality of the Hamiltonian
      integer :: H = 1
@@ -37,8 +40,9 @@
      !> Spin-orbit coupling
      logical :: SO = .false.
 
-     !> If .true. the spin-orbit is using the off-site implementation (else the on-site)
-     logical :: SO_off = .false.
+     !> Flavors of off-site implementation
+     logical :: SO_offsite = .false.
+     logical :: SO_onsite  = .false.
 
      ! Perhaps one could argue that one may
      ! associate a symmetry to the spin which
@@ -110,6 +114,10 @@
     use fdf, only : fdf_get, leqi, fdf_deprecated
     use alloc, only: re_alloc
 
+    use m_cite, only: add_citation
+    use files, only: slabel
+    use parallel, only: IONode
+
     character(len=32) :: opt
 
     ! Create pointer assignments...
@@ -127,14 +135,14 @@
     ! Time reversal symmetry
     TrSym  = .true.
 
-    ! All components of the 'spin' variable
-    ! is initially correct...
+    ! Initialize all components of the 'spin' variable
     spin%none = .false.
     spin%Col = .false.
     spin%NCol = .false.
     spin%SO = .false.
-    spin%SO_off = .false.
-    
+    spin%SO_offsite = .false.
+    spin%SO_onsite  = .false.
+
     ! Read in old flags (discouraged)
     spin%Col  = fdf_get('SpinPolarized',.false.)
     spin%NCol = fdf_get('NonCollinearSpin',.false.)
@@ -149,7 +157,7 @@
     if ( spin%SO ) then
        opt = 'spin-orbit'
     else if ( spin%NCol ) then
-       opt = 'non-colinear'
+       opt = 'non-collinear'
     else if ( spin%Col ) then
        opt = 'collinear'
     else
@@ -173,52 +181,41 @@
        
        spin%Col = .true.
        
-    else if ( leqi(opt, 'non-collinear') .or. leqi(opt, 'non-colinear') .or. &
+    else if ( leqi(opt, 'non-collinear') .or.  leqi(opt, 'non-colinear') .or. &
          leqi(opt, 'NC') .or. leqi(opt, 'N-C') ) then
        
        spin%NCol = .true.
        
     else if ( leqi(opt, 'spin-orbit') .or. leqi(opt, 'S-O') .or. &
-         leqi(opt, 'SOC') .or. leqi(opt, 'SO') ) then
-       ! Spin-orbit is the same as using the off-site implementation
-       
-       spin%SO = .true.
-       spin%SO_off = .true.
-
-    else if ( leqi(opt, 'spin-orbit+off') .or. leqi(opt, 'S-O+off') .or. &
-         leqi(opt, 'SOC+off') .or. leqi(opt, 'SO+off') ) then
-       
-       spin%SO = .true.
-       spin%SO_off = .true.
-
-    else if ( leqi(opt, 'spin-orbit+on') .or. leqi(opt, 'S-O+on') .or. &
-         leqi(opt, 'SOC+on') .or. leqi(opt, 'SO+on') ) then
-       
-       spin%SO = .true.
-       spin%SO_off = .false.
+         leqi(opt, 'SOC') .or. leqi(opt, 'SO') .or. &
+         leqi(opt, 'spin-orbit+offsite') .or. leqi(opt, 'S-O+offsite') .or. &
+         leqi(opt, 'SOC+offsite') .or. leqi(opt, 'SO+offsite') ) then
+       ! Spin-orbit defaults to the off-site implementation
+       
+       spin%SO = .true.
+       spin%SO_offsite = .true.
+
+    else if ( leqi(opt, 'spin-orbit+onsite') .or. leqi(opt, 'S-O+onsite') .or. &
+         leqi(opt, 'SOC+onsite') .or. leqi(opt, 'SO+onsite') ) then
+       
+       spin%SO = .true.
+       spin%SO_onsite = .true.
 
     else
        write(*,*) 'Unknown spin flag: ', trim(opt)
        call die('Spin: unknown flag, please assert the correct input.')
     end if
 
-    ! TODO once off-site is fully implemented
-    ! this should be removed to enable the off-site code
-    ! fully!
-    spin%SO_off = .false.
-
-    ! Note that, in what follows,
-    !   spinor_dim = min(h_spin_dim,2)
-    !   e_spin_dim = min(h_spin_dim,4)
-    !   nspin      = min(h_spin_dim,4)  ! Relevant for dhscf, local DM
-    !      should probably be called nspin_grid
-    !
-    ! so everything can be determined if h_spin_dim is known.
-    ! It is tempting to go back to the old 'nspin' overloading,
-    ! making 'nspin' mean again 'h_spin_dim'.
-    ! But this has to be done carefully, as some routines expect
-    ! an argument 'nspin' which really means 'spinor_dim' (like diagon),
-    ! and others (such as dhscf) expect 'nspin' to mean 'nspin_grid'.
+    if ( spin%SO_offsite ) then
+       call add_citation("10.1088/0953-8984/24/8/086005")
+    end if
+    if ( spin%SO_onsite ) then
+       call add_citation("10.1088/0953-8984/18/34/012")
+    end if
+    
+    ! These are useful for fine control beyond the old "nspin". Some routines expect
+    ! an argument 'nspin' which might really mean 'spin%spinor' (like diagon),
+    ! 'spin%Grid' (like dhscf), etc.
 
     if ( spin%SO ) then
        ! Spin-orbit case
@@ -235,9 +232,10 @@
        spin%Col  = .false.
        spin%NCol = .false.
        spin%SO   = .true.
+       ! off/on MUST already be set!
 
        ! should be moved...
-       TRSym      = .false.
+       TRSym = .false.    ! Cannot be changed !
 
     else if ( spin%NCol ) then
        ! Non-collinear case
@@ -250,13 +248,15 @@
        spin%spinor = 2
 
        ! Flags
-       spin%none = .false.
-       spin%Col  = .false.
-       spin%NCol = .true.
-       spin%SO   = .false.
+       spin%none       = .false.
+       spin%Col        = .false.
+       spin%NCol       = .true.
+       spin%SO         = .false.
+       spin%SO_onsite  = .false.
+       spin%SO_offsite = .false.
 
        ! should be moved...
-       TRSym      = .false.
+       TRSym = .false.    ! Cannot be changed !
 
     else if ( spin%Col ) then
        ! Collinear case
@@ -269,13 +269,17 @@
        spin%spinor = 2
 
        ! Flags
-       spin%none = .false.
-       spin%Col  = .true.
-       spin%NCol = .false.
-       spin%SO   = .false.
+       spin%none       = .false.
+       spin%Col        = .true.
+       spin%NCol       = .false.
+       spin%SO         = .false.
+       spin%SO_onsite  = .false.
+       spin%SO_offsite = .false.
 
        ! should be moved...
-       TRSym      = .true.
+       TRSym = .true.
+       ! Allow the user to choose not to have TRS
+       TRSym  = fdf_get('TimeReversalSymmetry',TrSym)
 
     else if ( spin%none ) then
        ! No spin configuration...
@@ -288,18 +292,20 @@
        spin%spinor = 1
 
        ! Flags
-       spin%none = .true.
-       spin%Col  = .false.
-       spin%NCol = .false.
-       spin%SO   = .false.
+       spin%none       = .true.
+       spin%Col        = .false.
+       spin%NCol       = .false.
+       spin%SO         = .false.
+       spin%SO_onsite  = .false.
+       spin%SO_offsite = .false.
 
        ! should be moved...
-       TRSym      = .true.
+       TRSym = .true.
+       ! Allow the user to choose not to have TRS
+       TRSym  = fdf_get('TimeReversalSymmetry',TrSym)
 
     end if
 
-    ! Get true time reversal symmetry
-    TRSym  = fdf_get('TimeReversalSymmetry',TrSym)
 
   contains
 
@@ -331,10 +337,10 @@
     if ( .not. IONode ) return
 
     if ( spin%SO ) then
-       if ( spin%SO_off ) then
-          opt = 'spin-orbit+off'
+       if ( spin%SO_offsite ) then
+          opt = 'spin-orbit+offsite'
        else
-          opt = 'spin-orbit+on'
+          opt = 'spin-orbit+onsite'
        end if
     else if ( spin%NCol ) then
        opt = 'non-collinear'

=== modified file 'Src/mixer.F'
--- Src/mixer.F	2018-04-17 13:06:00 +0000
+++ Src/mixer.F	2018-04-30 20:53:35 +0000
@@ -30,7 +30,7 @@
       use m_mixing_scf, only: scf_mix
       use m_mixing, only: mixing
       
-      use m_spin,           only: h_spin_dim, SpOrb
+      use m_spin,           only: spin 
       use parallel,         only: IONode
 
       implicit none
@@ -60,7 +60,7 @@
 
 
       ! Create residual function to minimize
-      allocate(F(maxnh,h_spin_dim))
+      allocate(F(maxnh,spin%H))
 
 !$OMP parallel default(shared), private(dtmp,i)
 
@@ -103,15 +103,15 @@
       ! Upon exit Xout contains the mixed quantity
       select case ( mix_spin )
       case ( MIX_SPIN_ALL )
-         call mixing( scf_mix, maxnh, h_spin_dim, Xin, F, Xout)
+         call mixing( scf_mix, maxnh, spin%H, Xin, F, Xout)
       case ( MIX_SPIN_SPINOR )
-         call mixing( scf_mix, maxnh, h_spin_dim, Xin, F, Xout,
+         call mixing( scf_mix, maxnh, spin%H, Xin, F, Xout,
      &        nsub=2)
       case ( MIX_SPIN_SUM )
-         call mixing( scf_mix, maxnh, h_spin_dim, Xin, F, Xout,
+         call mixing( scf_mix, maxnh, spin%H, Xin, F, Xout,
      &        nsub=1)
       case ( MIX_SPIN_SUM_DIFF )
-         call mixing( scf_mix, maxnh, h_spin_dim, Xin, F, Xout,
+         call mixing( scf_mix, maxnh, spin%H, Xin, F, Xout,
      &        nsub=2)
       end select
 
@@ -188,15 +188,15 @@
 
       if (muldeb .and. (.not. MixH) ) then 
          if (IONode) write (6,"(/a)") 'Using DM after mixing:'
-         if ( SpOrb .and. orbmoms) then
+         if ( spin%SO .and. orbmoms) then
             call moments( 1, na_u, no_u, maxnh, numh, listhptr,
      .           listh, S, Dscf, isa, lasto, iaorb, iphorb,
      .           indxuo )
          endif
          ! Call this unconditionally
          call mulliken( mullipop, na_u, no_u, maxnh,
-     &           numh, listhptr, listh, S, Dscf, isa,
-     &           lasto, iaorb, iphorb )
+     &          numh, listhptr, listh, S, Dscf, isa,
+     &          lasto, iaorb, iphorb )
       endif
 
       call timer( 'MIXER', 2 )

=== modified file 'Src/moments.F'
--- Src/moments.F	2018-04-17 13:06:00 +0000
+++ Src/moments.F	2018-04-30 20:53:35 +0000
@@ -296,4 +296,3 @@
       deallocate(l_orb)
 
       end subroutine moments
-

=== modified file 'Src/nlefsm.f'
--- Src/nlefsm.f	2016-11-04 14:10:50 +0000
+++ Src/nlefsm.f	2018-04-30 20:53:35 +0000
@@ -7,9 +7,11 @@
 !
       module m_nlefsm
 
+      use precision,     only : dp
+
       implicit none
 
-      public :: nlefsm
+      public :: nlefsm, nlefsm_SO_off
 
       private
 
@@ -53,7 +55,7 @@
 C                            hamiltonian matrix
 C integer listh(maxnh)     : Nonzero hamiltonian-matrix element column 
 C                            indexes for each matrix row
-C integer nspin            : Number of spin components of Dscf and H
+C integer nspin            : Number of spinor components (really min(nspin,2))
 C                            If computing only matrix elements, it
 C                            can be set to 1.
 C logical matrix_elements_only:
@@ -301,8 +303,6 @@
                  kg = kbproj_gindex(ks,koa)
                  call new_MATEL( 'S', kg, ig, xki(1:3,ina),
      &                Ski(ikb,nno), grSki(1:3,ikb,nno) )
-                 !  Maybe: Ski = epskb_sqrt * Ski
-                 !         grSki = epskb_sqrt * grSki
               enddo
 
            enddo ! loop over orbitals
@@ -330,7 +330,7 @@
              jo = listh(ind)
              listed(jo) = .true.
              if (.not. matrix_elements_only) then
-                do ispin = 1,nspin ! Both spins add up...
+                do ispin = 1, nspin ! Both spins add up... (nspin=spin%spinor here)
                    Di(jo) = Di(jo) + Dscf(ind,ispin)
                 enddo
              endif
@@ -381,7 +381,7 @@
           do j = 1,numh(iio)
              ind = listhptr(iio)+j
              jo = listh(ind)
-             do ispin = 1,nspin
+             do ispin = 1,nspin  ! Only diagonal parts
                 H(ind,ispin) = H(ind,ispin) + Vi(jo)
              enddo
              Vi(jo) = 0.0d0     ! See initial zero-out at top
@@ -441,4 +441,613 @@
       
       end subroutine nlefsm
 
+C nlefsm_SO_off calculates the KB elements to the total Hamiltonian 
+C (including the SO contribution)
+C when Off-Site Spin Orbit is included in the calculation 
+      subroutine nlefsm_SO_off( scell, nua, na, isa, xa, indxua,
+     .                      maxnh, maxnd, lasto, lastkb, iphorb, 
+     .                      iphKB, numd, listdptr, listd, numh, 
+     .                      listhptr, listh, nspin, Enl, Enl_offsiteSO, 
+     .                      fa, stress, H0 , H0_off,
+     &                      matrix_elements_only)
+
+
+C *********************************************************************
+C Calculates non-local (NL) pseudopotential contribution to total 
+C energy, atomic forces, stress and hamiltonian matrix elements.
+C Energies in Ry. Lengths in Bohr.
+C Writen by J.Soler and P.Ordejon, June 1997.
+C Modified by R. Cuadrado and J. I. Cerda for the 
+C off-site Spin-Orbit, January 2017.
+C **************************** INPUT **********************************
+C real*8  scell(3,3)       : Supercell vectors SCELL(IXYZ,IVECT)
+C integer nua              : Number of atoms in unit cell
+C integer na               : Number of atoms in supercell
+C integer isa(na)          : Species index of each atom
+C real*8  xa(3,na)         : Atomic positions in cartesian coordinates
+C integer indxua(na)       : Index of equivalent atom in unit cell
+C integer maxnh            : First dimension of H and listh
+C integer maxnd            : Maximum number of elements of the
+C                            density matrix
+C integer lasto(0:na)      : Position of last orbital of each atom
+C integer lastkb(0:na)     : Position of last KB projector of each atom
+C integer iphorb(no)       : Orbital index of each orbital in its atom,
+C                            where no=lasto(na)
+C integer iphKB(nokb)      : Index of each KB projector in its atom,
+C                            where nokb=lastkb(na)
+C integer numd(nuo)        : Number of nonzero elements of each row of the
+C                            density matrix
+C integer listdptr(nuo)    : Pointer to the start of each row (-1) of the
+C                            density matrix
+C integer listd(maxnd)     : Nonzero hamiltonian-matrix element column 
+C                            indexes for each matrix row
+C integer numh(nuo)        : Number of nonzero elements of each row of the
+C                            hamiltonian matrix
+C integer listhptr(nuo)    : Pointer to the start of each row (-1) of the
+C                            hamiltonian matrix
+C integer listh(maxnh)     : Nonzero hamiltonian-matrix element column 
+C                            indexes for each matrix row
+C integer nspin            : Number of spin_grid components (really min(nspin,4))
+C                            If computing only matrix elements, it
+C                            can be set to 1.
+C logical matrix_elements_only:
+C integer Dscf(maxnd,nspin): Density matrix. Not touched if computing
+C                            only matrix elements.
+C ******************* INPUT and OUTPUT *********************************
+C real*8 fa(3,na)          : NL forces (added to input fa)
+C real*8 stress(3,3)       : NL stress (added to input stress)
+C real*8 H0(maxnh,nspin)    : NL Hamiltonian (added to input H)
+C complex*16 H0_off(maxnh,nspin) : NL off-site Hamiltonian (added to input H)
+C **************************** OUTPUT *********************************
+C real*8 Enl               : NL energy
+C *********************************************************************
+C
+C  Modules
+C
+      use parallel,        only : Node, Nodes
+      use parallelsubs,    only : GetNodeOrbs, LocalToGlobalOrb
+      use parallelsubs,    only : GlobalToLocalOrb
+      use atmfuncs,        only : rcut, orb_gindex, kbproj_gindex
+      use atmfuncs,        only : epskb
+      use neighbour,       only : iana=>jan, r2ki=>r2ij, xki=>xij
+      use neighbour,       only : mneighb, reset_neighbour_arrays
+      use alloc,           only : re_alloc, de_alloc
+      use m_new_matel,     only : new_matel
+      use atm_types,       only: species_info, species
+      use sparse_matrices, only: Dscf, xijo
+     
+      integer, intent(in) ::
+     .   maxnh, na, maxnd, nspin, nua
+
+      integer, intent(in)  ::
+     .  indxua(na), iphKB(*), iphorb(*), isa(na),  
+     .  lasto(0:na), lastkb(0:na), listd(maxnd), listh(maxnh),
+     .  numd(*), numh(*), listdptr(*), listhptr(*)
+
+      real(dp), intent(in)    :: scell(3,3),  xa(3,na)
+      real(dp), intent(inout) :: fa(3,nua), stress(3,3)
+      real(dp), intent(inout) :: H0(maxnh)
+      complex(dp), intent(inout) :: H0_off(maxnh,4)
+                                           
+      real(dp), intent(out)   :: Enl, Enl_offsiteSO 
+      logical, intent(in)     :: matrix_elements_only
+
+      real(dp) ::   volcel
+      external ::   timer, volcel
+
+C Internal variables ................................................
+C maxno  = maximum number of basis orbitals overlapping a KB projector
+
+      integer, save ::  maxno = 2000
+  
+      integer
+     .  ia, ikb, ina, ind, ino, indt,
+     .  io, iio, ioa, is, ispin, ix, ig, kg,
+     .  j, jno, jo, jx, ka, ko, koa, ks, kua,
+     .  nkb, nna, nno, no, nuo, nuotot, maxkba
+
+      integer :: l, koa1, koa2, i 
+
+      integer, dimension(:), pointer :: iano, iono
+
+      real(dp)
+     .  Cijk, fik, rki, rmax, rmaxkb, rmaxo, 
+     .  volume, CVj, epsk(2), Vit  
+
+      real(dp)           :: r1(3), r2(3)
+
+      real(dp), dimension(:), pointer :: Di, Vi ! This is Vion
+      real(dp), dimension(:,:), pointer :: Ski, xno
+      real(dp), dimension(:,:,:), pointer :: grSki
+
+      complex(dp)          :: V_sot(2,2), F_so(3,2,2)
+      complex(dp), pointer :: V_so(:,:,:), Ds(:,:,:)
+
+
+      logical ::   within
+      logical, dimension(:), pointer ::  listed, listedall
+
+      complex(dp) :: E_offsiteSO(4)
+
+      type(species_info), pointer        :: spp
+
+      integer :: nd, ndn
+
+C ------------------------------------------------------------
+
+C Start time counte
+      call timer( 'nlefsm', 1 )
+
+C Find unit cell volume
+      volume = volcel( scell ) * nua / na
+
+C Find maximum range
+      rmaxo = 0.0d0
+      rmaxkb = 0.0d0
+      do ia = 1,na
+        is = isa(ia)
+        do ikb = lastkb(ia-1)+1,lastkb(ia)
+          ioa = iphKB(ikb)
+          rmaxkb = max( rmaxkb, rcut(is,ioa) )
+        enddo
+        do io = lasto(ia-1)+1,lasto(ia)
+          ioa = iphorb(io)
+          rmaxo = max( rmaxo, rcut(is,ioa) )
+        enddo
+      enddo
+      rmax = rmaxo + rmaxkb
+
+C Initialize arrays Di and Vi only once
+      no = lasto(na)
+      nuotot = lasto(nua)
+
+      call GetNodeOrbs(nuotot,Node,Nodes,nuo)
+
+C Allocate local memory
+
+      nullify( Vi ) ! This is Vion
+      call re_alloc( Vi, 1, no, 'Vi', 'nlefsm_SO_off' )
+      Vi(1:no) = 0.0_dp
+
+      nullify( listed )
+      call re_alloc( listed, 1, no, 'listed', 'nlefsm_SO_off' )
+      listed(1:no) = .false.
+      nullify( listedall )
+      call re_alloc( listedall, 1, no, 'listedall', 'nlefsm_SO_off' )
+      listedall(1:no) = .false.
+
+      allocate( V_so(2,2,no) )
+
+      if (.not. matrix_elements_only) then
+       nullify( Di ) 
+       call re_alloc( Di, 1, no, 'Di', 'nlefsm_SO_off' )
+       Di(1:no) = 0.0_dp
+       allocate( Ds(2,2,no) )
+      endif
+
+C Make list of all orbitals needed for this node
+
+      do io = 1,nuo
+        call LocalToGlobalOrb(io,Node,Nodes,iio)
+        listedall(iio) = .true.
+        do j = 1,numh(io)
+          jo = listh(listhptr(io)+j)
+          listedall(jo) = .true.
+        enddo
+      enddo
+
+C Find maximum number of KB projectors of one atom = maxkba
+      maxkba = 0
+      do ka = 1,na
+        nkb = lastkb(ka) - lastkb(ka-1)
+        maxkba = max(maxkba,nkb)
+      enddo
+
+C Allocate local arrays that depend on saved parameters
+      nullify( iano )
+      call re_alloc( iano, 1, maxno, 'iano', 'nlefsm_SO_off' )
+      nullify( iono )
+      call re_alloc( iono, 1, maxno, 'iono', 'nlefsm_SO_off' )
+      nullify( xno )
+      call re_alloc( xno, 1, 3, 1, maxno, 'xno',  'nlefsm_SO_off' )
+      nullify( Ski )
+      call re_alloc( Ski, 1, maxkba, 1, maxno, 'Ski','nlefsm_SO_off')
+      nullify( grSki )
+      call re_alloc( grSki, 1, 3, 1, maxkba, 1, maxno, 'grSki',
+     &               'nlefsm_SO_off' )
+
+C     Initialize neighb subroutine
+      call mneighb( scell, rmax, na, xa, 0, 0, nna )
+
+      nd= 0; ndn= 0
+      Enl = 0.0d0; E_offsiteSO(1:4)=dcmplx(0.0d0,0.0d0)
+      Enl_offsiteSO = 0.0d0
+C     Loop on atoms with KB projectors      
+      do ka = 1,na      ! Supercell atoms
+       kua = indxua(ka) ! Equivalent atom in the UC
+       ks = isa(ka)     ! Specie index of atom ka
+       nkb = lastkb(ka) - lastkb(ka-1) ! number of KB projs of atom ka
+
+C      Find neighbour atoms
+       call mneighb( scell, rmax, na, xa, ka, 0, nna )
+
+C      Find neighbour orbitals
+       Ski(:,:) = 0.0_dp
+       nno = 0; ; iano(:)=0; iono(:)=0
+       do ina = 1,nna  ! Neighbour atoms
+        ia = iana(ina) ! Atom index of ina (the neighbour to ka)
+        is = isa(ia)   ! Specie index of atom ia
+        rki = sqrt(r2ki(ina)) ! Square distance
+ 
+        do io = lasto(ia-1)+1,lasto(ia) ! Orbitals of atom ia
+
+C        Only calculate if needed locally
+
+        if (listedall(io)) then
+          ioa = iphorb(io)  ! Orbital index of orbital io in 
+C                             neighbour atom ina
+
+C         Find if orbital is within range
+          within = .false.
+          do ko = lastkb(ka-1)+1,lastkb(ka)
+           koa = iphKB(ko)
+           if ( rki .lt. rcut(is,ioa)+rcut(ks,koa) ) 
+     &            within = .true.
+          enddo
+
+C         Find overlap between neighbour orbitals and KB projectors
+          if (within) then
+C          Check maxno - if too small then increase array sizes
+           if (nno.eq.maxno) then
+            maxno = maxno + 100
+            call re_alloc( iano, 1, maxno, 'iano', 'nlefsm_SO_off',
+     &                    .true. )
+            call re_alloc( iono, 1, maxno, 'iono', 'nlefsm_SO_off',
+     &                    .true. )
+            call re_alloc( xno, 1, 3, 1, maxno,'xno','nlefsm_SO_off',
+     &                    .true. )
+            call re_alloc( Ski, 1, maxkba, 1, maxno, 'Ski',
+     &                    'nlefsm_SO_off', .true. )
+            call re_alloc( grSki, 1, 3, 1, maxkba, 1, maxno,
+     &                    'grSki', 'nlefsm_SO_off', .true. )
+           endif
+           nno = nno + 1  ! Number of neighbour orbitals
+           iono(nno) = io ! io orbital of atom ina (neighbour to ka)
+           iano(nno) = ia ! atom index of the neighbour ina
+           xno(1:3,nno) = xki(1:3,ina)
+           ikb = 0
+           do ko = lastkb(ka-1)+1,lastkb(ka) !Generic positions of kb's
+            ikb = ikb + 1
+            ioa = iphorb(io)
+            koa = iphKB(ko)            ! koa = -ikb 
+
+            if ( koa.ne.-ikb ) then
+             write(6,*) 'koa ERROR: koa,ikb=',koa,ikb
+             stop
+            endif
+            kg = kbproj_gindex(ks,koa)
+            ig = orb_gindex(is,ioa)
+            call new_MATEL( 'S', kg, ig, xki(1:3,ina),
+     &                     Ski(ikb,nno), grSki(1:3,ikb,nno) )
+           enddo
+          endif  ! Within
+
+         endif    
+        enddo    ! neighbour AO
+       enddo     ! neighbour atoms
+
+C----- Loop on neighbour orbitals
+       do ino = 1,nno
+        io = iono(ino)
+        ia = iano(ino)
+
+        call GlobalToLocalOrb(io,Node,Nodes,iio)
+
+        if (iio.gt.0) then
+C        Valid orbital
+         if (ia .le. nua) then
+          if (.not. matrix_elements_only) then
+           !Scatter density matrix row of orbital io
+           Ds(1:2,1:2,1:no) = dcmplx(0.0d0,0.0d0) 
+           do j = 1,numh(iio)
+            ind = listhptr(iio)+j  ! jptr
+            jo = listh(ind)       ! j
+            Di(jo) = 0.0_dp
+            do ispin = 1,min(2,nspin) ! Only diagonal parts 
+             Di(jo) = Di(jo) + Dscf(ind,ispin)
+            enddo
+            Ds(1,1,jo) = dcmplx(Dscf(ind,1), Dscf(ind,5))  ! D(ju,iu)
+            Ds(2,2,jo) = dcmplx(Dscf(ind,2), Dscf(ind,6))  ! D(jd,id)
+            Ds(1,2,jo) = dcmplx(Dscf(ind,3), Dscf(ind,4))  ! D(ju,id)
+            Ds(2,1,jo) = dcmplx(Dscf(ind,7),-Dscf(ind,8))  ! D(jd,iu)
+           enddo
+          endif
+
+C-------- Scatter filter of desired matrix elements
+          do j = 1,numh(iio)
+           jo = listh(listhptr(iio)+j)
+           listed(jo) = .true.
+          enddo
+
+C  Loading V_ion/V_so
+          Vi(1:no) = 0.0_dp
+          V_so(1:2,1:2,1:no) = dcmplx(0.0d0,0.0d0)
+
+C-------- Find matrix elements with other neighbour orbitals
+          do jno = 1,nno
+           jo = iono(jno)
+
+           if (listed(jo)) then
+
+C---------- Loop on KB projectors
+            ko  = lastkb(ka-1)
+            KB_loop: do
+             koa = -iphKB(ko+1)
+             spp => species(ks)
+             l = spp%pj_l(koa)
+
+c----------- Compute Vion
+             if ( l.eq.0 ) then
+              epsk(1) = epskb(ks,koa)
+              Vit = epsk(1) * Ski(koa,ino) * Ski(koa,jno)
+              Vi(jo) = Vi(jo) + Vit
+              if (.not. matrix_elements_only) then
+               Enl = Enl +  Di(jo) * Vit
+               CVj  = epsk(1) * Ski(koa,jno)
+               Cijk = 2.0_dp * Di(jo) * CVj
+               do ix = 1,3
+                fik = Cijk * grSki(ix,koa,ino)
+                fa(ix,ia)  = fa(ix,ia)  - fik
+                fa(ix,kua) = fa(ix,kua) + fik
+                do jx = 1,3
+                stress(jx,ix) = stress(jx,ix) +
+     &                           xno(jx,ino) * fik / volume
+                enddo
+               enddo
+              endif
+              ko = ko + 1 
+
+c----------- Compute Vion from j+/-1/2 and V_so
+             else
+              koa1 = -iphKB(ko+1)
+              koa2 = -iphKB(ko+2*(2*l+1))
+              epsk(1) = epskb(ks,koa1)
+              epsk(2) = epskb(ks,koa2)
+
+              call calc_Vj_offsiteSO( l, epsk, Ski(koa1:koa2,ino), 
+     &                       Ski(koa1:koa2,jno), grSki(:,koa1:koa2,ino),
+     &                       grSki(:,koa1:koa2,jno), Vit, V_sot, F_so )
+              Vi(jo) = Vi(jo) + Vit
+              V_so(1:2,1:2,jo)= V_so(1:2,1:2,jo) + V_sot(1:2,1:2)
+         
+c------------ Forces & SO contribution to E_NL
+              if (.not. matrix_elements_only) then     
+               Enl = Enl +  Di(jo) * Vit
+
+               E_offsiteSO(1) = E_offsiteSO(1) + V_sot(1,1)*Ds(1,1,jo) ! V(iu,ju)*D(ju,iu)
+               E_offsiteSO(2) = E_offsiteSO(2) + V_sot(2,2)*Ds(2,2,jo) ! V(id,jd)*D(jd,id)
+               E_offsiteSO(3) = E_offsiteSO(3) + V_sot(1,2)*Ds(2,1,jo) ! V(iu,jd)*D(jd,iu)
+               E_offsiteSO(4) = E_offsiteSO(4) + V_sot(2,1)*Ds(1,2,jo) ! V(id,ju)*D(ju,id)
+
+               ! These forces include both the 'ion' and 'SO' parts
+               do ix = 1,3
+                fik = 2.0_dp*dreal(Ds(1,1,jo)*F_so(ix,1,1) +
+     &                             Ds(2,2,jo)*F_so(ix,2,2) +
+     &                             Ds(2,1,jo)*F_so(ix,1,2) +
+     &                             Ds(1,2,jo)*F_so(ix,2,1) )
+                fa(ix,ia)  = fa(ix,ia)  - fik
+                fa(ix,kua) = fa(ix,kua) + fik
+                do jx = 1,3
+                 stress(jx,ix) = stress(jx,ix) +
+     &                           xno(jx,ino) * fik / volume
+                enddo
+               enddo
+              endif
+              ko = ko+2*(2*l+1)
+             endif
+             if ( ko.ge.lastkb(ka) ) exit KB_loop
+            enddo KB_loop
+
+           endif  ! listed
+          enddo ! jno orbitals
+
+C-------- Pick up contributions to H and restore Di and Vi
+          do j = 1,numh(iio)
+           ind = listhptr(iio)+j
+           jo = listh(ind)
+           H0(ind) = H0(ind) + Vi(jo)
+           H0_off(ind,1) = H0_off(ind,1) + V_so(1,1,jo)
+           H0_off(ind,2) = H0_off(ind,2) + V_so(2,2,jo)
+           H0_off(ind,3) = H0_off(ind,3) + V_so(1,2,jo)
+           H0_off(ind,4) = H0_off(ind,4) + V_so(2,1,jo)
+
+
+C     Careful with this Vi()
+           Vi(jo) = 0.0d0
+           listed(jo) = .false.
+          enddo
+         endif ! Atom in UC?
+
+        endif   ! iio .gt. 0 
+       enddo     ! ino AOs
+      enddo       ! atoms with KB projectors loop
+
+      if (.not. matrix_elements_only) then
+       Enl_offsiteSO = sum( dreal(E_offsiteSO(1:4)) )
+      endif
+
+C     Deallocate local memory
+
+      call reset_neighbour_arrays( )
+      call de_alloc( grSki, 'grSki', 'nlefsm_SO_off' )
+      call de_alloc( Ski, 'Ski', 'nlefsm_SO_off' )
+      call de_alloc( xno, 'xno', 'nlefsm_SO_off' )
+      call de_alloc( iono, 'iono', 'nlefsm_SO_off' )
+      call de_alloc( iano, 'iano', 'nlefsm_SO_off' )
+
+      call de_alloc( listedall, 'listedall', 'nlefsm_SO_off' )
+      call de_alloc( listed, 'listed', 'nlefsm_SO_off' )
+      call de_alloc( Vi, 'Vi', 'nlefsm_SO_off' )
+      deallocate( V_so )
+
+      if (.not. matrix_elements_only) then
+         call de_alloc( Di, 'Di', 'nlefsm_SO_off' )
+         deallocate( Ds )
+      endif
+
+      call timer( 'nlefsm', 2 )
+
+      end subroutine nlefsm_SO_off
+
+c-----------------------------------------------------------------------
+c
+!> Evaluates:
+!!   <i|V_NL|j>, where V_NL= Sum_{j,mj} |V,j,mj><V,j,mj|
+c
+c-----------------------------------------------------------------------
+      subroutine calc_Vj_offsiteSO( l, epskb, Ski, Skj, grSki, grSkj,
+     &                       V_ion, V_so, F_so )
+
+      implicit none
+
+      integer     , intent(in)  :: l
+      real(dp)    , intent(in)  :: epskb(2)
+      real(dp)    , intent(in)  :: Ski(-l:l,2), Skj(-l:l,2)
+      real(dp)    , intent(in)  :: grSki(3,-l:l,2), grSkj(3,-l:l,2)
+      real(dp)    , intent(out) :: V_ion
+      complex(dp) , intent(out) :: F_so(3,2,2)
+      complex(dp) , intent(out) :: V_so(2,2)
+
+
+      integer    :: J, ij, imj, m, is
+      real(dp)   :: aj, amj, al, a2l1, fac, facm,
+     &              epskpm, V_iont, cp, cm, facpm
+
+      real(dp)   :: cg(2*(2*l+1),2)
+      complex(dp):: u(-l:l,-l:l)
+      complex(dp):: SVi(2), SVj(2), grSVi(3,2)
+
+      external   :: die
+c-----------------------------------------------------------------------
+
+c---- set constants and factors
+      al   = dble(l)
+      a2l1 = dble( 2*l+1 )
+
+c---- load Clebsch-Gordan coefficients; cg(J,+-)
+      J = 0
+      cg(:,:) = 0.0_dp
+      do ij = 1, 2
+       aj = al + (2*ij-3)*0.5d0        ! j(ij=1)=l-1/2; j(ij=2)=l+1/2
+       facpm= (-1.0d0)**(aj-al-0.5d0)  ! +/- sign
+       do imj = 1, nint(2*aj)+1        ! Degeneracy for j
+        amj = -aj + dfloat(imj-1)      ! mj value
+        J = J+1                        ! (j,mj) index
+
+        cp = sqrt( (al+0.5d0+amj)/a2l1 )
+        cm = sqrt( (al+0.5d0-amj)/a2l1 )
+        if ( ij.eq. 1 ) then
+         cg(J,1) =  cm*facpm ! <j-|up>
+         cg(J,2) =  cp       ! <j-|down>
+        else
+         cg(J,1) =  cp*facpm ! <j+|up>
+         cg(J,2) =  cm       ! <j+|down>
+        endif
+       enddo
+      enddo
+
+c---- Ski(M)= <l,M|i> ; Si(m)= <l,m|i> = u(m,-M)*Ski(-M) + u(m,M)*Ski(M)
+      fac = 1.0d0/sqrt(2.0d0)
+      u(:,:) = cmplx(0.0d0,0.0d0)
+      u(0,0)= cmplx(1.0d0,0.0d0)
+      do m =  1, l
+       facm = fac*(-1.0d0)**m
+       u(-m,+M) = cmplx(1.0d0,0.0d0)*fac
+       u(-m,-M) = cmplx(0.0d0,1.0d0)*fac  ! J. Cerda
+       u(+m,+M) = cmplx(1.0d0,0.0d0)*facm
+       u(+m,-M) =-cmplx(0.0d0,1.0d0)*facm ! J. Cerda
+      enddo
+
+c---- Load V_so
+      V_so= cmplx(0.0d0,0.0d0); F_so= cmplx(0.0d0,0.0d0)
+      J = 0
+      do ij = 1, 2
+       aj = al + (2*ij-3)*0.5d0        ! j value
+       do imj = 1, nint(2*aj)+1        ! Degeneracy for j
+        amj = -aj + dfloat(imj-1)      ! mj value
+        J = J+1                        ! (j,mj) index
+
+        SVi(1:2)= cmplx(0.0d0,0.0d0); SVj(1:2)= cmplx(0.0d0,0.0d0)
+        grSVi(1:3,1:2)= cmplx(0.0d0,0.0d0)
+        do is = 1, 2  ! spin loop
+
+c        select correct m
+         if ( is.eq.1 ) then
+          m = nint(amj-0.5d0)    ! up   => m=mj-1/2
+         else
+          m = nint(amj+0.5d0)    ! down => m=mj+1/2
+         endif
+     
+         if ( iabs(m).le.l ) then
+          SVi(is)= Ski(+M,ij)*u(+m,M)
+          SVj(is)= Skj(+M,ij)*u(+m,M)
+          grSVi(1:3,is)= grSki(1:3,+M,ij)*u(+m,M)
+          if ( m.ne.0 ) then
+           SVi(is)= SVi(is) + Ski(-M,ij)*u(+m,-M)
+           SVj(is)= SVj(is) + Skj(-M,ij)*u(+m,-M)
+           grSVi(1:3,is)= grSVi(1:3,is) +
+     &                    grSki(1:3,-M,ij)*u(+m,-M)
+          endif
+          SVi(is) = SVi(is) * cg(J,is)
+          SVj(is) = SVj(is) * cg(J,is)
+ 
+          grSVi(1:3,is) = grSVi(1:3,is) * cg(J,is)
+         endif
+        enddo ! is
+
+c       up-up = <i,+|V,J><V,J|j,+>
+        V_so(1,1)  = V_so(1,1)  + SVi(1) * epskb(ij) * conjg(SVj(1))
+        F_so(:,1,1)= F_so(:,1,1)+ grSVi(:,1) * epskb(ij) * conjg(SVj(1))
+
+c       down-down = <i,-|V,J><V,J|j,->
+        V_so(2,2)  = V_so(2,2)  + SVi(2) * epskb(ij) * conjg(SVj(2))
+        F_so(:,2,2)= F_so(:,2,2)+ grSVi(:,2) * epskb(ij) * conjg(SVj(2))
+
+c       up-down = <i,+|V,J><V,J|j,->
+        V_so(1,2)  = V_so(1,2)  + SVi(1) * epskb(ij) * conjg(SVj(2))
+        F_so(:,1,2)= F_so(:,1,2)+ grSVi(:,1) * epskb(ij) * conjg(SVj(2))
+
+c       down-up= <i,-|V,J><V,J|j,+>
+        V_so(2,1)  = V_so(2,1)  + SVi(2) * epskb(ij) * conjg(SVj(1))
+        F_so(:,2,1)= F_so(:,2,1)+ grSVi(:,2) * epskb(ij) * conjg(SVj(1))
+
+       enddo ! mj
+      enddo ! ij
+
+cc--- debugging
+      if ( cdabs(V_so(1,2)+conjg(V_so(2,1))).gt.1.0d-4 ) then
+       write(6,'(a,2f12.6)') 'V_so(1,2)=',V_so(1,2)
+       write(6,'(a,2f12.6)') 'V_so(2,1)=',V_so(2,1)
+       call die('calc_Vj_LS: ERROR')
+      endif
+
+c---- substract out V_ion
+      epskpm = sqrt( epskb(1)*epskb(2) )
+      epskpm = sign(epskpm,epskb(1))
+ 
+      V_ion = 0.0d0
+      do M = -l, l
+       V_iont = ( l**2     * Ski(M,1)*epskb(1)*Skj(M,1)
+     &          + (l+1)**2 * Ski(M,2)*epskb(2)*Skj(M,2)
+     &          + l*(l+1)  * Ski(M,1)*epskpm  *Skj(M,2)
+     &          + l*(l+1)  * Ski(M,2)*epskpm  *Skj(M,1) )/(a2l1**2)
+       V_ion = V_ion + V_iont
+      enddo
+   
+      V_so(1,1) = V_so(1,1) - cmplx(1.0d0,0.0d0)*V_ion
+      V_so(2,2) = V_so(2,2) - cmplx(1.0d0,0.0d0)*V_ion
+
+      return
+      end subroutine calc_Vj_offsiteSO
+
       end module m_nlefsm

=== modified file 'Src/setup_H0.F'
--- Src/setup_H0.F	2017-06-23 19:49:18 +0000
+++ Src/setup_H0.F	2018-04-30 20:53:35 +0000
@@ -17,9 +17,12 @@
       
       USE siesta_options,  only: g2cut
       use sparse_matrices, only: H_kin_1D, H_vkb_1D
-      use sparse_matrices, only: H_so_2D
+      use sparse_matrices, only: H_so_on_2D, H_so_off_2D
       use sparse_matrices, only: Dscf
 
+      use m_nlefsm,        only: nlefsm_SO_off
+      use m_spin,          only: spin
+
       use sparse_matrices, only: listh, listhptr, numh, maxnh
       use siesta_geom
       use atmfuncs, only: uion
@@ -40,6 +43,7 @@
       use alloc, only: re_alloc, de_alloc
       use class_dSpData1D, only: val
       use class_dSpData2D, only: val
+      use class_zSpData2D, only: val
 
 #ifdef MPI
       use m_mpi_utils, only: globalize_sum
@@ -52,7 +56,16 @@
       real(dp) :: dummy_E
       integer  :: ia, is
 
-      real(dp), pointer :: H_val(:), H_so(:,:)
+      real(dp)    :: dummy_Eso
+      integer     :: ispin, i, j 
+      complex(dp) :: Dc
+#ifdef MPI
+      real(dp) :: buffer1
+#endif
+
+      real(dp), pointer :: H_val(:), H_so_on(:,:)
+      complex(dp), pointer :: H_so_off(:,:)
+
 
 #ifdef DEBUG
       call write_debug( '    PRE setup_H0' )
@@ -113,15 +126,59 @@
 !$OMP parallel workshare default(shared)
       H_val(:) = 0.0_dp
 !$OMP end parallel workshare
-      call nlefsm(scell, na_u, na_s, isa, xa, indxua, 
-     &            maxnh, maxnh, lasto, lastkb, iphorb, iphKB, 
-     &            numh, listhptr, listh, numh, listhptr, listh, 
-     &            1,
-     &            dummy_dm, dummy_E, dummy_fa, dummy_stress,
-     &            H_val,
-     &            matrix_elements_only=.true.) 
-
-
+
+      Eso = 0.0d0
+
+      if ( .not.spin%SO_offsite ) then
+       call nlefsm(scell, na_u, na_s, isa, xa, indxua, 
+     &             maxnh, maxnh, lasto, lastkb, iphorb, iphKB, 
+     &             numh, listhptr, listh, numh, listhptr, listh, 
+     &             1,
+     &             dummy_dm, dummy_E, dummy_fa, dummy_stress,
+     &             H_val,
+     &             matrix_elements_only=.true.) 
+      else
+        H_so_off => val(H_so_off_2D)
+        H_so_off = dcmplx(0._dp, 0._dp)
+        call nlefsm_SO_off(scell, na_u, na_s, isa, xa, indxua,
+     &                 maxnh, maxnh, lasto, lastkb, iphorb, iphKB,
+     &                 numh, listhptr, listh, numh, listhptr, listh,
+     &                 spin%Grid,
+     &                 dummy_E, dummy_Eso, dummy_fa,
+     &                 dummy_stress, H_val, H_so_off,
+     &                 matrix_elements_only=.true.)
+
+
+!
+!  Dc IS NOT the dense matrix, it is just a complex number 
+! (per each io index) used as an artifact to multiply the 
+! elements of the H_SO times the corresponding elements of 
+! DM in a such way that the result gives Re{Tr[H_SO*DM]}.
+!
+
+        do i = 1, maxnh
+
+!-------- Eso(u,u)
+          Dc = cmplx(Dscf(i,1),Dscf(i,5), dp)
+          Eso = Eso + real( H_so_off(i,1)*Dc, dp)
+!-------- Eso(d,d)
+          Dc = cmplx(Dscf(i,2),Dscf(i,6),dp)
+          Eso = Eso + real( H_so_off(i,2)*Dc, dp)
+!-------- Eso(u,d)
+          Dc = cmplx(Dscf(i,3),Dscf(i,4), dp)
+          Eso = Eso + real( H_so_off(i,4)*Dc, dp)
+!-------- Eso(d,u)
+          Dc = cmplx(Dscf(i,7),-Dscf(i,8), dp)
+          Eso = Eso + real( H_so_off(i,3)*Dc, dp)
+
+        enddo
+
+#ifdef MPI
+! Global reduction of Eso 
+      call globalize_sum(Eso,buffer1)
+      Eso = buffer1
+#endif
+      endif
 ! ..................
 
 ! If in the future the spin-orbit routine is able to compute
@@ -129,17 +186,14 @@
 ! computing forces and stresses, calling it in the first iteration
 ! should be enough
 !
-      if ( spin%SO ) then
-         H_so => val(H_so_2D)
+      if ( spin%SO_onsite ) then
+         H_so_on => val(H_so_on_2D)
 !$OMP parallel workshare default(shared)
-         H_so = 0._dp
+         H_so_on(:,:) = 0._dp
 !$OMP end parallel workshare
          call spinorb(no_u,no_l,iaorb,iphorb,isa,indxuo,
-     &        maxnh,numh,listhptr,listh,Dscf,H_so,Eso)
-      else
-         Eso = 0._dp
+     &        maxnh,numh,listhptr,listh,Dscf,H_so_on,Eso)
       end if
-      
 
 C     This will take care of possible changes to the mesh and atomic-related
 C     mesh structures for geometry changes
@@ -150,7 +204,7 @@
      &                 mscell, G2max, ntm,
      &                 maxnh, numh, listhptr, listh, datm,
      &                 dummy_fa, dummy_stress)
-  
+
       call timer('Setup_H0',2)
 
 #ifdef DEBUG

=== modified file 'Src/setup_hamiltonian.F'
--- Src/setup_hamiltonian.F	2018-03-30 23:16:29 +0000
+++ Src/setup_hamiltonian.F	2018-04-30 20:53:35 +0000
@@ -14,12 +14,14 @@
 
       USE siesta_options
       use sparse_matrices, only: H_kin_1D, H_vkb_1D
-      use sparse_matrices, only: H_ldau_2D, H_so_2D
+      use sparse_matrices, only: H_ldau_2D
+      use sparse_matrices, only: H_so_on_2D, H_so_off_2D
       use sparse_matrices, only: listh, listhptr, numh, maxnh
       use sparse_matrices, only: H, S, Hold
       use sparse_matrices, only: Dscf, Escf, xijo
       use class_dSpData1D,  only: val
       use class_dSpData2D,  only: val
+      use class_zSpData2D,  only: val
 
       use siesta_geom
       use atmfuncs, only: uion
@@ -40,7 +42,9 @@
       use parallel, only: Node
       use m_steps, only: istp
       use m_ntm
-      use m_spin,         only: spin
+
+      use m_spin,          only: spin 
+
       use m_dipol
       use alloc, only: re_alloc, de_alloc
       use m_gamma
@@ -68,7 +72,12 @@
       type(filesOut_t)    :: filesOut  ! blank output file names
       logical             :: use_rhog_in
 
-      real(dp), pointer   :: H_vkb(:), H_kin(:), H_ldau(:,:), H_so(:,:)
+      real(dp), pointer   :: H_vkb(:), H_kin(:), H_ldau(:,:)
+      real(dp), pointer   :: H_so_on(:,:)
+      complex(dp), pointer:: H_so_off(:,:)
+
+      complex(dp):: Dc
+      integer :: ind, i, j
 
 !------------------------------------------------------------------------- BEGIN
 
@@ -84,53 +93,53 @@
       do ispin = 1, spin%H
          do io = 1,maxnh
             Hold(io,ispin) = H(io,ispin)
-         enddo
-      enddo
+         end do
+      end do
 !$OMP end do
 
 !$OMP single
       H_kin => val(H_kin_1D)
       H_vkb => val(H_vkb_1D)
-      if ( spin%SO ) then
-         ! Sadly some compilers (g95), does
-         ! not allow bounds for pointer assignments :(
-         H_so => val(H_so_2D)
+      
+      if ( spin%SO_onsite ) then
+        ! Sadly some compilers (g95), does
+        ! not allow bounds for pointer assignments :(
+        H_so_on => val(H_so_on_2D)
+        
+      else if ( spin%SO_offsite ) then
+        H_so_off => val(H_so_off_2D)
+        
       end if
 !$OMP end single ! keep wait
 
-      ! We do not need to set the non-spinor components
-      ! For non-colinear they are set down below,
-      ! while for spin-orbit they are set to the H_so initial
-      ! spin-orbit.
-      
+      ! Initialize diagonal Hamiltonian
       do ispin = 1, spin%spinor
 !$OMP do
-          do io = 1,maxnh
-            H(io,ispin) = H_kin(io) + H_vkb(io)
-          end do
+        do io = 1,maxnh
+          H(io,ispin) = H_kin(io) + H_vkb(io)
+        end do
 !$OMP end do nowait
       end do
 
-      if ( spin%NCol ) then
-
-!$OMP do collapse(2)
-         do ispin = 3 , spin%H
-            do io = 1, maxnh
-               H(io,ispin) = 0._dp
-            end do
-         end do
-!$OMP end do nowait
-         
-      else if ( spin%SO ) then
-
-!$OMP do collapse(2)
-         do ispin = 3 , spin%H
-            do io = 1, maxnh
-               H(io,ispin) = H_so(io,ispin-2)
-            end do
-         end do
-!$OMP end do nowait
-
+      if ( spin%SO_onsite ) then
+!$OMP do collapse(2)
+        do ispin = 3 , spin%H
+          do io = 1,maxnh
+            H(io,ispin) = H_so_on(io,ispin-2)
+          end do
+        end do
+!$OMP end do nowait
+         
+      else
+
+!$OMP do collapse(2)
+        do ispin = 3 , spin%H
+          do io = 1,maxnh
+            H(io,ispin) = 0._dp
+          end do
+        end do
+!$OMP end do nowait
+         
       end if
          
 ! ..................
@@ -146,6 +155,7 @@
 !$OMP single
       Ekin = 0.0_dp
       Enl = 0.0_dp
+      Eso = 0.0_dp
 !$OMP end single ! keep wait
       
 !$OMP do collapse(2), reduction(+:Ekin,Enl)
@@ -157,21 +167,46 @@
       end do
 !$OMP end do nowait
 
-!$OMP single
-      Eso = 0._dp
-!$OMP end single
-      if ( spin%SO ) then
+!
+!  Dc IS NOT the dense matrix, it is just a complex number 
+! (per each io index) used as an artifact to multiply the 
+! elements of the H_SO times the corresponding elements of 
+! DM in a such way that the result gives Re{Tr[H_SO*DM]}.
+!
+
+      if ( spin%SO_offsite ) then
+        do io = 1, maxnh
+
+!-------- Eso(u,u)
+          Dc = cmplx(Dscf(io,1),Dscf(io,5), dp)
+          Eso = Eso + real( H_so_off(io,1)*Dc, dp)
+!-------- Eso(d,d)
+          Dc = cmplx(Dscf(io,2),Dscf(io,6), dp)
+          Eso = Eso + real( H_so_off(io,2)*Dc, dp)
+!-------- Eso(u,d)
+          Dc = cmplx(Dscf(io,3),Dscf(io,4), dp)
+          Eso = Eso + real( H_so_off(io,4)*Dc, dp)
+!-------- Eso(d,u)
+          Dc = cmplx(Dscf(io,7),-Dscf(io,8), dp)
+          Eso = Eso + real( H_so_off(io,3)*Dc, dp)
+
+        end do
+
+      else if ( spin%SO_onsite ) then
+        
 !$OMP do reduction(+:Eso)
          do io = 1, maxnh
-            Eso = Eso + H_so(io,1)*Dscf(io,7) + H_so(io,2)*Dscf(io,8)
-     .           + H_so(io,5)*Dscf(io,3) + H_so(io,6)*Dscf(io,4)
-     .           - H_so(io,3)*Dscf(io,5) - H_so(io,4)*Dscf(io,6)
+           Eso = Eso + H_so_on(io,1)*Dscf(io,7) +
+     &         H_so_on(io,2)*Dscf(io,8)+ H_so_on(io,5)*Dscf(io,3) +
+     &         H_so_on(io,6)*Dscf(io,4)- H_so_on(io,3)*Dscf(io,5) -
+     &         H_so_on(io,4)*Dscf(io,6)
          end do
 !$OMP end do nowait
-      end if
+
+       end if
 
 !$OMP end parallel
-      
+
 #ifdef MPI
       ! Global reduction of Ekin, Enl
       call globalize_sum(Ekin,buffer1)
@@ -191,7 +226,7 @@
 ! Hubbard term for LDA+U: energy, forces, stress and matrix elements ....
       if( switch_ldau ) then
         if ( spin%NCol ) then
-         call die('LDA+U cannot be used with non-colinear spin.')
+         call die('LDA+U cannot be used with non-collinear spin.')
         end if
         if ( spin%SO ) then
          call die('LDA+U cannot be used with spin-orbit coupling.')
@@ -243,6 +278,25 @@
      .            Exc, Dxc, dipol, stress, fal, stressl,
      .            use_rhog_in)
 
+      if ( spin%SO_offsite ) then
+
+! H(:, [5, 6]) are not updated in dhscf, see vmat for details.
+        
+!------- H(u,u)
+         H(:,1) = H(:,1) + real(H_so_off(:,1), dp)
+         H(:,5) =         dimag(H_so_off(:,1))
+!------- H(d,d)
+         H(:,2) = H(:,2) + real(H_so_off(:,2), dp)
+         H(:,6) =         dimag(H_so_off(:,2))
+!------- H(u,d)
+         H(:,3) = H(:,3) + real(H_so_off(:,3), dp)
+         H(:,4) = H(:,4) +dimag(H_so_off(:,3))
+!------- H(d,u)
+         H(:,7) = H(:,7) + real(H_so_off(:,4), dp)
+         H(:,8) = H(:,8) -dimag(H_so_off(:,4))
+
+      endif
+
       ! This statement will apply to iscf = 1, for example, when
       ! we do not use rhog_in. Rhog here is always the charge used to
       ! build H, that is, rhog_in.

=== modified file 'Src/siesta_options.F90'
--- Src/siesta_options.F90	2018-04-17 13:07:32 +0000
+++ Src/siesta_options.F90	2018-04-30 20:53:35 +0000
@@ -27,7 +27,6 @@
   ! -- pre 4.0 coordinate output logic -- to be implemented
   logical :: compat_pre_v4_dynamics      ! General switch
 
-  
   logical :: mix_scf_first ! Mix first SCF step?
   logical :: mix_charge    ! New: mix fourier components of rho
   logical :: mixH          ! Mix H instead of DM
@@ -234,7 +233,7 @@
   real(dp) :: total_spin    ! Total spin used in spin-polarized calculations
   real(dp) :: tt            ! Target temperature. Read in redata. Used in dynamics rout.
   real(dp) :: wmix          ! Mixing weight for DM in SCF iteration
-  real(dp) :: wmixkick       ! Mixing weight for DM in special 'kick' SCF steps
+  real(dp) :: wmixkick      ! Mixing weight for DM in special 'kick' SCF steps
 
   character(len=164) :: sname   ! System name, used to initialise read
 
@@ -250,5 +249,5 @@
   ! LUA-handle
   type(luaState) :: LUA
 #endif
-  
+
 END MODULE siesta_options

=== modified file 'Src/sparse_matrices.F'
--- Src/sparse_matrices.F	2016-08-18 10:36:36 +0000
+++ Src/sparse_matrices.F	2018-04-30 20:53:35 +0000
@@ -9,6 +9,7 @@
       use precision
       use class_dSpData1D
       use class_dSpData2D
+      use class_zSpData2D
       use class_Sparsity
       use class_OrbitalDistribution
       use class_Fstack_Pair_Geometry_dSpData2D
@@ -31,12 +32,16 @@
 
       real(dp), public, pointer :: xijo(:,:)=>null()
 
+      complex(dp), public, pointer :: H0_offsiteSO(:,:) => null()
 
       ! Pieces of H that do not depend on the SCF density matrix
       ! Formerly there was a single array H0 for this
       type(dSpData1D), public, save :: H_vkb_1D, H_kin_1D
       ! LDA+U and spin-orbit coupling Hamiltonian
-      type(dSpData2D), public, save :: H_ldau_2D, H_so_2D
+      type(dSpData2D), public, save :: H_ldau_2D, H_so_on_2D
+      ! Spin-orbit off-site
+      type(zSpData2D), public, save :: H_so_off_2D
+
 
       !  New interface data
       type(Sparsity), public, save :: sparse_pattern
@@ -54,7 +59,7 @@
       CONTAINS
 
       subroutine resetSparseMatrices( )
-      use alloc, only : de_alloc
+      use alloc,  only : de_alloc
 
       implicit none
 
@@ -65,7 +70,8 @@
       call delete( H_kin_1D )
       call delete( H_vkb_1D )
       call delete( H_ldau_2D )
-      call delete( H_so_2D )
+      call delete( H_so_on_2D )
+      call delete( H_so_off_2D )
 
       call delete( DM_2D ) ; nullify(Dscf)
       call delete( EDM_2D ) ; nullify(Escf)

=== modified file 'Src/spinorbit.f'
--- Src/spinorbit.f	2017-10-26 10:42:57 +0000
+++ Src/spinorbit.f	2018-04-30 20:53:35 +0000
@@ -38,7 +38,6 @@
       use precision, only: dp
       use atmfuncs
       use atm_types
-      use fdf, only: fdf_get
       use parallel, only: Node, Nodes
       use parallelsubs, only: LocalToGlobalOrb
 
@@ -53,7 +52,6 @@
 !     indexing technology)
 
       logical, save  :: vso_setup = .false.
-      real(dp), save :: so_strength = 1.0_dp
 
       integer, pointer, save   :: nr(:)
       real(dp), pointer, save  :: vso(:,:,:)
@@ -73,7 +71,6 @@
       use atmparams,       only: lmaxd
       use parallel,        only: Node
       use m_mpi_utils,     only: broadcast
-      use sys,             only: die
 
       integer  :: is, mx_nrval, li, ir, iup
       real(dp) :: a, b, rpb, ea
@@ -132,7 +129,8 @@
          enddo
          write(6,"(/)")
          if (.not. there_are_so_potentials) then
-            call die("No spin-orbit components for any species!!")
+            write(6,"(a)") "*** WARNING: No spin-orbit components " //
+     $                     "for any species... "
          endif
       endif
 
@@ -140,8 +138,6 @@
       call broadcast(r)
       call broadcast(drdi)
       
-      so_strength = fdf_get('SpinOrbitStrength',1.0_dp)
-
       end subroutine init_vso
 
 !------------------------------------------------
@@ -174,6 +170,7 @@
 C *********************************************************************
 C
       use m_mpi_utils, only: globalize_sum
+
       implicit none
 
 C Arguments
@@ -238,7 +235,7 @@
             
             call int_so_rad(is, li, joa, ioa, int_rad)
             call int_so_ang(li, mj, mi, int_ang(:))
-            Hso_ji(:)=so_strength*int_rad*int_ang(:)
+            Hso_ji(:) = int_rad * int_ang(:)
 
             H(ind,3) = H(ind,3) + Hso_ji(2)
             H(ind,4) = H(ind,4) + Hso_ji(3)

=== modified file 'Src/state_analysis.F'
--- Src/state_analysis.F	2018-04-19 08:46:09 +0000
+++ Src/state_analysis.F	2018-04-30 20:53:35 +0000
@@ -22,7 +22,7 @@
      &                          CartesianForce_to_ZmatForce
       use atomlist,      only : iaorb, iphorb, amass, no_u, lasto
       use atomlist,      only : indxuo
-      use m_spin,        only : nspin, SpOrb
+      use m_spin,        only : spin
       use m_fixed,       only : fixed
       use sparse_matrices
       use siesta_geom
@@ -167,7 +167,7 @@
       endif
 
 !     Population and moment analysis 
-      if ( SpOrb .and. orbmoms) then
+      if ( spin%SO .and. orbmoms) then
          call moments( 1, na_u, no_u, maxnh, numh, listhptr,
      .           listh, S, Dscf, isa, lasto, iaorb, iphorb,
      .           indxuo )

=== modified file 'Src/state_init.F'
--- Src/state_init.F	2018-04-15 14:49:37 +0000
+++ Src/state_init.F	2018-04-30 20:53:35 +0000
@@ -27,7 +27,8 @@
       use sparse_matrices, only: S   , S_1D
 
       use sparse_matrices, only: H_kin_1D, H_vkb_1D
-      use sparse_matrices, only: H_ldau_2D, H_so_2D
+      use sparse_matrices, only: H_ldau_2D
+      use sparse_matrices, only: H_so_on_2D, H_so_off_2D
 
       use sparse_matrices, only: sparse_pattern
       use sparse_matrices, only: block_dist, single_dist
@@ -98,6 +99,7 @@
       use class_Sparsity
       use class_dSpData1D
       use class_dSpData2D
+      use class_zSpData2D
       use class_dData2D
 #ifdef TEST_IO
       use m_test_io
@@ -138,7 +140,6 @@
       type(dData2D) :: tmp_2D
       real(dp) :: dummy_qspin(8)
 
-
 !------------------------------------------------------------------- BEGIN
       call timer( 'IterGeom', 1 )
 #ifdef DEBUG
@@ -339,6 +340,7 @@
       ! be higher than 1, hence we need to create "fake"
       ! containers and let the new<class> delete the old
       ! sparsity pattern
+
       nullify(numh,listhptr,listh)
       allocate(numh(no_l),listhptr(no_l))
       ! We do not need to allocate listh
@@ -562,11 +564,15 @@
          end if
       end if
       
-      if ( spin%SO ) then
-         write(oname,"(a,i0)") "H_so at geom step ", istep
-         call newdSpData2D(sparse_pattern,spin%H - 2,
-     &        block_dist,H_so_2D,name=oname)
-      end if
+      if ( spin%SO_onsite ) then
+        write(oname,"(a,i0)") "H_so (onsite) at geom step ", istep
+        call newdSpData2D(sparse_pattern,spin%H - 2,
+     &      block_dist,H_so_on_2D,name=oname)
+      else if ( spin%SO_offsite ) then
+        write(oname,"(a,i0)") "H_so (offsite) at geom step ", istep
+        call newzSpData2D(sparse_pattern,4,
+     &      block_dist,H_so_off_2D,name=oname)
+      endif
 
       write(oname,"(a,i0)") "S at geom step ", istep
       call newdSpData1D(sparse_pattern,block_dist,S_1D,name=oname)
@@ -585,6 +591,7 @@
 !     Initialize density matrix
       ! The resizing of Dscf is done inside new_dm
       call new_DM(auxchanged, DM_history, DM_2D, EDM_2D)
+
       Dscf => val(DM_2D)
       Escf => val(EDM_2D)
       if (spin%H > 1) call print_spin(dummy_qspin)

=== added directory 'Tests/FePt_soc'
=== added file 'Tests/FePt_soc/FePt_soc.fdf'
--- Tests/FePt_soc/FePt_soc.fdf	1970-01-01 00:00:00 +0000
+++ Tests/FePt_soc/FePt_soc.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,65 @@
+#
+#  -- Caution: quality parameters set artificially low  !
+#
+SystemName              FePt distorted bulk structure -- soc test 
+SystemLabel             FePt_soc
+NumberOfAtoms           2
+NumberOfSpecies         2
+%block Chemical_Species_label
+        1       26      Fe_fept_SOC
+        2       78      Pt_fept_SOC
+%endblock Chemical_Species_label
+
+Spin SO
+
+%block DM.InitSpin
+ 1  +2.   90.  90.
+ 2  +1.   90.  90.
+%endblock DM.InitSpin
+
+LatticeConstant  1.0 Ang
+%block LatticeVectors
+    2.793068700    0.000000000    0.000000000
+    0.000000000    2.70           0.000000000
+    0.000000000    0.000000000    3.792000000
+%endblock LatticeVectors
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+%block AtomicCoordinatesAndAtomicSpecies
+    1.396535500    1.45           0.000000000  1   
+    0.000000000    0.000000000    1.896000000  2   
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+%Block PAO.Basis
+ Fe_fept_SOC   2
+   n=4   0   2   P
+    0.0   0.0
+   n=3   2   2
+    0.0   0.0
+ Pt_fept_SOC   2
+   n=6   0   2   P
+    0.0   0.0
+   n=5   2   2
+    0.0   0.0
+%EndBlock PAO.Basis
+
+%block kgrid_Monkhorst_Pack
+   3   0   0    0.0
+   0   5   0    0.0
+   0   0   5    0.0
+%endblock kgrid_Monkhorst_Pack
+
+xc.functional           GGA
+xc.authors              PBE
+MeshCutoff              400. Ry
+SolutionMethod          diagon
+
+DM.Tolerance            1.0E-4
+MaxSCFIterations        6
+scf-must-converge F
+DM.MixingWeight         0.01
+DM.NumberPulay          4
+scf-mix-spin all
+
+MullikenInSCF  T
+WriteForces            T

=== added file 'Tests/FePt_soc/FePt_soc.pseudos'
--- Tests/FePt_soc/FePt_soc.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/FePt_soc/FePt_soc.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+Fe_fept_SOC
+Pt_fept_SOC

=== added file 'Tests/FePt_soc/makefile'
--- Tests/FePt_soc/makefile	1970-01-01 00:00:00 +0000
+++ Tests/FePt_soc/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=FePt_soc
+include ../test.mk

=== modified file 'Tests/Makefile'
--- Tests/Makefile	2018-04-25 11:39:40 +0000
+++ Tests/Makefile	2018-04-30 20:53:35 +0000
@@ -28,7 +28,7 @@
 #          for a serial run.
 #
 #          It is also possible to have separate working directories,
-#          by using the a "label". For example:
+#          by using a "label". For example:
 #
 #          make label=finer 
 #     
@@ -75,12 +75,8 @@
 #  sih-pexsi-spin
 #  TranSiesta-TBTrans
 
-# These tests are extremely time consuming and
-# should only be runned sometimes
-# Currently they may be runned individually.
-#	SOC_FePt_xx SOC_FePt_xz SOC_FePt_zy SOC_FePt_zz
-tests_soc  = SOC_Pt2_xx SOC_Pt2_xz SOC_Pt2_zy SOC_Pt2_zz
-tests_soc += SOC_FePt_xx SOC_FePt_xz SOC_FePt_zy SOC_FePt_zz
+# SOC tests
+tests_soc = ge_soc_bands Pt_dimer_soc FePt_soc
 
 # Tests only applicable for LUA
 tests_lua  = lua_si111 lua_h2o

=== added directory 'Tests/More_SOC_Examples'
=== added file 'Tests/More_SOC_Examples/README'
--- Tests/More_SOC_Examples/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,10 @@
+These are examples of computation of the Magnetic Anisotropy Energy (MAE)
+and other checks for two systems: a Pt dimer, and a FePt bulk structure.
+
+Note that to run the tests automatically you have to move a specific directory
+one level up in the hierarchy:
+
+  mv offsite_SOC_FePt_xz ..
+  cd ../offsite_SOC_FePt_xz
+  make
+

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_FePt_xx'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xx/FePtxx.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xx/FePtxx.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xx/FePtxx.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,72 @@
+SystemName              FePt-L1_0
+SystemLabel             FePtxx
+NumberOfAtoms           2
+NumberOfSpecies         2
+%block Chemical_Species_label
+        1       26      Fe_fept_SOC
+        2       78      Pt_fept_SOC
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+ Fe_fept_SOC   2
+   n=4   0   2   P
+    0.0   0.0
+   n=3   2   2
+    0.0   0.0
+ Pt_fept_SOC   2
+   n=6   0   2   P
+    0.0   0.0
+   n=5   2   2
+    0.0   0.0
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+%block LatticeVectors
+    3.792000000    0.000000000    0.000000000
+    0.000000000    2.793068700    0.000000000
+    0.000000000    0.000000000    2.793068700
+%endblock LatticeVectors
+LatticeConstant  1.0 Ang
+
+%block AtomicCoordinatesAndAtomicSpecies
+    0.000000000   1.396535500    1.396535500   1   
+    1.896000000   0.000000000    0.000000000   2   
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+Spin SO
+
+%block DM.InitSpin
+ 1  +2.   90.  0.
+ 2  +1.   90.  0.
+%endblock DM.InitSpin
+
+
+%block kgrid_Monkhorst_Pack
+  11   0   0    0.0
+   0  21   0    0.0
+   0   0  21    0.0
+%endblock kgrid_Monkhorst_Pack
+
+xc.functional           GGA
+xc.authors              PBE
+MeshCutoff              1400. Ry
+SolutionMethod          diagon
+ElectronicTemperature   1 meV
+
+DM.Tolerance            1.0E-5
+MaxSCFIterations        600
+DM.MixingWeight         0.005
+DM.NumberPulay          8
+DM.UseSaveDM            T
+DM.MixSCF1  F
+DM.NumberKick           25
+#Diag.DivideAndConquer   F
+MixHamiltonian T
+
+MullikenInSCF  T
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xx/FePtxx.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xx/FePtxx.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xx/FePtxx.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+Fe_fept_SOC
+Pt_fept_SOC

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xx/README'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xx/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xx/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,33 @@
+Ramon Cuadrado del Burgo - March 7, 2016   
+
+  The present test performs fully-relativistic simulation of FePt-L1_0 structure 
+placed with the (001) axis along X and Z, X-alignment and Z-alignment, respectively. 
+By means of the calculation of the total selfconsistent energy it is possible to 
+obtain the magnetic anisotropy energy (MAE) checking in addition that:
+
+ * X-alignment:
+
+   1)  The total SC energy for x-axis spin orientation corresponds to the 
+      lower value: E_x < E_y = E_z => x-axis = Easy axis    
+
+   2)  The total SC energy for y-axis and z-axis spin orientation have the 
+      same value and they are larger than x-axis SC total energy: 
+      E_y = E_z => y/z-axis = Hard axes         
+
+ * Z-alignment:
+
+   1)  The total SC energy for z-axis spin orientation corresponds to the 
+      lower value: E_z < E_x = E_y => z-axis = Easy axis    
+
+   2)  The total SC energy for x-axis and y-axis spin orientation have the 
+      same value and they are larger than z-axis SC total energy: 
+      E_x = E_y => x/y-axis = Hard axes         
+
+ For comparison of (X/Z)-alignment, the total SC energies and MAEs of both 
+configurations have to have the same values when the system is rotated in space:
+
+    a) E_x(X-alignment) = E_z(Z-alignment)
+       E_y(X-alignment) = E_x(Z-alignment)   
+       E_z(X-alignment) = E_y(Z-alignment)   
+ 
+    b) MAE(X-alignment) = E_x - E_(y/z) <=> MAE(Z-alignment) = E_z - E_(x/y)

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xx/makefile'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xx/makefile	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xx/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=FePtxx
+include ../test.mk

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/FePtxx.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/FePtxx.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/FePtxx.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,74 @@
+SystemName              SOC FePt X-alignment x
+SystemLabel             FePtxx
+
+Spin    SO
+%block DM.InitSpin
+ 1  +2.  90.  0.
+ 2  +1.  90.  0.
+%endblock DM.InitSpin
+
+NumberOfAtoms           2
+NumberOfSpecies         2
+%block Chemical_Species_label
+        1       26      Fe
+        2       78      Pt
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+ Fe   2
+   n=4   0   2   P
+    0.0   0.0
+   n=3   2   2
+    0.0   0.0
+ Pt   2
+   n=6     0       2    P
+    0.00000     0.00000
+   n=5     2       2
+    0.00000     0.00000
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+LatticeConstant  1.0 Ang
+%block LatticeVectors
+    3.792000000    0.000000000    0.000000000
+    0.000000000    2.793068700    0.000000000
+    0.000000000    0.000000000    2.793068700
+%endblock LatticeVectors
+
+%block AtomicCoordinatesAndAtomicSpecies
+   0.000000000   1.396535500    1.396535500    1
+   1.896000000   0.000000000    0.000000000    2
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+%block kgrid_Monkhorst_Pack
+  11   0   0    0.0
+   0  21   0    0.0
+   0   0  21    0.0
+%endblock kgrid_Monkhorst_Pack
+
+XC.functional           GGA
+XC.authors              PBE
+
+MeshCutoff              1400. Ry
+
+SolutionMethod          diagon
+
+ElectronicTemperature    1 meV
+
+DM.Tolerance            1.0E-5
+
+MaxSCFIterations        1000
+
+DM.MixingWeight         0.005
+DM.NumberPulay          8
+DM.UseSaveDM            T
+DM.NumberKick           25
+DM.MixSCF1  F
+
+WriteMullikenPop       1
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/FePtxx.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/FePtxx.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/FePtxx.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+Fe_SOC
+Pt_SOC

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/README'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,33 @@
+Ramon Cuadrado del Burgo - March 7, 2016   
+
+  The present test performs fully-relativistic simulation of FePt-L1_0 structure 
+placed with the (001) axis along X and Z, X-alignment and Z-alignment, respectively. 
+By means of the calculation of the total selfconsistent energy it is possible to 
+obtain the magnetic anisotropy energy (MAE) checking in addition that:
+
+ * X-alignment:
+
+   1)  The total SC energy for x-axis spin orientation corresponds to the 
+      lower value: E_x < E_y = E_z => x-axis = Easy axis    
+
+   2)  The total SC energy for y-axis and z-axis spin orientation have the 
+      same value and they are larger than x-axis SC total energy: 
+      E_y = E_z => y/z-axis = Hard axes         
+
+ * Z-alignment:
+
+   1)  The total SC energy for z-axis spin orientation corresponds to the 
+      lower value: E_z < E_x = E_y => z-axis = Easy axis    
+
+   2)  The total SC energy for x-axis and y-axis spin orientation have the 
+      same value and they are larger than z-axis SC total energy: 
+      E_x = E_y => x/y-axis = Hard axes         
+
+ For comparison of (X/Z)-alignment, the total SC energies and MAEs of both 
+configurations have to have the same values when the system is rotated in space:
+
+    a) E_x(X-alignment) = E_z(Z-alignment)
+       E_y(X-alignment) = E_x(Z-alignment)   
+       E_z(X-alignment) = E_y(Z-alignment)   
+ 
+    b) MAE(X-alignment) = E_x - E_(y/z) <=> MAE(Z-alignment) = E_z - E_(x/y)

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/makefile'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/makefile	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xx/offsite_SOC_FePt_xx/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=FePtxx
+include ../test.mk

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_FePt_xz'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xz/FePtxz.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xz/FePtxz.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xz/FePtxz.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,72 @@
+SystemName              FePt-L1_0
+SystemLabel             FePtxz
+NumberOfAtoms           2
+NumberOfSpecies         2
+%block Chemical_Species_label
+        1       26      Fe_fept_SOC
+        2       78      Pt_fept_SOC
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+ Fe_fept_SOC   2
+   n=4   0   2   P
+    0.0   0.0
+   n=3   2   2
+    0.0   0.0
+ Pt_fept_SOC   2
+   n=6   0   2   P
+    0.0   0.0
+   n=5   2   2
+    0.0   0.0
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+%block LatticeVectors
+    3.792000000    0.000000000    0.000000000
+    0.000000000    2.793068700    0.000000000
+    0.000000000    0.000000000    2.793068700
+%endblock LatticeVectors
+LatticeConstant  1.0 Ang
+
+%block AtomicCoordinatesAndAtomicSpecies
+    0.000000000   1.396535500    1.396535500   1   
+    1.896000000   0.000000000    0.000000000   2   
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+Spin SO
+
+#%block DM.InitSpin
+# 1  +2.   0.  0.
+# 2  +1.   0.  0.
+#%endblock DM.InitSpin
+
+
+%block kgrid_Monkhorst_Pack
+  11   0   0    0.0
+   0  21   0    0.0
+   0   0  21    0.0
+%endblock kgrid_Monkhorst_Pack
+
+xc.functional           GGA
+xc.authors              PBE
+MeshCutoff              1400. Ry
+SolutionMethod          diagon
+ElectronicTemperature   1 meV
+
+DM.Tolerance            1.0E-5
+MaxSCFIterations        600
+DM.MixingWeight         0.005
+DM.NumberPulay          8
+DM.UseSaveDM            T
+DM.MixSCF1  F
+DM.NumberKick           25
+#Diag.DivideAndConquer   F
+MixHamiltonian T
+
+MullikenInSCF  T
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xz/FePtxz.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xz/FePtxz.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xz/FePtxz.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+Fe_fept_SOC
+Pt_fept_SOC

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xz/README'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xz/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xz/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,33 @@
+Ramon Cuadrado del Burgo - March 7, 2016   
+
+  The present test performs fully-relativistic simulation of FePt-L1_0 structure 
+placed with the (001) axis along X and Z, X-alignment and Z-alignment, respectively. 
+By means of the calculation of the total selfconsistent energy it is possible to 
+obtain the magnetic anisotropy energy (MAE) checking in addition that:
+
+ * X-alignment:
+
+   1)  The total SC energy for x-axis spin orientation corresponds to the 
+      lower value: E_x < E_y = E_z => x-axis = Easy axis    
+
+   2)  The total SC energy for y-axis and z-axis spin orientation have the 
+      same value and they are larger than x-axis SC total energy: 
+      E_y = E_z => y/z-axis = Hard axes         
+
+ * Z-alignment:
+
+   1)  The total SC energy for z-axis spin orientation corresponds to the 
+      lower value: E_z < E_x = E_y => z-axis = Easy axis    
+
+   2)  The total SC energy for x-axis and y-axis spin orientation have the 
+      same value and they are larger than z-axis SC total energy: 
+      E_x = E_y => x/y-axis = Hard axes         
+
+ For comparison of (X/Z)-alignment, the total SC energies and MAEs of both 
+configurations have to have the same values when the system is rotated in space:
+
+    a) E_x(X-alignment) = E_z(Z-alignment)
+       E_y(X-alignment) = E_x(Z-alignment)   
+       E_z(X-alignment) = E_y(Z-alignment)   
+ 
+    b) MAE(X-alignment) = E_x - E_(y/z) <=> MAE(Z-alignment) = E_z - E_(x/y)

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xz/makefile'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xz/makefile	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xz/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=FePtxz
+include ../test.mk

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/FePtxz.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/FePtxz.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/FePtxz.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,70 @@
+SystemName              SOC FePt X-alignment z
+SystemLabel             FePtxz
+
+Spin    SO
+
+NumberOfAtoms           2
+NumberOfSpecies         2
+%block Chemical_Species_label
+        1       26      Fe
+        2       78      Pt
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+ Fe   2
+   n=4   0   2   P
+    0.0   0.0
+   n=3   2   2
+    0.0   0.0
+ Pt   2
+   n=6     0       2    P
+    0.00000     0.00000
+   n=5     2       2
+    0.00000     0.00000
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+LatticeConstant  1.0 Ang
+%block LatticeVectors
+    3.792000000    0.000000000    0.000000000
+    0.000000000    2.793068700    0.000000000
+    0.000000000    0.000000000    2.793068700
+%endblock LatticeVectors
+
+%block AtomicCoordinatesAndAtomicSpecies
+   0.000000000   1.396535500    1.396535500    1
+   1.896000000   0.000000000    0.000000000    2
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+%block kgrid_Monkhorst_Pack
+  11   0   0    0.5
+   0  21   0    0.0
+   0   0  21    0.0
+%endblock kgrid_Monkhorst_Pack
+
+XC.functional           GGA
+XC.authors              PBE
+
+MeshCutoff              1400. Ry
+
+SolutionMethod          diagon
+
+ElectronicTemperature    1 meV
+
+DM.Tolerance            1.0E-5
+
+MaxSCFIterations        1000
+
+DM.MixingWeight         0.005
+DM.NumberPulay          8
+DM.UseSaveDM            T
+DM.NumberKick           25
+DM.MixSCF1  F
+
+WriteMullikenPop       1
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/FePtxz.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/FePtxz.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/FePtxz.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+Fe_SOC
+Pt_SOC

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/README'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,33 @@
+Ramon Cuadrado del Burgo - March 7, 2016   
+
+  The present test performs fully-relativistic simulation of FePt-L1_0 structure 
+placed with the (001) axis along X and Z, X-alignment and Z-alignment, respectively. 
+By means of the calculation of the total selfconsistent energy it is possible to 
+obtain the magnetic anisotropy energy (MAE) checking in addition that:
+
+ * X-alignment:
+
+   1)  The total SC energy for x-axis spin orientation corresponds to the 
+      lower value: E_x < E_y = E_z => x-axis = Easy axis    
+
+   2)  The total SC energy for y-axis and z-axis spin orientation have the 
+      same value and they are larger than x-axis SC total energy: 
+      E_y = E_z => y/z-axis = Hard axes         
+
+ * Z-alignment:
+
+   1)  The total SC energy for z-axis spin orientation corresponds to the 
+      lower value: E_z < E_x = E_y => z-axis = Easy axis    
+
+   2)  The total SC energy for x-axis and y-axis spin orientation have the 
+      same value and they are larger than z-axis SC total energy: 
+      E_x = E_y => x/y-axis = Hard axes         
+
+ For comparison of (X/Z)-alignment, the total SC energies and MAEs of both 
+configurations have to have the same values when the system is rotated in space:
+
+    a) E_x(X-alignment) = E_z(Z-alignment)
+       E_y(X-alignment) = E_x(Z-alignment)   
+       E_z(X-alignment) = E_y(Z-alignment)   
+ 
+    b) MAE(X-alignment) = E_x - E_(y/z) <=> MAE(Z-alignment) = E_z - E_(x/y)

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/makefile'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/makefile	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_xz/offsite_SOC_FePt_xz/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=FePtxz
+include ../test.mk

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_FePt_zy'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zy/FePtzy.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zy/FePtzy.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zy/FePtzy.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,72 @@
+SystemName              FePt-L1_0
+SystemLabel             FePtzy
+NumberOfAtoms           2
+NumberOfSpecies         2
+%block Chemical_Species_label
+        1       26      Fe_fept_SOC
+        2       78      Pt_fept_SOC
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+ Fe_fept_SOC   2
+   n=4   0   2   P
+    0.0   0.0
+   n=3   2   2
+    0.0   0.0
+ Pt_fept_SOC   2
+   n=6   0   2   P
+    0.0   0.0
+   n=5   2   2
+    0.0   0.0
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+%block LatticeVectors
+    2.793068700    0.000000000    0.000000000
+    0.000000000    2.793068700    0.000000000
+    0.000000000    0.000000000    3.792000000
+%endblock LatticeVectors
+LatticeConstant  1.0 Ang
+
+%block AtomicCoordinatesAndAtomicSpecies
+    1.396535500    1.396535500    0.000000000  1   
+    0.000000000    0.000000000    1.896000000  2   
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+Spin SO
+
+%block DM.InitSpin
+ 1  +2.   90.  90.
+ 2  +1.   90.  90.
+%endblock DM.InitSpin
+
+
+%block kgrid_Monkhorst_Pack
+  21   0   0    0.0
+   0  21   0    0.0
+   0   0  11    0.0
+%endblock kgrid_Monkhorst_Pack
+
+xc.functional           GGA
+xc.authors              PBE
+MeshCutoff              1400. Ry
+SolutionMethod          diagon
+ElectronicTemperature   1 meV
+
+DM.Tolerance            1.0E-5
+MaxSCFIterations        600
+DM.MixingWeight         0.005
+DM.NumberPulay          8
+DM.UseSaveDM            T
+DM.MixSCF1  F
+DM.NumberKick           25
+#Diag.DivideAndConquer   F
+MixHamiltonian T
+
+MullikenInSCF  T
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zy/FePtzy.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zy/FePtzy.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zy/FePtzy.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+Fe_fept_SOC
+Pt_fept_SOC

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zy/README'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zy/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zy/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,33 @@
+Ramon Cuadrado del Burgo - March 7, 2016   
+
+  The present test performs fully-relativistic simulation of FePt-L1_0 structure 
+placed with the (001) axis along X and Z, X-alignment and Z-alignment, respectively. 
+By means of the calculation of the total selfconsistent energy it is possible to 
+obtain the magnetic anisotropy energy (MAE) checking in addition that:
+
+ * X-alignment:
+
+   1)  The total SC energy for x-axis spin orientation corresponds to the 
+      lower value: E_x < E_y = E_z => x-axis = Easy axis    
+
+   2)  The total SC energy for y-axis and z-axis spin orientation have the 
+      same value and they are larger than x-axis SC total energy: 
+      E_y = E_z => y/z-axis = Hard axes         
+
+ * Z-alignment:
+
+   1)  The total SC energy for z-axis spin orientation corresponds to the 
+      lower value: E_z < E_x = E_y => z-axis = Easy axis    
+
+   2)  The total SC energy for x-axis and y-axis spin orientation have the 
+      same value and they are larger than z-axis SC total energy: 
+      E_x = E_y => x/y-axis = Hard axes         
+
+ For comparison of (X/Z)-alignment, the total SC energies and MAEs of both 
+configurations have to have the same values when the system is rotated in space:
+
+    a) E_x(X-alignment) = E_z(Z-alignment)
+       E_y(X-alignment) = E_x(Z-alignment)   
+       E_z(X-alignment) = E_y(Z-alignment)   
+ 
+    b) MAE(X-alignment) = E_x - E_(y/z) <=> MAE(Z-alignment) = E_z - E_(x/y)

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zy/makefile'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zy/makefile	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zy/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=FePtzy
+include ../test.mk

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/FePtzy.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/FePtzy.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/FePtzy.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,68 @@
+SystemName              FePt-L1_0
+SystemLabel             FePtzy
+NumberOfAtoms           2
+NumberOfSpecies         2
+%block Chemical_Species_label
+        1       26      Fe
+        2       78      Pt
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+ Fe   2
+   n=4   0   2   P
+    0.0   0.0
+   n=3   2   2
+    0.0   0.0
+ Pt   2
+   n=6   0   2   P
+    0.0   0.0
+   n=5   2   2
+    0.0   0.0
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+%block LatticeVectors
+    2.793068700    0.000000000    0.000000000
+    0.000000000    2.793068700    0.000000000
+    0.000000000    0.000000000    3.792000000
+%endblock LatticeVectors
+LatticeConstant  1.0 Ang
+
+%block AtomicCoordinatesAndAtomicSpecies
+    1.396535500    1.396535500    0.000000000  1   
+    0.000000000    0.000000000    1.896000000  2   
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+Spin SO
+%block DM.InitSpin
+ 1  +2.  90.  90.
+ 2  +1.  90.  90.
+%endblock DM.InitSpin
+
+%block kgrid_Monkhorst_Pack
+  21   0   0    0.0
+   0  21   0    0.0
+   0   0  11    0.0
+%endblock kgrid_Monkhorst_Pack
+
+xc.functional           GGA
+xc.authors              PBE
+MeshCutoff              1400. Ry
+SolutionMethod          diagon
+ElectronicTemperature   1 meV
+
+DM.Tolerance            1.0E-5
+MaxSCFIterations        600
+DM.MixingWeight         0.005
+DM.NumberPulay          8
+DM.UseSaveDM            T
+DM.MixSCF1  F
+DM.NumberKick           25
+
+WriteMullikenPop       1
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/FePtzy.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/FePtzy.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/FePtzy.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+Fe_SOC
+Pt_SOC

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/README'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,33 @@
+Ramon Cuadrado del Burgo - March 7, 2016   
+
+  The present test performs fully-relativistic simulation of FePt-L1_0 structure 
+placed with the (001) axis along X and Z, X-alignment and Z-alignment, respectively. 
+By means of the calculation of the total selfconsistent energy it is possible to 
+obtain the magnetic anisotropy energy (MAE) checking in addition that:
+
+ * X-alignment:
+
+   1)  The total SC energy for x-axis spin orientation corresponds to the 
+      lower value: E_x < E_y = E_z => x-axis = Easy axis    
+
+   2)  The total SC energy for y-axis and z-axis spin orientation have the 
+      same value and they are larger than x-axis SC total energy: 
+      E_y = E_z => y/z-axis = Hard axes         
+
+ * Z-alignment:
+
+   1)  The total SC energy for z-axis spin orientation corresponds to the 
+      lower value: E_z < E_x = E_y => z-axis = Easy axis    
+
+   2)  The total SC energy for x-axis and y-axis spin orientation have the 
+      same value and they are larger than z-axis SC total energy: 
+      E_x = E_y => x/y-axis = Hard axes         
+
+ For comparison of (X/Z)-alignment, the total SC energies and MAEs of both 
+configurations have to have the same values when the system is rotated in space:
+
+    a) E_x(X-alignment) = E_z(Z-alignment)
+       E_y(X-alignment) = E_x(Z-alignment)   
+       E_z(X-alignment) = E_y(Z-alignment)   
+ 
+    b) MAE(X-alignment) = E_x - E_(y/z) <=> MAE(Z-alignment) = E_z - E_(x/y)

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/makefile'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/makefile	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zy/offsite_SOC_FePt_zy/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=FePtzy
+include ../test.mk

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_FePt_zz'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zz/FePtzz.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zz/FePtzz.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zz/FePtzz.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,72 @@
+SystemName              FePt-L1_0
+SystemLabel             FePtzz
+NumberOfAtoms           2
+NumberOfSpecies         2
+%block Chemical_Species_label
+        1       26      Fe_fept_SOC
+        2       78      Pt_fept_SOC
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+ Fe_fept_SOC   2
+   n=4   0   2   P
+    0.0   0.0
+   n=3   2   2
+    0.0   0.0
+ Pt_fept_SOC   2
+   n=6   0   2   P
+    0.0   0.0
+   n=5   2   2
+    0.0   0.0
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+%block LatticeVectors
+    2.793068700    0.000000000    0.000000000
+    0.000000000    2.793068700    0.000000000
+    0.000000000    0.000000000    3.792000000
+%endblock LatticeVectors
+LatticeConstant  1.0 Ang
+
+%block AtomicCoordinatesAndAtomicSpecies
+    1.396535500    1.396535500    0.000000000  1   
+    0.000000000    0.000000000    1.896000000  2   
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+Spin SO
+
+#%block DM.InitSpin
+# 1  +2.   90.  0.
+# 2  +1.   90.  0.
+#%endblock DM.InitSpin
+
+
+%block kgrid_Monkhorst_Pack
+  21   0   0    0.0
+   0  21   0    0.0
+   0   0  11    0.0
+%endblock kgrid_Monkhorst_Pack
+
+xc.functional           GGA
+xc.authors              PBE
+MeshCutoff              1400. Ry
+SolutionMethod          diagon
+ElectronicTemperature   1 meV
+
+DM.Tolerance            1.0E-5
+MaxSCFIterations        600
+DM.MixingWeight         0.005
+DM.NumberPulay          8
+DM.UseSaveDM            T
+DM.MixSCF1  F
+DM.NumberKick           25
+#Diag.DivideAndConquer   F
+MixHamiltonian T
+
+MullikenInSCF  T
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zz/FePtzz.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zz/FePtzz.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zz/FePtzz.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+Fe_fept_SOC
+Pt_fept_SOC

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zz/README'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zz/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zz/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,33 @@
+Ramon Cuadrado del Burgo - March 7, 2016   
+
+  The present test performs fully-relativistic simulation of FePt-L1_0 structure 
+placed with the (001) axis along X and Z, X-alignment and Z-alignment, respectively. 
+By means of the calculation of the total selfconsistent energy it is possible to 
+obtain the magnetic anisotropy energy (MAE) checking in addition that:
+
+ * X-alignment:
+
+   1)  The total SC energy for x-axis spin orientation corresponds to the 
+      lower value: E_x < E_y = E_z => x-axis = Easy axis    
+
+   2)  The total SC energy for y-axis and z-axis spin orientation have the 
+      same value and they are larger than x-axis SC total energy: 
+      E_y = E_z => y/z-axis = Hard axes         
+
+ * Z-alignment:
+
+   1)  The total SC energy for z-axis spin orientation corresponds to the 
+      lower value: E_z < E_x = E_y => z-axis = Easy axis    
+
+   2)  The total SC energy for x-axis and y-axis spin orientation have the 
+      same value and they are larger than z-axis SC total energy: 
+      E_x = E_y => x/y-axis = Hard axes         
+
+ For comparison of (X/Z)-alignment, the total SC energies and MAEs of both 
+configurations have to have the same values when the system is rotated in space:
+
+    a) E_x(X-alignment) = E_z(Z-alignment)
+       E_y(X-alignment) = E_x(Z-alignment)   
+       E_z(X-alignment) = E_y(Z-alignment)   
+ 
+    b) MAE(X-alignment) = E_x - E_(y/z) <=> MAE(Z-alignment) = E_z - E_(x/y)

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zz/makefile'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zz/makefile	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zz/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=FePtzz
+include ../test.mk

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/FePtzz.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/FePtzz.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/FePtzz.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,64 @@
+SystemName              FePt-L1_0
+SystemLabel             FePtzz
+NumberOfAtoms           2
+NumberOfSpecies         2
+%block Chemical_Species_label
+        1       26      Fe
+        2       78      Pt
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+ Fe   2
+   n=4   0   2   P
+    0.0   0.0
+   n=3   2   2
+    0.0   0.0
+ Pt   2
+   n=6   0   2   P
+    0.0   0.0
+   n=5   2   2
+    0.0   0.0
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+%block LatticeVectors
+    2.793068700    0.000000000    0.000000000
+    0.000000000    2.793068700    0.000000000
+    0.000000000    0.000000000    3.792000000
+%endblock LatticeVectors
+LatticeConstant  1.0 Ang
+
+%block AtomicCoordinatesAndAtomicSpecies
+    1.396535500    1.396535500    0.000000000  1   
+    0.000000000    0.000000000    1.896000000  2   
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+Spin SO
+
+%block kgrid_Monkhorst_Pack
+  21   0   0    0.0
+   0  21   0    0.0
+   0   0  11    0.0
+%endblock kgrid_Monkhorst_Pack
+
+xc.functional           GGA
+xc.authors              PBE
+MeshCutoff              1400. Ry
+SolutionMethod          diagon
+ElectronicTemperature   1 meV
+
+DM.Tolerance            1.0E-5
+MaxSCFIterations        600
+DM.MixingWeight         0.005
+DM.NumberPulay          8
+DM.UseSaveDM            T
+DM.MixSCF1  F
+DM.NumberKick           25
+
+WriteMullikenPop       1
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/FePtzz.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/FePtzz.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/FePtzz.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+Fe_SOC
+Pt_SOC

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/README'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,33 @@
+Ramon Cuadrado del Burgo - March 7, 2016   
+
+  The present test performs fully-relativistic simulation of FePt-L1_0 structure 
+placed with the (001) axis along X and Z, X-alignment and Z-alignment, respectively. 
+By means of the calculation of the total selfconsistent energy it is possible to 
+obtain the magnetic anisotropy energy (MAE) checking in addition that:
+
+ * X-alignment:
+
+   1)  The total SC energy for x-axis spin orientation corresponds to the 
+      lower value: E_x < E_y = E_z => x-axis = Easy axis    
+
+   2)  The total SC energy for y-axis and z-axis spin orientation have the 
+      same value and they are larger than x-axis SC total energy: 
+      E_y = E_z => y/z-axis = Hard axes         
+
+ * Z-alignment:
+
+   1)  The total SC energy for z-axis spin orientation corresponds to the 
+      lower value: E_z < E_x = E_y => z-axis = Easy axis    
+
+   2)  The total SC energy for x-axis and y-axis spin orientation have the 
+      same value and they are larger than z-axis SC total energy: 
+      E_x = E_y => x/y-axis = Hard axes         
+
+ For comparison of (X/Z)-alignment, the total SC energies and MAEs of both 
+configurations have to have the same values when the system is rotated in space:
+
+    a) E_x(X-alignment) = E_z(Z-alignment)
+       E_y(X-alignment) = E_x(Z-alignment)   
+       E_z(X-alignment) = E_y(Z-alignment)   
+ 
+    b) MAE(X-alignment) = E_x - E_(y/z) <=> MAE(Z-alignment) = E_z - E_(x/y)

=== added file 'Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/makefile'
--- Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/makefile	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_FePt_zz/offsite_SOC_FePt_zz/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=FePtzz
+include ../test.mk

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_Pt2_xx'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/Pt2xx.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/Pt2xx.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/Pt2xx.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,64 @@
+SystemName              Pt2 x-alignment x
+SystemLabel             Pt2xx
+
+Spin   SO
+
+%block DM.InitSpin
+ 1  +1.  90.  0.
+ 2  +1.  90.  0.
+%endblock DM.InitSpin
+
+NumberOfAtoms           2
+NumberOfSpecies         1
+%block Chemical_Species_label
+        1       78      Pt_pt2_SOC
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+Pt_pt2_SOC       2          
+ n=6   0    2   P   1  
+   7.158    6.085
+   1.000    1.000
+ n=5   2    2         
+   5.044    3.098
+   1.000    1.000
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+LatticeConstant  20.0 Ang
+%block LatticeVectors
+     1.00      .00      .00
+      .00     1.00      .00
+      .00      .00     1.00
+%endblock LatticeVectors
+
+%block AtomicCoordinatesAndAtomicSpecies
+     -1.19940   0.00000   0.00000    1
+      1.19940   0.00000   0.00000    1
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+XC.functional           GGA
+XC.authors              PBE
+
+MeshCutoff              500. Ry
+
+SolutionMethod          diagon
+
+ElectronicTemperature   1 meV
+DM.Tolerance            1.0E-6
+MaxSCFIterations        1000
+DM.MixingWeight         0.005
+DM.MixSCF1  F
+MixHamiltonian  T
+DM.NumberPulay          6
+DM.NumberKick           25
+DM.UseSaveDM   T
+
+WriteMullikenPop       1
+WriteEigenvalues       T
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/Pt2xx.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/Pt2xx.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/Pt2xx.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,1 @@
+Pt_pt2_SOC

=== added file 'Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/README'
--- Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,34 @@
+Ramon Cuadrado del Burgo - March 7, 2016   
+
+  The present test performs fully-relativistic simulation of Pt2 dimer 
+placed along X and Z axes, X-alignment and Z-alignment, respectively. 
+By means of the calculation of the total selfconsistent energy it is 
+possible to obtain the magnetic anisotropy energy (MAE) checking in 
+addition that:
+
+ * X-alignment:
+
+   1)  The total SC energy for x-axis spin orientation corresponds to the 
+      lower value: E_x < E_y = E_z => x-axis = Easy axis    
+
+   2)  The total SC energy for y-axis and z-axis spin orientation have the 
+      same value and they are larger than x-axis SC total energy: 
+      E_y = E_z => y/z-axis = Hard axes         
+
+ * Z-alignment:
+
+   1)  The total SC energy for z-axis spin orientation corresponds to the 
+      lower value: E_z < E_x = E_y => z-axis = Easy axis    
+
+   2)  The total SC energy for x-axis and y-axis spin orientation have the 
+      same value and they are larger than z-axis SC total energy: 
+      E_x = E_y => x/y-axis = Hard axes         
+
+ For comparison of (X/Z)-alignment, the total SC energies and MAEs of both 
+configurations have to have the same values when the dimer is rotated in space:
+
+    a) E_x(X-alignment) = E_z(Z-alignment)
+       E_y(X-alignment) = E_x(Z-alignment)   
+       E_z(X-alignment) = E_y(Z-alignment)   
+ 
+    b) MAE(X-alignment) = E_x - E_(y/z) <=> MAE(Z-alignment) = E_z - E_(x/y)

=== added file 'Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/makefile'
--- Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/makefile	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_Pt2_xx/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=Pt2xx
+include ../test.mk

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_Pt2_xz'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/Pt2xz.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/Pt2xz.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/Pt2xz.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,59 @@
+SystemName              Pt2 x-alignment z
+SystemLabel             Pt2xz
+
+Spin   SO
+
+NumberOfAtoms           2
+NumberOfSpecies         1
+%block Chemical_Species_label
+        1       78      Pt_pt2_SOC
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+Pt_pt2_SOC          2          
+ n=6   0    2   P   1  
+   7.158    6.085
+   1.000    1.000
+ n=5   2    2         
+   5.044    3.098
+   1.000    1.000
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+LatticeConstant  20.0 Ang
+%block LatticeVectors
+     1.00      .00      .00
+      .00     1.00      .00
+      .00      .00     1.00
+%endblock LatticeVectors
+
+%block AtomicCoordinatesAndAtomicSpecies
+     -1.19940   0.00000   0.00000    1
+      1.19940   0.00000   0.00000    1
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+XC.functional           GGA
+XC.authors              PBE
+
+MeshCutoff              500. Ry
+
+SolutionMethod          diagon
+
+ElectronicTemperature   1 meV
+DM.Tolerance            1.0E-6
+MaxSCFIterations        1000
+DM.MixingWeight         0.005
+DM.MixSCF1  F
+MixHamiltonian  T
+DM.NumberPulay          6
+DM.NumberKick           25
+DM.UseSaveDM   T
+
+WriteMullikenPop       1
+WriteEigenvalues       T
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/Pt2xz.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/Pt2xz.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/Pt2xz.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,1 @@
+Pt_pt2_SOC

=== added file 'Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/README'
--- Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/README	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/README	2018-04-30 20:53:35 +0000
@@ -0,0 +1,34 @@
+Ramon Cuadrado del Burgo - March 7, 2016   
+
+  The present test performs fully-relativistic simulation of Pt2 dimer 
+placed along X and Z axes, X-alignment and Z-alignment, respectively. 
+By means of the calculation of the total selfconsistent energy it is 
+possible to obtain the magnetic anisotropy energy (MAE) checking in 
+addition that:
+
+ * X-alignment:
+
+   1)  The total SC energy for x-axis spin orientation corresponds to the 
+      lower value: E_x < E_y = E_z => x-axis = Easy axis    
+
+   2)  The total SC energy for y-axis and z-axis spin orientation have the 
+      same value and they are larger than x-axis SC total energy: 
+      E_y = E_z => y/z-axis = Hard axes         
+
+ * Z-alignment:
+
+   1)  The total SC energy for z-axis spin orientation corresponds to the 
+      lower value: E_z < E_x = E_y => z-axis = Easy axis    
+
+   2)  The total SC energy for x-axis and y-axis spin orientation have the 
+      same value and they are larger than z-axis SC total energy: 
+      E_x = E_y => x/y-axis = Hard axes         
+
+ For comparison of (X/Z)-alignment, the total SC energies and MAEs of both 
+configurations have to have the same values when the dimer is rotated in space:
+
+    a) E_x(X-alignment) = E_z(Z-alignment)
+       E_y(X-alignment) = E_x(Z-alignment)   
+       E_z(X-alignment) = E_y(Z-alignment)   
+ 
+    b) MAE(X-alignment) = E_x - E_(y/z) <=> MAE(Z-alignment) = E_z - E_(x/y)

=== added file 'Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/makefile'
--- Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/makefile	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_Pt2_xz/makefile	2018-04-30 20:53:35 +0000
@@ -0,0 +1,2 @@
+name=Pt2xz
+include ../test.mk

=== added directory 'Tests/More_SOC_Examples/offsite_SOC_Pt2_zy'
=== added file 'Tests/More_SOC_Examples/offsite_SOC_Pt2_zy/Pt2zy.fdf'
--- Tests/More_SOC_Examples/offsite_SOC_Pt2_zy/Pt2zy.fdf	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_Pt2_zy/Pt2zy.fdf	2018-04-30 20:53:35 +0000
@@ -0,0 +1,64 @@
+SystemName              Pt2 z-alignment y
+SystemLabel             Pt2zy
+
+Spin   SO
+
+%block DM.InitSpin
+ 1  +1.   90.  90.
+ 2  +1.   90.  90.
+%endblock DM.InitSpin
+
+NumberOfAtoms           2
+NumberOfSpecies         1
+%block Chemical_Species_label
+        1       78      Pt_pt2_SOC
+%endblock Chemical_Species_label
+
+PAO.EnergyShift  100 meV
+PAO.SplitNorm      0.15
+PAO.OldStylePolOrbs F
+Restricted.Radial.Grid  F
+%Block PAO.Basis
+Pt_pt2_SOC          2          
+ n=6   0    2   P   1  
+   7.158    6.085
+   1.000    1.000
+ n=5   2    2         
+   5.044    3.098
+   1.000    1.000
+%EndBlock PAO.Basis
+
+AtomicCoordinatesFormat NotScaledCartesianAng
+LatticeConstant  20.0 Ang
+%block LatticeVectors
+     1.00      .00      .00
+      .00     1.00      .00
+      .00      .00     1.00
+%endblock LatticeVectors
+
+%block AtomicCoordinatesAndAtomicSpecies
+      0.00000   0.00000   -1.19940   1
+      0.00000   0.00000    1.19940   1
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+XC.functional           GGA
+XC.authors              PBE
+
+MeshCutoff              500. Ry
+
+SolutionMethod          diagon
+
+ElectronicTemperature   1 meV
+DM.Tolerance            1.0E-6
+MaxSCFIterations        1000
+DM.MixingWeight         0.005
+DM.MixSCF1  F
+MixHamiltonian  T
+DM.NumberPulay          6
+DM.UseSaveDM   T
+DM.NumberKick           25
+
+WriteMullikenPop       1
+WriteEigenvalues       T
+WriteForces            T
+WriteCoorStep          T

=== added file 'Tests/More_SOC_Examples/offsite_SOC_Pt2_zy/Pt2zy.pseudos'
--- Tests/More_SOC_Examples/offsite_SOC_Pt2_zy/Pt2zy.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/More_SOC_Examples/offsite_SOC_Pt2_zy/Pt2zy.pseudos	2018-04-30 20:53:35 +0000
@@ -0,0 +1,1 @@
+Pt_pt2_SOC