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

[Nick Papior]

Transferred H0_offsiteSO to H_so_off_2D as a class-object.

Fixed some possible inconsistencies in the conversion between
complex and real/imaginary part without kind specifications.

Clarified usage in setup_hamiltonian. Now the structure has been
re-assigned as prior to the OSSO implementation.

Transferred H0_offsiteSO to H_so_off_2D as a class-object.

Fixed some possible inconsistencies in the conversion between
complex and real/imaginary part without kind specifications.

Clarified usage in setup_hamiltonian. Now the structure has been
re-assigned as prior to the OSSO implementation.

[Alberto Garcia (albertog)]

Great. Thanks!

review: Approve

=== modified file 'Src/compute_energies.F90'
--- Src/compute_energies.F90	2018-04-23 07:36:55 +0000
+++ Src/compute_energies.F90	2018-04-23 11:10:59 +0000
@@ -48,8 +48,7 @@
                                  rmaxo, no_l
       use m_ntm,           only: ntm
 
-      use m_spin,          only: spin 
-      use sparse_matrices, only: H0_offsiteSO
+      use m_spin,          only: spin
       use parallel,        only: IONode 
 
       use m_dipol,         only: dipol
@@ -300,10 +299,15 @@
       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)
@@ -356,41 +360,37 @@
       Eso = 0._dp
 
       if ( spin%SO_offsite ) then
+         H_so_off => val(H_so_off_2D)
 
-         ! The computation of the trace is different here, as H0_offsiteSO has
+         ! 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),kind=dp)
-            Hc = H0_offsiteSO(io,1)
-          Eso = Eso + real( Hc*Dc )
+          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),kind=dp)
-            Hc = H0_offsiteSO(io,2)
-          Eso = Eso + real( Hc*Dc )
+          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),kind=dp)
-            Hc = H0_offsiteSO(io,4)
-          Eso = Eso + real( Hc*Dc )
+          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),kind=dp)
-            Hc = H0_offsiteSO(io,3)
-          Eso = Eso + real( Hc*Dc )
-
-        enddo
-
-      endif
-
-      if ( spin%SO_onsite ) then
+          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

=== modified file 'Src/final_H_f_stress.F'
--- Src/final_H_f_stress.F	2018-04-20 09:41:55 +0000
+++ Src/final_H_f_stress.F	2018-04-23 11:10:59 +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, nlefsm_offsiteSO
+      use m_nlefsm,     only: nlefsm, nlefsm_SO_off
       use m_overfsm,    only: overfsm
       use m_kinefsm,    only: kinefsm
       use m_naefs,      only: naefs
@@ -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
@@ -176,13 +179,14 @@
        call globalize_sum(Enl,buffer1)
        Enl = buffer1
 #endif
-      else 
-       call nlefsm_offsiteSO(scell, na_u, na_s, isa, xa, indxua,
+      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,
+     &                stressl, H_tmp, H_so_off,
      &                matrix_elements_only=.false.)
 #ifdef MPI
 ! Global reduction of energy terms
@@ -193,7 +197,7 @@
 #endif
       endif
 
-      if ( spin%SO_onsite) then
+      if ( spin%SO_onsite ) then
          call spinorb(no_u,no_l,iaorb,iphorb,isa,indxuo,
      &        maxnh,numh,listhptr,listh,Dscf,H_tmp(:,3:),Eso)
       endif
@@ -292,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/nlefsm.f'
--- Src/nlefsm.f	2018-04-23 09:26:04 +0000
+++ Src/nlefsm.f	2018-04-23 11:10:59 +0000
@@ -8,11 +8,10 @@
       module m_nlefsm
 
       use precision,     only : dp
-      use sparse_matrices, only: H0_offsiteSO
 
       implicit none
 
-      public :: nlefsm, nlefsm_offsiteSO, calc_Vj_offsiteSO
+      public :: nlefsm, nlefsm_SO_off
 
       private
 
@@ -445,11 +444,12 @@
 C nlefsm_offsiteSO 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_offsiteSO( scell, nua, na, isa, xa, indxua,
+      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 , matrix_elements_only)
+     .                      fa, stress, H0 , H0_off,
+     &                      matrix_elements_only)
 
 
 C *********************************************************************
@@ -496,7 +496,8 @@
 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 H(maxnh,nspin)    : NL Hamiltonian (added to input H)
+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 *********************************************************************
@@ -525,7 +526,8 @@
 
       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) 
+      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
@@ -853,10 +855,10 @@
            ind = listhptr(iio)+j
            jo = listh(ind)
            H0(ind) = H0(ind) + Vi(jo)
