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}, 7999 (2006); L
> Fern\'andez--Seivane
> > +and Jaime Ferrer, Phys. Rev. Lett. {\bf 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 {\bf 24}, 086005
> (2012);
>
> same here
>
> > +``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, {\bf 108},
>
> same here
>
> > +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, H$^{\sigma\sigma'}_{\mu\nu}$~(DM$^{\sigma\sigma'}_{\mu\nu}$)
>
> I would just put everything in $$, otherwise do \mathbf H
>
> >  [$\sigma,\sigma'$=$\uparrow,\downarrow$] will be $\neq$0. 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 H$^{\sigma\sigma'}_{\mu\nu}$~(and
>
> Everything in $$
>
> > +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 %% only supported by the off--site SO
> > +    %% formalism *** Why ?
> > +  %
> > +  %%% *** 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
> > +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 <10$^{-5}$eV), \fdf{ElectronicTemperature},
> > -\textbf{k}-point sampling and high values of
> > +selfconsistency~(typically between 10$^{-3}$eV -- 10$^{-4}$eV),
>
> everything in $$
>
> > +\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,
> > +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 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.
> > -
> > +into account is the mixing scheme employed. You are encouraged to use
> > +\fdf{SCF.Mix} \fdf*{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.
> >
> >  \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.
> >
> >  \end{fdflogicalF}
> >
> > +% *** The following items should not be relevant.
> > +
> > +%%For the off--site SO approximation some mandatory flags have to set
> > +%%up to \fdffalse\ in the fdf file:
> > +
> > +%%\begin{fdflogicalT}{PAO!OldStylePolOrbs}
> > +
> > +%%  By default is set up to \fdftrue, however it has to set up to
> \fdffalse.
>
> delete this
>
> > +
> > +%%\end{fdflogicalT}
> > +
> > +%%\begin{fdflogicalT}{DM.MixSCF1}
> > +
> > +%%  By default is set up to \fdftrue, however it has to set up to
> \fdffalse.
> > +
> > +%%\end{fdflogicalT}
> > +
> > +%%\begin{fdflogicalT}{Restricted.Radial.Grid}
> > +
> > +%%  By default is set up to \fdftrue, however it has to set up to
> \fdffalse.
> > +
> > +%%\end{fdflogicalT}
> > +
> > +
> >
> >  \subsection{The self-consistent-field loop}
> >
> >
> > === modified file 'Src/atom.F'
> > --- Src/atom.F        2018-01-02 09:06:31 +0000
> > +++ Src/atom.F        2018-04-20 09:43:43 +0000
> > @@ -1840,7 +1842,7 @@
> >       .           'GHOST: WARNING: Ghost state for L =', l
> >              ighost=1
> >            else
> > -            write(6,'(a,i3)') 'GHOST: No ghost state for L =',l
> > +            !write(6,'(a,i3)') 'GHOST: No ghost state for L =',l
>
> commented out?
>
> >            endif
> >
> >          elseif (dkbcos.eq.0.0d0) then
> >
> > === modified file 'Src/bands.F'
> > --- Src/bands.F       2017-10-10 19:27:53 +0000
> > +++ Src/bands.F       2018-04-20 09:43:43 +0000
> > @@ -366,7 +366,8 @@
> >        use atmfuncs,     only : symfio, cnfigfio, labelfis, nofis
> >        use writewave,    only : wfs_filename
> >
> > -      use m_spin,       only: NoMagn, SPpol, NonCol, SpOrb
> > +!      use m_spin,       only: NoMagn, SPpol, NonCol, SpOrb
>
> remove comment
>
> > +      use m_spin,       only: spin
> >
> >        use m_diag_option, only: ParallelOverK, Serial
> >        use m_diag, only: diag_init
> >
> > === modified file 'Src/compute_energies.F90'
> > --- Src/compute_energies.F90  2017-12-21 15:49:49 +0000
> > +++ Src/compute_energies.F90  2018-04-20 09:43:43 +0000
> > @@ -170,36 +239,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) )      ! more
> clear: -4 here...
>
> Yes, what abuot the sign here?
>
> > +          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) )  ! ... and
> -4 here, as in Ebs
> > +          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) )  ! ...
> and -4 here.
> > +          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
> >
> > === modified file 'Src/dfscf.f'
> > --- Src/dfscf.f       2016-04-05 09:48:40 +0000
> > +++ Src/dfscf.f       2018-04-20 09:43:43 +0000
> > @@ -70,8 +71,12 @@
> >        use parallel,      only: Nodes, Node
> >        use sys,           only: die
> >        use parallelsubs,  only: GlobalToLocalOrb
> > -      use m_spin,        only: NonCol, SpOrb
> > -      use m_spin,        only: h_spin_dim, spinor_dim
> > +!      use m_spin,        only: NonCol, SpOrb
> > +!      use m_spin,        only: h_spin_dim, spinor_dim
> > +
> > +! CC RC  Added  for the offSpOrb
>
> clean comments
>
> > +      use m_spin,        only: spin
> > +! CC RC  Added  for the offSpOrb
> >
> >        implicit none
> >
> > @@ -276,15 +278,13 @@
> >                  jb = ibuff(j)
> >                  D(ib,jb,:) = Dscf(ind,:)
> >  ! JFR           ! We cannot assume symmetry (ib,jb) -- (jb,ib)
> > -                if ( SpOrb ) then
> > -                   D(jb,ib,1) = Dscf(ind,1)
> > -                   D(jb,ib,2) = Dscf(ind,2)
> > -                   D(jb,ib,3) = Dscf(ind,7)
> > -                   D(jb,ib,4) = Dscf(ind,8)
> > -                   D(jb,ib,5) = -Dscf(ind,5)
> > -                   D(jb,ib,6) = -Dscf(ind,6)
> > -                   D(jb,ib,7) = Dscf(ind,3)
> > -                   D(jb,ib,8) = Dscf(ind,4)
> > +                if ( spin%SO ) then
> > +                   D(jb,ib,1) =  Dscf(ind,1)
> > +                   D(jb,ib,2) =  Dscf(ind,2)
> > +                   D(jb,ib,3) =  Dscf(ind,7)
> > +                   D(jb,ib,4) =  Dscf(ind,8)
> > +                   D(jb,ib,7) =  Dscf(ind,3)
> > +                   D(jb,ib,8) =  Dscf(ind,4)
>
> what about signs here?
>
> >                  else
> >                     D(jb,ib,:) = Dscf(ind,:)
> >                  endif
> > @@ -354,13 +354,14 @@
> >          ! done with nspin=4 (grid_nspin)?
> >          V(1:nsp,1:nspin) = Vscf(1:nsp,ip,1:nspin)
> >
> > +!! CC RC  Added  for the offSpOrb
>
> clean comments
>
> >  C     Factor two for nondiagonal elements for non-collinear spin
> > -        if ( NonCol .or. SpOrb ) then
> > -           V(1:nsp,3:4) = 2.0_dp * V(1:nsp,3:4)
> > -           if ( SpOrb ) then
> > -              V(1:nsp,7:8) = V(1:nsp,3:4)
> > -           end if
> > -        end if
> > +        if( spin%SO ) then
> > +         V(1:nsp,7:8) = V(1:nsp,3:4)
> > +        elseif( spin%NCol ) then
> > +         V(1:nsp,3:4) = 2.0_dp * V(1:nsp,3:4)
> > +        endif
> > +!! CC RC  Added  for the offSpOrb
> >
> >  C     Loop on first orbital of mesh point
> >          do ic = 1,nc
> >
> > === modified file 'Src/m_spin.F90'
> > --- Src/m_spin.F90    2018-02-27 21:05:45 +0000
> > +++ Src/m_spin.F90    2018-04-20 09:43:43 +0000
> > @@ -173,40 +178,42 @@
> >
> >         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.
> > +       ! Spin-orbit defaults to the off-site implementation
> > +
> > +       spin%SO = .true.
> > +       spin%SO_offsite = .true.
> > +
> > +    else if ( leqi(opt, 'spin-orbit+offsite') .or. leqi(opt,
> 'S-O+offsite') .or. &
> > +         leqi(opt, 'SOC+offsite') .or. leqi(opt, 'SO+offsite') ) then
> > +
>
> I think these two if-statements should be merged
>
> > +       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.
> > -
> > +    if (spin%SO_offsite) then
> > +       call add_citation("10.1088/0953-8984/24/8/086005")
> > +    endif
> > +    if (spin%SO_onsite) then
> > +       call add_citation("10.1088/0953-8984/18/34/012")
> > +    endif
> > +
> >      ! Note that, in what follows,
> >      !   spinor_dim = min(h_spin_dim,2)
> >      !   e_spin_dim = min(h_spin_dim,4)
> >
> > === modified file 'Src/rhoofd.F90'
> > --- Src/rhoofd.F90    2016-09-13 13:22:20 +0000
> > +++ Src/rhoofd.F90    2018-04-20 09:43:43 +0000
> > @@ -304,7 +329,17 @@
> >             iphi = iphorb(i)
> >
> >  !          Retrieve phi values
> > -           Clocal(:,ic) = dsqrt(2._dp) * phia(iphi,:)
> > +           if ( spin%SO_offsite ) then
> > +              Clocal(:,ic) = phia(iphi,:)
> > +!            if(ip.lt.10) then
> > +!             write(6,*) ' ip/nc/Clocal(1,:) = ', ip, nc, Clocal(1,:)
> > +!            else
> > +!             stop 'Stopping in rhoofd...'
>
> comments?
>
> > +!            endif
> > +
> > +           else
> > +              Clocal(:,ic) = dsqrt(2._dp) * phia(iphi,:)
> > +           end if
> >
> >  !          Loop on second orbital of mesh point
> >             do jc = 1, ic - 1
> > @@ -343,7 +387,18 @@
> >             ilc(ic) = il
> >
> >  !          Retrieve phi values
> > -           Clocal(:,ic) = dsqrt(2._dp) * phi(:,imp)
> > +           if ( spin%SO_offsite ) then
> > +              Clocal(:,ic) = phi(:,imp)
> > +!            if(ip.lt.5) then
> > +!             write(6,*) ' ip/nc/Clocal(1,:) = ', ip, nc, Clocal(1,:)
> > +!            else
> > +!             stop 'Stopping in rhoofd...'
> > +!            endif
>
> comments?
>
> > +
> > +
> > +           else
> > +              Clocal(:,ic) = dsqrt(2._dp) * phi(:,imp)
> > +           end if
> >
> >  !          Loop on second orbital of mesh point
> >             do jc = 1, ic - 1
> >
> > === modified file 'Src/state_analysis.F'
> > --- Src/state_analysis.F      2018-04-19 08:46:09 +0000
> > +++ Src/state_analysis.F      2018-04-20 09:43:43 +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 : SpOrb
>
> change to spin
>
> >        use m_fixed,       only : fixed
> >        use sparse_matrices
> >        use siesta_geom
>
>
> --
> https://code.launchpad.net/~albertog/siesta/merge-OSSO/+merge/343501
> You are the owner of lp:~albertog/siesta/merge-OSSO.
>