Merge lp:~nickpapior/siesta/4.1-wfs into lp:~albertog/siesta/4.1-wfs

[Nick Papior]

Implemented tSpin in bands/writewave and fixed a few output options

Using tSpin cleans up a bit some of the *magic* variables used
here and there. Actually the use of tSpin suggests that we need
a new variable in tSpin to clarify whether it is NCol .or. SO
which could help in certain places.

The output files are now a bit weird when comparing format
for NC.
- ioeig uses no_u * 2, 1
- bands uses no_u, 4 for WFS file
- bands uses no_u * 2, 1 for .bands
- wwave uses no_u, 4 for WFSX

It is not clear to me why they are different, and/or
if they *should* be different.
I would prefer if we can use the same convention all-over.
This may require changes all-around.

I don't have any preference. Probably I would prefer:

   no_u, spin%H

which clarifies exactly which type of calculation it was.
This is probably opposite of my previous perception.

Also, I have cleaned up print_spin which did lots of
duplicate work for SO calculations.

Once this is in, I can approve the merge.

I think it would be best if we could clean up this once and for all.
Let me know what you think!

[Alberto Garcia (albertog)]

I agree with your tSpin and print_spin changes. They do make the code more clear. Thanks.

The new tSpin variable (equivalent to Ncol.or.SO), and other enhancements, can maybe be defined in the 4.2 branch.

Regarding the output files headings:

  In bands and ioeig, no_u*2 refers to the number of states, and 1 to the
  number of spin blocks (=1 since there is no spin distinction).

  In the wavefunction files, no_u refers to the number of basis orbitals, and
  4 to the number of expansions of length no_u needed at each band (4 for a complex
  spinor). The number of states in the file is given by nk and the following records
  for each k-point.

One possible source of fragility is that the "spinor_dim second-dimension" trick in ioeig and bands might not work when 'neigwanted' is not no_u. This could be also handled more robustly in 4.2, after the merge, and back-ported if deemed necessary for some use cases.

review: Approve

[Alberto Garcia (albertog)]

More on the bands/EIG/WFS headings:

I agree that they could be improved, but doing so at this time would delay the pressing merge of the NC/SO wfs functionality. Besides, just changing a small set of numbers by another small set of numbers is probably not going to do a lot towards self-documentation of the files. To get out of the historical constraints we probably need a more thorough rethink of the files' contents.

[Nick Papior (nickpapior)]

For EIG files etc. Agreed, it isn't totally necessary. I just wanted this on the table since it is rather confusing ;)
Once you know what they mean, it is easy to bypass the problems they yield.

=== modified file 'Src/Makefile'
--- Src/Makefile	2019-06-20 13:50:42 +0000
+++ Src/Makefile	2019-08-08 11:56:55 +0000
@@ -1346,7 +1346,7 @@
 write_subs.o: m_steps.o m_stress.o m_ts_global_vars.o parallel.o precision.o
 write_subs.o: siesta_cml.o siesta_geom.o siesta_options.o units.o zmatrix.o
 writewave.o: alloc.o atmfuncs.o atomlist.o densematrix.o diag.o diag_option.o
-writewave.o: get_kpoints_scale.o kpoint_grid.o m_spin.o parallel.o
+writewave.o: get_kpoints_scale.o kpoint_grid.o m_spin.o m_spin.o parallel.o
 writewave.o: parallelsubs.o precision.o siesta_geom.o sys.o units.o
 xml.o: precision.o
 zm_broyden_optim.o: m_broyddj_nocomm.o parallel.o precision.o sys.o units.o
 
=== modified file 'Src/bands.F'
--- Src/bands.F	2019-02-24 16:24:33 +0000
+++ Src/bands.F	2019-08-08 11:56:55 +0000
@@ -296,7 +296,7 @@
       end subroutine initbands
 
 