-           H0_offsiteSO(ind,1) = H0_offsiteSO(ind,1) + V_so(1,1,jo)
-           H0_offsiteSO(ind,2) = H0_offsiteSO(ind,2) + V_so(2,2,jo)
-           H0_offsiteSO(ind,3) = H0_offsiteSO(ind,3) + V_so(1,2,jo)
-           H0_offsiteSO(ind,4) = H0_offsiteSO(ind,4) + V_so(2,1,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()
@@ -894,7 +896,7 @@
 
       call timer( 'nlefsm', 2 )
 
-      end subroutine nlefsm_offsiteSO
+      end subroutine nlefsm_SO_off
 
 c-----------------------------------------------------------------------
 c

=== modified file 'Src/setup_H0.F'
--- Src/setup_H0.F	2018-04-19 12:47:26 +0000
+++ Src/setup_H0.F	2018-04-23 11:10:59 +0000
@@ -17,11 +17,10 @@
       
       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 sparse_matrices, only: H0_offsiteSO
-      use m_nlefsm,        only: nlefsm_offsiteSO
+      use m_nlefsm,        only: nlefsm_SO_off
       use m_spin,          only: spin
 
       use sparse_matrices, only: listh, listhptr, numh, maxnh
@@ -44,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
@@ -57,13 +57,15 @@
       integer  :: ia, is
 
       real(dp)    :: dummy_Eso
-      integer     :: io, ispin, i, j 
-      complex(dp) :: Hc, Dc
+      integer     :: ispin, i, j 
+      complex(dp) :: Dc
 #ifdef MPI
-      real(dp)            :: buffer1
+      real(dp) :: buffer1
 #endif
 
-      real(dp), pointer :: H_val(:), H_so(:,:)
+      real(dp), pointer :: H_val(:), H_so_on(:,:)
+      complex(dp), pointer :: H_so_off(:,:)
+
 
 #ifdef DEBUG
       call write_debug( '    PRE setup_H0' )
@@ -136,13 +138,14 @@
      &             H_val,
      &             matrix_elements_only=.true.) 
       else
-        H0_offsiteSO=dcmplx(0.0d0,0.0d0)
-        call nlefsm_offsiteSO(scell, na_u, na_s, isa, xa, indxua,
+        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,
+     &                 dummy_stress, H_val, H_so_off,
      &                 matrix_elements_only=.true.)
 
 
@@ -153,24 +156,20 @@
 ! DM in a such way that the result gives Re{Tr[H_SO*DM]}.
 !
 
-        do io = 1, maxnh
+        do i = 1, maxnh
 
 !-------- Eso(u,u)
-            Dc = cmplx(Dscf(io,1),Dscf(io,5),kind=dp)
-            Hc = H0_offsiteSO(io,1)
-          Eso = Eso + real( Hc*Dc )
+          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(io,2),Dscf(io,6),kind=dp)
-            Hc = H0_offsiteSO(io,2)
-          Eso = Eso + real( Hc*Dc )
+          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(io,3),Dscf(io,4),kind=dp)
-            Hc = H0_offsiteSO(io,4)
-          Eso = Eso + real( Hc*Dc )
+          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(io,7),-Dscf(io,8),kind=dp)
-            Hc = H0_offsiteSO(io,3)
-          Eso = Eso + real( Hc*Dc )
+          Dc = cmplx(Dscf(i,7),-Dscf(i,8), dp)
+          Eso = Eso + real( H_so_off(i,3)*Dc, dp)
 
         enddo
 
@@ -188,12 +187,12 @@
 ! should be enough
 !
       if ( spin%SO_onsite ) then
-         H_so => val(H_so_2D)
+         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)
+     &        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

=== modified file 'Src/setup_hamiltonian.F'
--- Src/setup_hamiltonian.F	2018-04-20 09:41:55 +0000
+++ Src/setup_hamiltonian.F	2018-04-23 11:10:59 +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
@@ -42,7 +44,6 @@
       use m_ntm
 
       use m_spin,          only: spin 
-      use sparse_matrices, only: H0_offsiteSO
 
       use m_dipol
       use alloc, only: re_alloc, de_alloc
@@ -71,9 +72,11 @@
       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):: Hc, Dc
+      complex(dp):: Dc
       integer :: ind, i, j
 
 !------------------------------------------------------------------------- BEGIN