-      subroutine bands( no_s, h_spin_dim, spinor_dim,
+      subroutine bands( no_s, spin,
      .                  no_u, no_l,maxnh, maxk,
      .                  numh, listhptr, listh, H, S, ef, xij, indxuo,
      .                  writeb, nk, kpoint, ek, occtol, getPSI )
@@ -306,11 +306,11 @@
 C Written by J.Soler, August 1997 and August 1998.
 C Initialisation moved into a separate routine, JDG Jan 2000.
 C WFS options by A. Garcia, April 2012
+C Cleaned for tSpin by N. Papior, August 2019
 
 C **************************** INPUT **********************************
 C integer no_s                : Number of basis orbitals in supercell
-C integer h_spin_dim          : 
-C integer spinor_dim          : 
+C type(tSpin) spin            : Containing all spin-information
 C integer maxnh               : Maximum number of orbitals interacting
 C                               with any orbital
 C integer maxk                : Last dimension of kpoint and ek
@@ -320,7 +320,7 @@
 C                               hamiltonian matrix
 C integer listh(maxlh)        : Nonzero hamiltonian-matrix element
 C                               column indexes for each matrix row
-C real*8  H(maxnh,h_spin_dim)      : Hamiltonian in sparse form
+C real*8  H(maxnh,spin%H)     : Hamiltonian in sparse form
 C real*8  S(maxnh)            : Overlap in sparse form
 C real*8  ef                  : Fermi energy
 C real*8  xij(3,maxnh)        : Vectors between orbital centers (sparse)
@@ -338,7 +338,7 @@
 C real*8  kpoint(3,maxk)      : k point vectors
 C real*8  occtol              : Occupancy threshold for DM build
 C *************************** OUTPUT **********************************
-C real*8  ek(no_u,spinor_dim,maxk) : Eigenvalues
+C real*8  ek(no_u,spin%spinor,maxk) : Eigenvalues
 C *************************** UNITS ***********************************
 C Lengths in atomic units (Bohr).
 C k vectors in reciprocal atomic units.
@@ -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 t_spin, only: tSpin
 
       use m_diag_option, only: ParallelOverK, Serial
       use m_diag, only: diag_init
@@ -374,13 +374,13 @@
 
       implicit          none
 
-      integer  ::   h_spin_dim
-      integer  ::   maxk, maxnh, spinor_dim, no_u, no_l, nk, no_s,
+      type(tSpin), intent(in) :: spin
+      integer  ::   maxk, maxnh, no_u, no_l, nk, no_s,
      .              indxuo(no_s), listh(maxnh), numh(no_l),
      .              listhptr(*)
       logical  ::   writeb
-      real(dp) ::   ef, ek(no_u,spinor_dim,maxk),
-     .              H(maxnh,h_spin_dim), kpoint(3,maxk),
+      real(dp) ::   ef, ek(no_u,spin%spinor,maxk),
+     .              H(maxnh,spin%H), kpoint(3,maxk),
      .              S(maxnh), xij(3,maxnh), occtol
       logical, intent(in) :: getPSI
 
@@ -394,8 +394,7 @@
       logical, parameter :: fixspin = .false.
 
       integer :: ik, il, io, ispin, iu, iu_wfs, iuo, naux, nhs, j
-      integer :: nspin_flag  ! For wfs output
-      logical :: SaveParallelOverK, non_coll
+      logical :: SaveParallelOverK
       
       real(dp)
      .  Dnew, qs(2), e1, e2, efs(2), emax, emin, Enew, eV, qk, qtot,
@@ -424,15 +423,12 @@
 
 C     Allocate local arrays - only aux is relevant here
 
-      non_coll = (h_spin_dim >= 4)
-      if ( non_coll ) then
-         nspin_flag = 4
-         nhs  = 2 * (2*no_u) * (2*no_l)
+      if ( spin%Grid == 4 ) then
+         nhs = 2 * (2*no_u) * (2*no_l)
       else
-         nspin_flag = h_spin_dim
          nhs = 2 * no_u*no_l
       endif
-      naux  = 2*no_u*5
+      naux = 2*no_u*5
       call allocDenseMatrix(nhs, nhs, nhs)
       call re_alloc( aux, 1, naux, 'aux', 'bands' )
 
@@ -445,7 +441,7 @@
         rewind (iu_wfs)
 
         write(iu_wfs) nk, .false. ! nk, Gamma, same file-format in WFS as for Gamma-point
-        write(iu_wfs) nspin_flag
+        write(iu_wfs) spin%Grid
         write(iu_wfs) no_u
         write(iu_wfs) (iaorb(j),labelfis(isa(iaorb(j))),
      .            iphorb(j), cnfigfio(isa(iaorb(j)),iphorb(j)),
@@ -455,7 +451,7 @@
       endif
 
 C Find the band energies
-      if (NoMagn .or. SPpol) then
+      if ( spin%Grid < 4 ) then
 C fixspin and qs are not used in diagk, since getD=.false. ...
         qs(1) = 0.0_dp
         qs(2) = 0.0_dp
@@ -472,7 +468,7 @@
            call diag_init()
         end if
 
-        call diagk( h_spin_dim, no_l, no_s, spinor_dim, 
+        call diagk( spin%H, no_l, no_s, spin%spinor, 
      .              maxnh, maxnh, 
      .              no_u, numh, listhptr, listh, numh, listhptr,
      .              listh, H, S, getD, getPSI, fixspin, qtot, qs, temp,
@@ -484,7 +480,7 @@
         ParallelOverK = SaveParallelOverK
         if ( ParallelOverK ) Serial = .true.
 
-      elseif ( NonCol ) then
+      elseif ( spin%NCol ) then
         call diag2k(no_l, no_s, maxnh, maxnh, no_u,
      .              numh, listhptr, listh, numh, listhptr,
      .              listh, H, S, getD, getPSI, qtot, temp, e1, e2,
@@ -493,7 +489,7 @@
      .              Haux, Saux, psi, Haux, Saux, aux,
      .              no_u, occtol, 1, no_u )
 
-      elseif ( SpOrb ) then
+      elseif ( spin%SO ) then
         call diag3k(no_l, no_s, maxnh, maxnh, no_u, numh,
      .              listhptr, listh, numh, listhptr, listh,
      .              H, S, getD, getPSI, qtot, temp, e1, e2, xij,
@@ -502,7 +498,7 @@
      .              Haux, Saux, psi, Haux, Saux, aux,
      .              no_u, occtol, 1, no_u )
       else
-        call die( 'bands: ERROR: incorrect value of h_spin_dim')
+        call die( 'bands: ERROR: incorrect value of spin%H')
       endif
 
 C Write bands
@@ -525,7 +521,7 @@
      .      path = path + sqrt( (kpoint(1,ik)-kpoint(1,ik-1))**2 +
      .                          (kpoint(2,ik)-kpoint(2,ik-1))**2 +
      .                          (kpoint(3,ik)-kpoint(3,ik-1))**2 )
-          do ispin = 1,spinor_dim
+          do ispin = 1, spin%spinor
             do io = 1, no_u
               emax = max( emax, ek(io,ispin,ik) )
               emin = min( emin, ek(io,ispin,ik) )
@@ -537,10 +533,10 @@
         write(iu,*) emin/eV, emax/eV
 
 C Write eigenvalues
-        if ( non_coll) then
+        if ( spin%Grid == 4 ) then
           write(iu,*) 2*no_u, 1, nk     ! A single spin channel, double number of bands
         else 
-          write(iu,*) no_u, spinor_dim, nk
+          write(iu,*) no_u, spin%spinor, nk
         endif
 
         ! For non-collinear calculations, the ek array will contain the first no_u
@@ -554,11 +550,11 @@
      .                            (kpoint(3,ik)-kpoint(3,ik-1))**2 )
             write(iu,'(f10.6,10f12.4,/,(10x,10f12.4))')
      .         path, ((ek(io,ispin,ik)/eV,io=1,no_u),
-     .         ispin=1,spinor_dim)
+     .         ispin=1,spin%spinor)
           else
             write(iu,'(3f9.5,8f12.4,/,(27x,8f12.4))')
      .        kpoint(1,ik), kpoint(2,ik), kpoint(3,ik),
-     .        ((ek(io,ispin,ik)/eV,io=1,no_u),ispin=1,spinor_dim)
+     .        ((ek(io,ispin,ik)/eV,io=1,no_u),ispin=1,spin%spinor)
           endif
         enddo
 
 
=== modified file 'Src/print_spin.F90'
--- Src/print_spin.F90	2019-03-29 13:10:59 +0000
+++ Src/print_spin.F90	2019-08-08 11:56:55 +0000
@@ -2,7 +2,6 @@
   !
   ! Prints spin in output and CML files
   !
-  use m_spin,          only: nspin  ! This is spin%grid (1, 2, or 4)
   use m_spin,          only: spin
   use atomlist,        only: no_l
   use sparse_matrices, only: listhptr, numh
@@ -16,79 +15,73 @@
 
   implicit none
   
-  real(dp), intent(out)  :: qspin(4)
+  real(dp), intent(out)  :: qspin(spin%Grid)
 
   integer  :: ispin, io, j, ind
   real(dp) :: qaux
   real(dp) :: Stot        ! Total spin magnitude
   real(dp) :: Svec(3)     ! Total spin vector
-  real(dp) :: dm_aux(4)   ! Auxiliary array to hermitify Dscf for SOC
 #ifdef MPI
-  real(dp) :: qtmp(4)
+  real(dp) :: qtmp(spin%Grid)
 #endif
 
-  if (nspin < 2) RETURN
-
-  if (spin%SO) then
+  qspin(:) = 0.0_dp
+
+  if ( spin%Grid < 2 ) return
+
+  if ( spin%SO ) then
      
-     do ispin = 1,nspin
-        qspin(ispin) = 0.0_dp
-        do io = 1,no_l
-           do j = 1,numh(io)
-              ind = listhptr(io)+j
-              ! In the SOC case, hermitify Dscf
-              dm_aux(1:4) = dscf(ind,1:4)
-              dm_aux(3) = 0.5_dp*(dm_aux(3)+dscf(ind,7))
-              dm_aux(4) = 0.5_dp*(dm_aux(4)+dscf(ind,8))
-              qspin(ispin) = qspin(ispin) + dm_aux(ispin) * S(ind)
-           enddo
-        enddo
-     enddo
+    do io = 1,no_l
+      do j = 1,numh(io)
+        ind = listhptr(io)+j
+        ! In the SOC case, hermitify Dscf
+        qspin(1:2) = qspin(1:2) + dscf(ind,1:2) * S(ind)
+        qspin(3) = qspin(3) + 0.5_dp*(dscf(ind,3)+dscf(ind,7)) * S(ind)
+        qspin(4) = qspin(4) + 0.5_dp*(dscf(ind,4)+dscf(ind,8)) * S(ind)
+      end do
+    end do
 
   else
      
-     do ispin = 1,nspin
-        qspin(ispin) = 0.0_dp
-        do io = 1,no_l
-           do j = 1,numh(io)
-              ind = listhptr(io)+j
-              qspin(ispin) = qspin(ispin) + Dscf(ind,ispin) * S(ind)
-           enddo
-        enddo
-     enddo
+    do io = 1,no_l
+      do j = 1,numh(io)
+        ind = listhptr(io)+j
+        qspin(:) = qspin(:) + Dscf(ind,:) * S(ind)
+      end do
+    end do
      
   endif
 #ifdef MPI
-     !       Global reduction of spin components
-     call globalize_sum(qspin(1:nspin),qtmp(1:nspin))
-     qspin(1:nspin) = qtmp(1:nspin)
+  ! Global reduction of spin components
+  call globalize_sum(qspin(1:spin%Grid),qtmp(1:spin%Grid))
+  qspin(1:spin%Grid) = qtmp(1:spin%Grid)
 #endif
 
-     if (nspin .eq. 2) then
-
-        if (IOnode) then
-           Svec(1) = 0.0_dp
-           Svec(2) = 0.0_dp
-           Svec(3) = qspin(1) - qspin(2)
-           Stot = Svec(3)
-           write(6,'(5x,a,f10.5,2f10.1,f10.5)') 'spin moment: S , {S} = ', Stot, Svec
-        endif
-        if (cml_p) call cmlAddProperty(xf=mainXML,            &
-             value=qspin(1)-qspin(2), dictref='siesta:stot', &
-             units='siestaUnits:spin')
-
-     elseif (nspin .eq. 4) then
-
-        call spnvec( nspin, qspin, qaux, Stot, Svec )
-        if (IONode) then
-           write(6,'(5x,a,4f10.5)') 'spin moment: S , {S} = ', Stot, Svec
-           if (cml_p) then
-              call cmlAddProperty(xf=mainXML, value=Stot,  &
-                   dictref='siesta:stot', units='siestaUnits:spin')
-              call cmlAddProperty(xf=mainXML, value=Svec,  &
-                   dictref='siesta:svec', units='siestaUnits:spin')
-           endif !cml_p
-        endif
-     endif
+  ! We are left with printing out to stdout
+  if ( .not. IONode ) return
+  
+  if ( spin%Grid == 2) then
+
+    Svec(1) = 0.0_dp
+    Svec(2) = 0.0_dp
+    Svec(3) = qspin(1) - qspin(2)
+    Stot = Svec(3)
+    write(6,'(5x,a,f10.5,2f10.1,f10.5)') 'spin moment: S , {S} = ', Stot, Svec
+    if (cml_p) call cmlAddProperty(xf=mainXML,            &
+        value=qspin(1)-qspin(2), dictref='siesta:stot', &
+        units='siestaUnits:spin')
+    
+  else if ( spin%Grid == 4 ) then
+
+    call spnvec( spin%Grid, qspin, qaux, Stot, Svec )
+    write(6,'(5x,a,4f10.5)') 'spin moment: S , {S} = ', Stot, Svec
+    if (cml_p) then
+      call cmlAddProperty(xf=mainXML, value=Stot,  &
+          dictref='siesta:stot', units='siestaUnits:spin')
+      call cmlAddProperty(xf=mainXML, value=Svec,  &
+          dictref='siesta:svec', units='siestaUnits:spin')
+    end if !cml_p
+
+  end if
 
 end subroutine print_spin
 
=== modified file 'Src/siesta_analysis.F'
--- Src/siesta_analysis.F	2019-02-27 08:38:58 +0000
+++ Src/siesta_analysis.F	2019-08-08 11:56:55 +0000
@@ -46,8 +46,7 @@
       use m_energies
       use m_steps, only: final
       use m_ntm
-      use m_spin,         only: nspin, spinor_dim, h_spin_dim
-      use m_spin,         only: SpOrb, NonCol, SPpol, NoMagn
+      use m_spin,         only: spin
       use m_dipol
       use m_eo
       use m_forces,          only: fa
@@ -68,10 +67,10 @@
 
       real(dp)          :: dummy_str(3,3)
       real(dp)          :: dummy_strl(3,3)
-      real(dp)          :: qspin(4)         ! Local
+      real(dp)          :: qspin(spin%Grid)         ! Local
 
-      real(dp)          :: polxyz(3,nspin)  ! Autom., small
-      real(dp)          :: polR(3,nspin)    ! Autom., small
+      real(dp)          :: polxyz(3,spin%Grid)  ! Autom., small
+      real(dp)          :: polR(3,spin%Grid)    ! Autom., small
       real(dp)          :: qaux
       real(dp), pointer :: ebk(:,:,:)       ! Band energies
 
@@ -127,7 +126,7 @@
 
       if (fdf_get("Read-H-from-file",.false.)) then
          if (SIESTA_worker) then
-            call read_spmatrix(maxnh, no_l, h_spin_dim, numh,
+            call read_spmatrix(maxnh, no_l, spin%H, numh,
      .           listhptr, listh, H, found, userfile="H_IN")
             if (.not. found) call die("Could not find H_IN")
             current_ef = ef
@@ -137,7 +136,7 @@
 
 #ifdef SIESTA__PEXSI
       if (fdf_get("PEXSI.DOS",.false.)) then
-         call pexsi_dos(no_u, no_l, spinor_dim,
+         call pexsi_dos(no_u, no_l, spin%spinor,
      $        maxnh, numh, listhptr, listh, H, S, qtot, ef)
          
       endif
@@ -263,6 +262,17 @@
 !        latest DM, although this is probably too much).
 !        See the logic in siesta_forces.
 
+!     Setup the number of states depending on the used method
+!     For polarized calculations we use no_u and separate spin-channels.
+!     For NC/SO we have mixed spin-channels and separating makes no
+!     sense.
+      if ( spin%Grid == 4 ) then
+!       We are NC/SO
+        max_n_states = no_u * 2
+      else
+        max_n_states = no_u
+      end if
+
 !     Find and print wavefunctions at selected k-points
 !**   This uses H,S, and xa
       if (nwk.gt.0) then
@@ -270,7 +280,7 @@
         if (IONode) print "(a)",
      $              "Writing WFSX for selected k-points in "
      $                           // trim(wfs_filename)
-        call wwave( no_s, h_spin_dim, spinor_dim, no_u, no_l, maxnh, 
+        call wwave( no_s, spin, no_u, no_l, maxnh, 
      &              nwk,
      &              numh, listhptr, listh, H, S, Ef, xijo, indxuo,
      &              gamma_wavefunctions, nwk, wfk, no_u, occtol )
@@ -278,13 +288,10 @@
 
 
 !**   This uses H,S, and xa
-      if (write_coop) then
+      if ( write_coop ) then
         ! Output the wavefunctions for the kpoints in the SCF set
         ! Note that we allow both a band number and an energy range
         ! The user is responsible for using appropriate values.
-         max_n_states = no_u
-         if (h_spin_dim >= 4) max_n_states = 2*no_u
-         
         wfs_band_min = fdf_get("WFS.BandMin",1)
         wfs_band_max = fdf_get("WFS.BandMax",max_n_states)
         call setup_wfs_list(nkpnt,max_n_states,
@@ -294,7 +301,7 @@
          wfs_filename = trim(slabel)//".fullBZ.WFSX"
          if (IONode) print "(a)", "Writing WFSX for COOP/COHP in " 
      $                           // trim(wfs_filename)
-         call wwave( no_s, h_spin_dim, spinor_dim, no_u, no_l, maxnh, 
+         call wwave( no_s, spin, no_u, no_l, maxnh, 
      .        nkpnt,
      .        numh, listhptr, listh, H, S, Ef, xijo, indxuo,
      .        gamma_SCF, nkpnt, kpoint, no_u, occtol)
@@ -303,85 +310,101 @@
 !     Compute bands
 !**   This uses H,S, and xa
       nullify( ebk )
-      call re_alloc( ebk, 1, no_u, 1, spinor_dim, 1, maxbk,
-     &               'ebk', 'siesta_analysis' )
-
-      if (nbk.gt.0) then
+      if ( spin%Grid == 4 ) then
+        call re_alloc( ebk, 1, no_u*2, 1, 1, 1, maxbk,
+     &      'ebk', 'siesta_analysis' )
+      else
+        call re_alloc( ebk, 1, no_u, 1, spin%spinor, 1, maxbk,
+     &      'ebk', 'siesta_analysis' )
+      end if
+      
+      if ( nbk > 0 ) then
         if (IONode) print "(a)", "Computing bands..."
         getPSI = fdf_get('WFS.Write.For.Bands', .false.)
         if (getPSI) then
-           wfs_filename = trim(slabel)//".bands.WFSX"
-         if (IONode) print "(a)", "Writing WFSX for bands in "
-     $                           // trim(wfs_filename)
-         max_n_states = no_u
-         if (h_spin_dim >= 4) max_n_states = 2*no_u
-           wfs_band_min = fdf_get("WFS.BandMin",1)
-           wfs_band_max = fdf_get("WFS.BandMax",max_n_states)
-           call setup_wfs_list(nbk,max_n_states,
-     $          wfs_band_min,wfs_band_max,
-     $          use_scf_weights=.false.,
-     $          use_energy_window=.false.)
-        endif
-        call bands( no_s, h_spin_dim, spinor_dim, no_u, no_l, maxnh, 
+          wfs_filename = trim(slabel)//".bands.WFSX"
+          if (IONode) print "(a)", "Writing WFSX for bands in "
+     $        // trim(wfs_filename)
+          wfs_band_min = fdf_get("WFS.BandMin",1)
+          wfs_band_max = fdf_get("WFS.BandMax",max_n_states)
+          call setup_wfs_list(nbk,max_n_states,
+     $        wfs_band_min,wfs_band_max,
+     $        use_scf_weights=.false.,
+     $        use_energy_window=.false.)
+        end if
+        call bands( no_s, spin, no_u, no_l, maxnh, 
      .              maxbk,
      .              numh, listhptr, listh, H, S, Ef, xijo, indxuo,
-     .              .true., nbk, bk, ebk, occtol, getPSI )
-        if (IOnode) then
+     .      .true., nbk, bk, ebk, occtol, getPSI )
+
+        if ( IOnode ) then
           if ( writbk ) then
-            write(6,'(/,a,/,a4,a12)')
-     &       'siesta: Band k vectors (Bohr**-1):', 'ik', 'k'
+            write(6,'(/,a,/,a4,tr1,a12)')
+     &          'siesta: Band k vectors (Bohr**-1):', 'ik', 'k'
             do ik = 1,nbk
-              write(6,'(i4,3f12.6)') ik, (bk(ix,ik),ix=1,3)
-            enddo
-          endif
+              write(6,'(i4,3(tr1,f12.6))') ik, (bk(ix,ik),ix=1,3)
+            end do
+          end if
         
           if ( writb ) then
-            write(6,'(/,a,/,a4,a3,a7)')
-     &       'siesta: Band energies (eV):', 'ik', 'is', 'eps'
-            do ispin = 1, spinor_dim
+            if ( spin%Grid == 4 ) then
+              write(6,'(/,a,/,a4,tr1,a9)')
+     &            'siesta: Band energies (eV):', 'ik', 'eps'
               do ik = 1,nbk
-                write(6,'(i4,i3,10f7.2)')
-     &            ik, ispin, (ebk(io,ispin,ik)/eV,io=1,min(10,no_u))
-                if (no_u.gt.10) write(6,'(7x,10f7.2)')
-     &              (ebk(io,ispin,ik)/eV,io=11,no_u)
+                write(6,'(i4,10(tr1,f9.4))')
+     &              ik, (ebk(io,1,ik)/eV,io=1,min(10,2*no_u))
+                if ( 2*no_u > 10 ) write(6,'(4x,10(tr1,f9.4))')
+     &              (ebk(io,1,ik)/eV,io=11,no_u*2)
+              end do
+
+            else
+              write(6,'(/,a,/,a4,a3,tr1,a9)')
+     &            'siesta: Band energies (eV):', 'ik', 'is', 'eps'
+              do ispin = 1, spin%spinor
+                do ik = 1,nbk
+                  write(6,'(i4,i3,10(tr1,f9.4))')
+     &                ik, ispin, (ebk(io,ispin,ik)/eV,io=1,min(10,no_u))
+                  if ( no_u > 10 ) write(6,'(7x,10(tr1,f9.4))')
+     &                (ebk(io,ispin,ik)/eV,io=11,no_u)
+                enddo
               enddo
-            enddo
+            end if
           endif
+          
         endif
       endif
 
 !     Print eigenvalues
-      if (IOnode .and. writeig) then
-        if ((isolve.eq.SOLVE_DIAGON .or.
-     &       ((isolve.eq.SOLVE_MINIM) .and. minim_calc_eigenvalues))
-     &      .and. no_l.lt.1000) then
-          if ( h_spin_dim <= 2 ) then
-            write(6,'(/,a,/,a4,a3,a7)')
-     &       'siesta: Eigenvalues (eV):', 'ik', 'is', 'eps'
-            do ik = 1,nkpnt
-              do ispin = 1,spinor_dim
-                write(6,'(i4,i3,10f7.2)')
+      if ( IOnode .and. writeig .and. no_l < 1000 .and. 
+     &    (isolve.eq.SOLVE_DIAGON .or.
+     &    (isolve.eq.SOLVE_MINIM .and. minim_calc_eigenvalues)) ) then
+
+        if ( spin%Grid == 4 ) then
+          write(6,'(/,a)') 'siesta: Eigenvalues (eV):'
+          do ik = 1,nkpnt
+            write(6,'(a,i6)') 'ik =', ik
+            write(6,'(10(tr1,f9.4))') (eo(io,1,ik)/eV,io=1,2*neigwanted)
+          end do
+
+        else
+          write(6,'(/,a,/,a4,a3,tr1,a9)')
+     &        'siesta: Eigenvalues (eV):', 'ik', 'is', 'eps'
+          do ik = 1,nkpnt
+            do ispin = 1, spin%spinor
+              write(6,'(i4,i3,10(tr1,f9.4))')
      &            ik,ispin,(eo(io,ispin,ik)/eV,io=1,min(10,neigwanted))
-                if (no_u.gt.10) write(6,'(7x,10f7.2)')
+              if (no_u.gt.10) write(6,'(7x,10(tr1,f9.4))')
      &            (eo(io,ispin,ik)/eV,io=11,neigwanted)
-              enddo
-            enddo
-          else
-            write(6,'(/,a)') 'siesta: Eigenvalues (eV):'
-            do ik = 1,nkpnt
-              write(6,'(a,i6)') 'ik =', ik
-              write(6,'(10f7.2)')
-     &          ((eo(io,ispin,ik)/eV,io=1,neigwanted),ispin=1,2)
-            enddo
-          endif
-          write(6,'(a,f15.6,a)') 'siesta: Fermi energy =', ef/eV, ' eV'
-        endif
+            end do
+          end do
+        end if
+        write(6,'(a,f15.6,a)') 'siesta: Fermi energy =', ef/eV, ' eV'
       endif
 
       if (((isolve.eq.SOLVE_DIAGON).or.
      &     ((isolve.eq.SOLVE_MINIM).and.minim_calc_eigenvalues))
      &     .and.IOnode)
-     &     call ioeig(eo,ef,neigwanted,nspin,nkpnt,no_u,spinor_dim,
+     &     call ioeig(eo,ef,neigwanted,spin%H,nkpnt,no_u,spin%spinor,
      &     nkpnt,kpoint, kweight)
 
 
@@ -437,7 +460,7 @@
 !     Calculation of the bulk polarization using the Berry phase
 !     formulas by King-Smith and Vanderbilt
       if (nkpol.gt.0 .and. .not.bornz) then
-        if (NonCol .or. SpOrb) then
+        if ( spin%NCol .or. spin%SO ) then
           if (IOnode) then
             write(6,'(/a)')
      .       'siesta_analysis: bulk polarization implemented only for'
@@ -446,7 +469,7 @@
           endif
         else
          call KSV_pol(na_u, na_s, xa_last, rmaxo, scell_last, 
-     &               ucell_last,no_u, no_l, no_s, nspin, qspin, 
+     &               ucell_last,no_u, no_l, no_s, spin%Grid, qspin, 
      &               maxnh, nkpol, numh, listhptr, listh, 
      &               H, S, xijo, indxuo, isa, iphorb, 
      &               iaorb, lasto, shape,
@@ -456,7 +479,7 @@
 
 !     Calculation of the optical conductivity
       call optical(na_u, na_s, xa_last, scell_last, ucell_last,
-     &             no_u, no_l, no_s, nspin, qspin,
+     &             no_u, no_l, no_s, spin%Grid, qspin,
      &             maxnh, numh, listhptr, listh, H, S, xijo,
      $             indxuo, ebk, ef, temp,
      &             isa, iphorb, iphKB, lasto, lastkb, shape )
@@ -486,7 +509,7 @@
         g2max = g2cut
         dummy_str  = 0.0
         dummy_strl = 0.0
-        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_last, indxua, 
      &              ntm, 0, 0, 0, filesOut, 
      &              maxnh, numh, listhptr, listh, Dscf, Datm,
 
=== modified file 'Src/writewave.F'
--- Src/writewave.F	2019-06-17 08:37:13 +0000
+++ Src/writewave.F	2019-08-08 11:56:55 +0000
@@ -40,7 +40,7 @@
       USE alloc, only: re_alloc
       USE sys, only: die
       USE atomlist, only: no_u
-      USE m_spin, only: h_spin_dim
+      USE m_spin, only: spin
 
       implicit none
 
@@ -51,11 +51,11 @@
       nullify(wfk)
 
       nwk = 0
-      if (h_spin_dim < 4) then
-         max_n_states = no_u
+      if ( spin%Grid == 4 ) then
+        max_n_states = 2 * no_u
       else
-         max_n_states = 2*no_u
-      endif
+        max_n_states = no_u
+      end if
       call initwave( max_n_states, nwk, wfk )
 
       if (nwk .eq. 0) then
@@ -391,7 +391,7 @@
 
       end subroutine initwave
 
-      subroutine wwave( no, nspin, maxspn, maxo, maxuo, maxnh, 
+      subroutine wwave( no, spin, maxo, maxuo, maxnh, 
      .                  maxk,
      .                  numh, listhptr, listh, H, S, ef, xij, indxuo,
      .                  gamma, nk, kpoint, nuotot, occtol)
@@ -399,10 +399,10 @@
 C Finds wavefunctions at selected k-points.
 C Written by P. Ordejon, June 2003
 C from routine 'bands' written by J.M.Soler
+C Changed to use tSpin, N. Papior, August 2019
 C **************************** INPUT **********************************
 C integer no                  : Number of basis orbitals
-C integer nspin  = h_spin_dim : Number of spin components of H, D
-C integer maxspn = spinor_dim : Second dimension of ek
+C tSpin spin                  : Contains spin information
 C integer maxo                : First dimension of ek
 C integer maxuo               : Second dimension of H and S
 C integer maxnh               : Maximum number of orbitals interacting  
@@ -414,7 +414,7 @@
 C                               hamiltonian matrix
 C integer listh(maxlh)        : Nonzero hamiltonian-matrix element  
 C                               column indexes for each matrix row
-C real*8  H(maxnh,nspin)      : Hamiltonian in sparse form
+C real*8  H(maxnh,spin%H)     : Hamiltonian in sparse form
 C real*8  S(maxnh)            : Overlap in sparse form
 C real*8  ef                  : Fermi energy
 C real*8  xij(3,maxnh)        : Vectors between orbital centers (sparse)
@@ -446,7 +446,7 @@
       use atmfuncs,     only : symfio, cnfigfio, labelfis, nofis
       use atomlist,     only : iaorb, iphorb
       use siesta_geom,  only : isa
-      use m_spin,       only : NoMagn, SPpol, NonCol, SpOrb
+      use t_spin, only: tSpin
       use units, only : eV
 
 #ifdef MPI
@@ -457,12 +457,12 @@
 
       implicit          none
 
+      type(tSpin), intent(in) :: spin
       logical, intent(in) :: Gamma
       integer           maxk, maxnh, maxo, maxuo, nk, no, 
-     .                  h_spin_dim, spinor_dim, nspin, maxspn,
      .                  nuotot, indxuo(no), listh(maxnh), numh(*), 
      .                  listhptr(*)
-      real(dp)          ef, H(:,:), kpoint(3,maxk), 
+      real(dp)          ef, H(maxnh,spin%H), kpoint(3,maxk), 
      .                  S(maxnh), xij(3,maxnh), occtol
 
       external          io_assign, io_close, memory
@@ -488,11 +488,7 @@
       data Dnew, Enew, e1, e2, qk, qtot, temp, wk /8*0.d0/
 
       logical :: SaveParallelOverK
-      integer :: nspin_flag
       
-      h_spin_dim=size(H,dim=2)
-      nspin_flag = h_spin_dim
-      if (h_spin_dim == 8) nspin_flag = 4
 
 C Get local number of orbitals
 #ifdef MPI
@@ -517,15 +513,15 @@
 C Check spin
 
 C Allocate local arrays - only aux is relevant here
-      if ( h_spin_dim >= 4 ) then
+      if ( spin%Grid == 4 ) then
          nhs = 2 * (2*nuotot) * (2*nuo)
       else
          nhs = 2 * nuotot * nuo
       endif
       call allocDenseMatrix(nhs, nhs, nhs) 
 
-      allocate(ek(nuotot,maxspn,nk))
-      call memory('A','D',maxspn*nk*nuotot,'writewave')
+      allocate(ek(nuotot,spin%spinor,nk))
+      call memory('A','D',nuotot*spin%spinor*nk,'writewave')
       naux = 2*nuotot*5
       allocate(aux(naux))
       call memory('A','D',naux,'writewave')
@@ -544,12 +540,12 @@
           write(6,'(a)') 'writewave: Wave Functions Coefficients'
           write(6,*)
           write(6,'(a26,2x,i6)') 'Number of k-points = ', nk
-          write(6,'(a26,2x,i6)') 'Number of Spins = ', nspin_flag
+          write(6,'(a26,2x,i6)') 'Number of Spins = ', spin%Grid
           write(6,'(a26,2x,i6)') 'Number of basis orbs = ',nuotot
           write(6,*)
         endif
         write(iu) nk, gamma
-        write(iu) nspin_flag   ! This will be 1, 2, or 4
+        write(iu) spin%Grid   ! This will be 1, 2, or 4
         write(iu) nuotot
         write(iu) (iaorb(j),labelfis(isa(iaorb(j))),
      .            iphorb(j), cnfigfio(isa(iaorb(j)),iphorb(j)),
@@ -567,9 +563,9 @@
 
 C Call appropriate diagonalization routine
 
-      if (NoMagn .or. SPpol) then
+      if ( spin%Grid <= 2 ) then
        if (gamma) then
-         call diagg( h_spin_dim, nuo, maxuo, maxnh, maxnh, 
+         call diagg( spin%H, nuo, maxuo, maxnh, maxnh, 
      .               maxo, numh, listhptr, listh, numh, listhptr, 
      .               listh, H, S, getD, getPSI,
      .               fixspin, qtot, qs, temp,
@@ -589,7 +585,7 @@
              call diag_init()
           end if
 
-         call diagk( h_spin_dim, nuo, no, maxspn, maxnh, maxnh,
+         call diagk( spin%H, nuo, no, spin%spinor, maxnh, maxnh,
      .               maxo, numh, listhptr, listh, numh, listhptr, 
      .               listh, H, S, getD, getPSI,
      .               fixspin, qtot, qs, temp,
@@ -603,7 +599,7 @@
          
        endif
 
-      elseif (NonCol) then
+      elseif ( spin%NCol) then
         if (gamma) then
           call diag2g( nuo, no, maxnh, maxnh, maxo, numh,
      .                 listhptr, listh, numh, listhptr, listh,
@@ -620,7 +616,7 @@
      .                 Haux, Saux, psi, Haux, Saux, aux,
      .                 nuotot, occtol, 1, nuotot )
         endif
-      elseif (SpOrb) then
+      elseif ( spin%SO ) then
         if (gamma) then
           call diag3g( nuo, no, maxnh, maxnh, maxo, numh,
      .                 listhptr, listh, numh, listhptr, listh,
 
=== modified file 'version.info'
--- version.info	2019-06-20 13:50:42 +0000
+++ version.info	2019-08-08 11:56:55 +0000
@@ -1,1 +1,1 @@
-siesta-4.1-1074--wfs-1115
+siesta-4.1-1074--wfs-1115--nicpa-1