@@ -90,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_onsite ) then
-         ! Sadly some compilers (g95), does
-         ! not allow bounds for pointer assignments :(
-         H_so => val(H_so_2D)
+        ! 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
 
-      if ( .not. spin%SO_offsite ) then
-       do ispin = 1, spin%spinor
-         if (ispin .le. 2) then
-!$OMP do
-           do io = 1,maxnh
-             H(io,ispin) = H_kin(io) + H_vkb(io)
-           end do
-!$OMP end do nowait
-         else
-!$OMP do
-           do io = 1,maxnh
-             H(io,ispin) = 0.0_dp
-           end do
-!$OMP end do nowait
-         end if
-       end do
-      end if
+      ! 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
+!$OMP end do nowait
+      end do
 
       if ( spin%SO_onsite ) 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
-      else if ( spin%NCol ) then
-!$OMP do
-         do io = 1,maxnh
-            H(io,3) = 0._dp
-            H(io,4) = 0._dp
-         end do
-!$OMP end do nowait
+        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
          
 ! ..................
@@ -171,40 +174,36 @@
 ! DM in a such way that the result gives Re{Tr[H_SO*DM]}.
 !
 
-      if( spin%SO_offsite ) then
-
+      if ( spin%SO_offsite ) then
         do io = 1, maxnh
 
 !-------- Eso(u,u)
-            Dc = cmplx(Dscf(io,1),Dscf(io,5),kind=dp)
-            Hc = H0_offsiteSO(io,1)
-            Eso = Eso + real( Hc*Dc )
+          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),kind=dp)
-            Hc = H0_offsiteSO(io,2)
-            Eso = Eso + real( Hc*Dc )
+          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),kind=dp)
-            Hc = H0_offsiteSO(io,4)
-            Eso = Eso + real( Hc*Dc )
+          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),kind=dp)
-            Hc = H0_offsiteSO(io,3)
-            Eso = Eso + real( Hc*Dc )
-
-        enddo
-
-      endif
-
-      if ( spin%SO_onsite ) then
+          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
 
@@ -224,10 +223,6 @@
 !     Non-SCF part of total energy
       call update_E0()
 
-      if ( spin%SO_offsite ) then
-       H(:,:) = 0.0_dp
-      endif
-
 ! Hubbard term for LDA+U: energy, forces, stress and matrix elements ....
       if( switch_ldau ) then
         if ( spin%NCol ) then
@@ -285,24 +280,20 @@
 
       if ( spin%SO_offsite ) then
 
-        do is = 1 , spin%spinor
-          do io = 1, maxnh
-            H(io,is) = H(io,is) + H_kin(io) + H_vkb(io)
-          enddo
-        enddo
-
-C------- H(u,u)
-         H(:,1) = H(:,1) + real(H0_offsiteSO(:,1))
-         H(:,5) =          imag(H0_offsiteSO(:,1))
-C------- H(d,d)
-         H(:,2) = H(:,2) + real(H0_offsiteSO(:,2))
-         H(:,6) =          imag(H0_offsiteSO(:,2))
-C------- H(u,d)
-         H(:,3) = H(:,3) + real(H0_offsiteSO(:,3))
-         H(:,4) = H(:,4) + imag(H0_offsiteSO(:,3))
-C------- H(d,u)
-         H(:,7) = H(:,7) + real(H0_offsiteSO(:,4))
-         H(:,8) = H(:,8) - imag(H0_offsiteSO(:,4))
+! 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
 

=== modified file 'Src/sparse_matrices.F'
--- Src/sparse_matrices.F	2018-04-19 12:47:26 +0000
+++ Src/sparse_matrices.F	2018-04-23 11:10:59 +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
@@ -37,7 +38,10 @@
       ! 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
@@ -56,7 +60,6 @@
 
       subroutine resetSparseMatrices( )
       use alloc,  only : de_alloc
-      use m_spin, only : spin
 
       implicit none
 
@@ -67,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)
@@ -76,10 +80,6 @@
       call delete( xij_2D ) ; nullify(xijo)
 
       call delete( DM_history )
-
-      if( spin%SO_offsite ) then 
-         call de_alloc( H0_offsiteSO, 'H0_offsiteSO', 'sparseMat' )
-      endif
       
       ! Using MixH is bad as several utilities
       ! are not relying on FoX, better to leave out

=== modified file 'Src/state_init.F'
--- Src/state_init.F	2018-04-19 12:47:26 +0000
+++ Src/state_init.F	2018-04-23 11:10:59 +0000
@@ -26,10 +26,9 @@
       use sparse_matrices, only: xijo, xij_2D
       use sparse_matrices, only: S   , S_1D
 
-      use sparse_matrices, only: H0_offsiteSO
-
       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
@@ -100,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
@@ -564,15 +564,14 @@
          end if
       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_2D,name=oname)
-      end if
-
-      if ( spin%SO_offsite ) then
-         call re_alloc(H0_offsiteSO, 1,maxnh, 1,4, 'H0_offsiteSO',
-     $                 'state_init')
+      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

=== modified file 'version.info'
--- version.info	2018-04-23 09:26:04 +0000
+++ version.info	2018-04-23 11:10:59 +0000
@@ -1,2 +1,2 @@
-trunk-688--merge-OSSO-700
+trunk-688--merge-OSSO-700--class-1