Merge lp:~nickpapior/siesta/4.1-streamline-file-dim-nc-soc into lp:siesta/4.1

[Nick Papior]

Changed spin-index integers from NC/SOC == 4, to NC == 4, SOC == 8

This is to streamline the files across and enables utilities
to distinguish between NC and SOC calculations.
Such distinguishing is generally not needed, but it may be
useful in some cases.

Two test cases are added:
- bi2se3-stm-denchar
- bi2se3-sts-denchar

where I have tested this against 4.1 for stm and denchar
executable.

Some of the files contains the nspin index which generally was limited to 4.

However, to ensure an equal file basis across the files we now allow nspin == 8 for SOC
calculations.

This is to enable utilities to differentiate between NC and SOC calculations which may be
useful in some situations.

Out of consistency it is better to unify this across all files to ensure correct handling.

[Alberto Garcia (albertog)]

I agree with most changes, except those mentioned in the comments below. One is cosmetic, the other I think it leads to wrong behavior.

Beyond these, I have some other cosmetic comments in the code:

- Change 'adviced' to 'advised' (this comes from earlier commits)
- In the new output of Eig2DOS, a message appropriate for nspin=1 is missing. It would print
        "# Eigenvalues calculated from a spin-polarized calculation"

review: Needs Fixing

[lp:~nickpapior/siesta/4.1-streamline-file-dim-nc-soc]

1121. By Nick Papior on 2019-09-26

Backported a change in Denchar-wavofr which was wrong

Also fixed NC/SOC NumberOfEigenstates outputs for eigenvalues.
This simply required a streamlining of data in the ioeig routine
which was very simple.

[Nick Papior (nickpapior)]

I have backported the wrong-doings!

1. wavofr now is the same as before (please do double check)
2. I have reverted Eig2DOS and ioeig to write out bands. I also added a comment
3. I have fixed the NumberOfEigenstates output for ioeig in case of NC/SOC calculations and tested this on the mix_broyden test (NC).

I will fix adviced on the 4.1 branch (I can only find it in unrelated codes).

[Nick Papior (nickpapior)]

Added comments to comments ;)

[Alberto Garcia (albertog)]

OK.

I have added a few more comments around the "ravel" call in ioeig.

Also, I have included a check for nspin>2 in the wfs2wfsx converters, as they will not work in that case. There should be no WFS or WFS.nc files with nspin>2...

The changes are in my working directory. I will do the merge directly myself after pulling later
changes to 4.1. That is why I say "comment only" now.

[Alberto Garcia (albertog)]

Approved and merged.

review: Approve

=== modified file 'Src/ioeig.f'
--- Src/ioeig.f	2019-02-25 14:21:14 +0000
+++ Src/ioeig.f	2019-09-26 07:10:24 +0000
@@ -20,13 +20,13 @@
       implicit          none
 
       integer,  intent(in) :: no     ! no_u: number of orbitals in unit cell
-      integer,  intent(in) :: nspin  ! 'nspin_grid': 1, 2, or 4
+      integer,  intent(in) :: nspin  ! 'nspin_grid': 1, 2, 4 or 8
       integer,  intent(in) :: nk
       integer,  intent(in) :: maxo   ! no_u again
       integer,  intent(in) :: nspinor
       integer,  intent(in) :: maxk
       real(dp), intent(in) :: ef
-      real(dp), intent(in) :: eo(maxo, nspinor, maxk)
+      real(dp), intent(in), target :: eo(maxo, nspinor, maxk)
       real(dp), intent(in) :: kpoints(3,nk)
       real(dp), intent(in) :: kweights(nk)
       
@@ -34,6 +34,7 @@
 
 c Internal 
       integer           ik, iu, io, is
+      real(dp), pointer :: eok(:)
 
       character(len=label_length+4) :: fname
 c -------------------------------------------------------------------
@@ -44,20 +45,25 @@
       open( iu, file=fname, form='formatted', status='unknown' )      
 
       write(iu,"(e17.9)") ef/eV
+      ! The output corresponds to the number of bands.
+      ! Thus it shouldn't be confused with the number of orbitals,
+      ! although they are related!
       if ( nspin > nspinor ) then
-        ! NC/SOC case: a single block of wavefunctions (without spin tag)
-        ! nspin (4) will be used here to signal NC/SOC to utilities such as 
-        ! Eig2DOS
-        write(iu,"(tr1,i10,i2,tr1,i10)")   no*2, nspin, nk
+        write(iu,"(tr1,i10,tr1,i0,tr1,i10)") 2*no, nspin, nk
       else
-        write(iu,"(tr1,i10,i2,tr1,i10)")   no, min(nspin,2), nk
+        ! This will always be nspin == nspinor with max(nspinor) == 2
+        write(iu,"(tr1,i10,tr1,i0,tr1,i10)") no, nspin, nk
       end if
       do ik = 1,nk
-        ! This statement will correctly print eo in the NC/SOC case,
-        ! as the eigenvalues take up two complete 'io' columns,
-        ! and nspinor = 2. 
-        write(iu,"(i10,10(tr1,e17.9),/,(tr10,10(tr1,e17.9)))")
-     .          ik, ((eo(io,is,ik)/eV,io=1,no),is=1,nspinor)
+        if ( nspin > nspinor ) then
+          ! ensure we catch users doing neigwanted calculations
+          call ravel(maxo * nspinor, eo(1,1,ik), eok)
+          write(iu,"(i10,10(tr1,e17.9),/,(tr10,10(tr1,e17.9)))")
+     .        ik, (eok(io)/eV,io=1,no*nspinor)
+        else
+          write(iu,"(i10,10(tr1,e17.9),/,(tr10,10(tr1,e17.9)))")
+     .        ik, ((eo(io,is,ik)/eV,io=1,no),is=1,nspinor)
+        end if
       enddo
 
       call io_close( iu )
@@ -74,12 +80,13 @@
            call cmlStartPropertyList(mainXML, dictRef='siesta:kpt_band')
            do ik = 1, nk
              call cmlAddKPoint(xf=mainXML, coords=kpoints(:, ik), 
-     .             weight=kweights(ik))
+     .           weight=kweights(ik))
+             ! ensure we catch users doing neigwanted calculations
+             call ravel(maxo * nspinor, eo(1,1,ik), eok)
              call cmlAddProperty(xf=mainXML,
-     .           value=reshape(eo(1:no,1:nspinor,ik)/eV,
-     .           (/no*nspinor/)),
-     .            dictRef='siesta:eigenenergies',
-     .            units='siestaUnits:ev')
+     .           value=eok(no*nspinor)/eV,
+     .           dictRef='siesta:eigenenergies',
+     .           units='siestaUnits:ev')
            end do
            call cmlEndPropertyList(mainXML)
          else
@@ -105,5 +112,14 @@
          end if
          call cmlEndPropertyList(mainXML)
       endif
+
+      contains
+
+      subroutine ravel(n, eo, eop)
+      integer, intent(in) :: n
+      real(dp), intent(in), target :: eo(n)
+      real(dp), pointer :: eop(:)
+      eop => eo(:)
+      end subroutine
       
       end subroutine ioeig
 
=== modified file 'Src/writewave.F'
--- Src/writewave.F	2019-08-08 11:55:07 +0000
+++ Src/writewave.F	2019-09-26 07:10:24 +0000
@@ -540,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 = ', spin%Grid
+          write(6,'(a26,2x,i6)') 'Number of Spins = ', spin%H
           write(6,'(a26,2x,i6)') 'Number of basis orbs = ',nuotot
           write(6,*)
         endif
         write(iu) nk, gamma
-        write(iu) spin%Grid   ! This will be 1, 2, or 4
+        write(iu) spin%H   ! This will be 1, 2, 4 or 8 (to signal NC from SOC)
         write(iu) nuotot
         write(iu) (iaorb(j),labelfis(isa(iaorb(j))),
      .            iphorb(j), cnfigfio(isa(iaorb(j)),iphorb(j)),
@@ -738,7 +738,7 @@
           write(6,*)
           write(6,'(a22,2x,i6,2x,3f10.6)') 'k-point = ',ik, k(1:3)
           if (non_coll) then
-             write(6,'(a22,2x,i6)') 'Spin component (non-collinear)= ',0
+             write(6,'(a22,2x)') 'Spinors are mixed'
           else
              write(6,'(a22,2x,i6)') 'Spin component = ',ispin
           endif
 
=== added directory 'Tests/bi2se3-stm-denchar'
=== added file 'Tests/bi2se3-stm-denchar/Makefile'
--- Tests/bi2se3-stm-denchar/Makefile	1970-01-01 00:00:00 +0000
+++ Tests/bi2se3-stm-denchar/Makefile	2019-09-26 07:10:24 +0000
@@ -0,0 +1,6 @@
+#
+# Single-test makefile
+#
+name=bi2se3-stm-denchar
+#
+include ../test.mk
 
=== added file 'Tests/bi2se3-stm-denchar/README'
--- Tests/bi2se3-stm-denchar/README	1970-01-01 00:00:00 +0000
+++ Tests/bi2se3-stm-denchar/README	2019-09-26 07:10:24 +0000
@@ -0,0 +1,9 @@
+This example also involves denchar and stm.
+
+To complete this example you also need to
+run the following things:
+
+1. Go into work directory
+2. Create a soft-link using: ln -s bi2se3-stm-denchar.bands.WFSX bi2se3-stm-denchar.WFSX
+3. Run denchar < ../bi2se3-stm-denchar.fdf
+4. Run stm < ../bi2se3-stm-denchar.fdf
 
=== added file 'Tests/bi2se3-stm-denchar/bi2se3-stm-denchar.fdf'
--- Tests/bi2se3-stm-denchar/bi2se3-stm-denchar.fdf	1970-01-01 00:00:00 +0000
+++ Tests/bi2se3-stm-denchar/bi2se3-stm-denchar.fdf	2019-09-26 07:10:24 +0000
@@ -0,0 +1,185 @@
+#
+# Three quintuple layers of Bi2Se3 to test the calculation of the spin texture.
+# This calculation will produce HSX and WFSX files that can be processed by
+# the Util/COOP/spin_texture utility program.
+#
+# (thanks to Roberto Robles)
+#
+SystemName   bi2se3-stm-denchar
+SystemLabel  bi2se3-stm-denchar
+
+# Geom
+NumberOfAtoms       15
+NumberOfSpecies      2
+%block ChemicalSpeciesLabel
+   1  83 Bi
+   2  34 Se
+%endblock ChemicalSpeciesLabel
+LatticeConstant   4.13800000 Ang
+%block LatticeVectors
+     1.0000000000000000    0.0000000000000000    0.0000000000000000
+    -0.5000000000000000    0.8660254037844386    0.0000000000000000
+     0.0000000000000000    0.0000000000000000   11.7544706113098112
+%endblock LatticeVectors
+AtomicCoordinatesFormat Fractional
+%block AtomicCoordinatesAndAtomicSpecies
+  0.6666666666666667  0.3333333333333333  0.2650624783220366   1   1  Bi
+  0.3333333333333333  0.6666666666666667  0.3423936279930260   1   2  Bi
+  0.3333333333333333  0.6666666666666667  0.4613344252566094   1   3  Bi
+  0.0000000000000000  0.0000000000000000  0.5386655576677342   1   4  Bi
+  0.0000000000000000  0.0000000000000000  0.6576063590530268   1   5  Bi
+  0.6666666666666667  0.3333333333333333  0.7349374737996754   1   6  Bi
+  0.3333333333333333  0.6666666666666667  0.2287521714264216   2   7  Se
+  0.0000000000000000  0.0000000000000000  0.3037280662035826   2   8  Se
+  0.6666666666666667  0.3333333333333333  0.3787039395991627   2   9  Se
+  0.0000000000000000  0.0000000000000000  0.4250241177721818   2  10  Se
+  0.6666666666666667  0.3333333333333333  0.5000000000000000   2  11  Se
+  0.3333333333333333  0.6666666666666667  0.5749758645633420   2  12  Se
+  0.6666666666666667  0.3333333333333333  0.6212960521574118   2  13  Se
+  0.3333333333333333  0.6666666666666667  0.6962719343852299   2  14  Se
+  0.0000000000000000  0.0000000000000000  0.7712478166130481   2  15  Se
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+# Basis 
+PAO.EnergyShift  100 meV   ! Extension of AOs (smaller => more extended)
+PAO.SplitNorm      0.15
+PAO.SplitNormH      0.50
+PAO.SoftDefault      T
+%block PAO.Basis
+ Se   3
+   n=4   0   2   ! n,l,#zetas
+    0.0   0.0    ! default
+   n=4   1   2
+    0.0   0.0
+   n=3   2   1
+    0.0
+ Bi   3  
+   n=6   0   2
+    0.0   0.0
+   n=6   1   2
+    0.0   0.0
+   n=5   2   1
+    0.0
+%endblock PAO.Basis
+
+xc.functional         GGA          # Exchange-correlation functional
+xc.authors            PBE          # Exchange-correlation version
+
+%block DM.InitSpin
+ 1  0.  0. 0.
+ 2  0.  0. 0.
+ 3  0.  0. 0.
+ 4  0.  0. 0.
+ 5  0.  0. 0.
+ 6  0.  0. 0.
+ 7  0.  0. 0.
+ 8  0.  0. 0.
+ 9  0.  0. 0.
+10  0.  0. 0.
+11  0.  0. 0.
+12  0.  0. 0.
+13  0.  0. 0.
+14  0.  0. 0.
+15  0.  0. 0.
+%endblock DM.InitSpin
+
+SolutionMethod  diagon
+
+%block kgrid_Monkhorst_Pack
+   3   0   0   0.0
+   0   3   0   0.0
+   0   0   1   0.0
+%endblock kgrid_Monkhorst_Pack
+
+ ################ SCF options #############
+SolutionMethod  diagon
+ElectronicTemperature   0.0010 eV
+MaxSCFIterations   300
+DM.NumberPulay   3
+DM.MixingWeight     0.0100
+MeshCutoff        300.0000 Ry
+DM.Tolerance        1.E-4
+SCF.H.Converge  F
+UseSaveData     F
+DM.UseSaveDM    T
+SpinOrbit       T
+SaveElectrostaticPotential true
+
+# ---- SPIN TEXTURE
+#
+# We use a 'BandPoints' block to generate the states at chosen k-points
+# These keywords are necessary to generate the HSX and WFSX files
+#
+SaveHS T
+WFS.Write.For.Bands T      
+#
+# These keywords could be used to limit the number of bands stored in the WFSX file
+
+Wfs.band.min 190
+Wfs.band.max 230
+
+BandLinesScale      pi/a
+%block BandPoints
+  0.0 0.0 0.0 1
+  0.2 0.1 0.0 1.
+%endblock BandPoints
+
+
+
+#  *************************************************************************
+#                     Input variables for Denchar
+#  (besides SystemLabel, NumberOfSpecies and ChemicalSpecies, defined above)
+#  *************************************************************************
+
+Write.Denchar true
+Denchar.TypeOfRun      3D
+
+Denchar.PlotCharge          F
+Denchar.PlotWaveFunctions   T
+
+Denchar.CoorUnits       bohr           # Format for coordinate of the points
+                                       # Bohr
+                                       # Ang
+
+Denchar.DensityUnits   Ele/bohr**3     # Units of Charge Density
+                                       # Ele/bohr**3
+                                       # Ele/Ang**3
+                                       # Ele/UnitCell
+
+Denchar.MinX           -0.10  bohr      # Minimum coord of the window in X-dir
+Denchar.MaxX            4.0  bohr      # Maximum coord of the window in X-dir
+
+Denchar.MinY           -0.1  bohr      # Minimum coord of the window in Y-dir
+Denchar.MaxY            5.0  bohr      # Maximum coord of the window in Y-dir
+
+Denchar.MinZ            23.0  bohr      # Maximum coord of the window in X-dir
+Denchar.MaxZ            71.0  bohr      # Maximum coord of the window in X-dir
+
+Denchar.NumberPointsX    30            # Number of points in X-axis
+Denchar.NumberPointsY    30            # Number of points in Y-axis
+Denchar.NumberPointsZ    30            # Number of points in Y-axis
+
+
+#  *************************************************************************
+#                     Input variables for STM
+#  (besides SystemLabel, NumberOfSpecies and ChemicalSpecies, defined above)
+#  *************************************************************************
+STM.DensityUnits   Ele/bohr**3   # Units of Charge Density
+                                     # Ele/bohr**3
+                                     # Ele/Ang**3
+                                     # Ele/UnitCell
+
+STM.MinZ         37.0 Ang     # Minimum coordinate of the window in Z-dir
+STM.MaxZ         40.0 Ang     # Maximum coordinate of the window in Z-dir
+STM.RefZ         37.0 Ang     # Reference plane to propagate wf into vacuum
+STM.VacZ         40.0 Ang     # Reference plane to compute vacuum potential
+
+STM.Emin         -6.3 eV
+STM.Emax         -1.0 eV
+
+STM.NumberPointsX       48         # Number of points in X-axis
+STM.NumberPointsY       48         # Number of points in Y-axis
+STM.NumberPointsZ       48         # Number of points in Z-axis
+
+STM.NumberCellsX         1
+STM.NumberCellsY         1
 
=== added file 'Tests/bi2se3-stm-denchar/bi2se3-stm-denchar.pseudos'
--- Tests/bi2se3-stm-denchar/bi2se3-stm-denchar.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/bi2se3-stm-denchar/bi2se3-stm-denchar.pseudos	2019-09-26 07:10:24 +0000
@@ -0,0 +1,1 @@
+Bi Se
 
=== added directory 'Tests/bi2se3-sts-denchar'
=== added file 'Tests/bi2se3-sts-denchar/Makefile'
--- Tests/bi2se3-sts-denchar/Makefile	1970-01-01 00:00:00 +0000
+++ Tests/bi2se3-sts-denchar/Makefile	2019-09-26 07:10:24 +0000
@@ -0,0 +1,6 @@
+#
+# Single-test makefile
+#
+name=bi2se3-sts-denchar
+#
+include ../test.mk
 
=== added file 'Tests/bi2se3-sts-denchar/README'
--- Tests/bi2se3-sts-denchar/README	1970-01-01 00:00:00 +0000
+++ Tests/bi2se3-sts-denchar/README	2019-09-26 07:10:24 +0000
@@ -0,0 +1,9 @@
+This example also involves denchar and stm.
+
+To complete this example you also need to
+run the following things:
+
+1. Go into work directory
+2. Create a soft-link using: ln -s bi2se3-sts-denchar.bands.WFSX bi2se3-sts-denchar.WFSX
+3. Run denchar < ../bi2se3-sts-denchar.fdf
+4. Run stm < ../bi2se3-sts-denchar.fdf
 
=== added file 'Tests/bi2se3-sts-denchar/bi2se3-sts-denchar.fdf'
--- Tests/bi2se3-sts-denchar/bi2se3-sts-denchar.fdf	1970-01-01 00:00:00 +0000
+++ Tests/bi2se3-sts-denchar/bi2se3-sts-denchar.fdf	2019-09-26 07:10:24 +0000
@@ -0,0 +1,182 @@
+#
+# Three quintuple layers of Bi2Se3 to test the calculation of the spin texture.
+# This calculation will produce HSX and WFSX files that can be processed by
+# the Util/COOP/spin_texture utility program.
+#
+# (thanks to Roberto Robles)
+#
+SystemName   bi2se3-sts-denchar
+SystemLabel  bi2se3-sts-denchar
+
+# Geom
+NumberOfAtoms       15
+NumberOfSpecies      2
+%block ChemicalSpeciesLabel
+   1  83 Bi
+   2  34 Se
+%endblock ChemicalSpeciesLabel
+LatticeConstant   4.13800000 Ang
+%block LatticeVectors
+     1.0000000000000000    0.0000000000000000    0.0000000000000000
+    -0.5000000000000000    0.8660254037844386    0.0000000000000000
+     0.0000000000000000    0.0000000000000000   11.7544706113098112
+%endblock LatticeVectors
+AtomicCoordinatesFormat Fractional
+%block AtomicCoordinatesAndAtomicSpecies
+  0.6666666666666667  0.3333333333333333  0.2650624783220366   1   1  Bi
+  0.3333333333333333  0.6666666666666667  0.3423936279930260   1   2  Bi
+  0.3333333333333333  0.6666666666666667  0.4613344252566094   1   3  Bi
+  0.0000000000000000  0.0000000000000000  0.5386655576677342   1   4  Bi
+  0.0000000000000000  0.0000000000000000  0.6576063590530268   1   5  Bi
+  0.6666666666666667  0.3333333333333333  0.7349374737996754   1   6  Bi
+  0.3333333333333333  0.6666666666666667  0.2287521714264216   2   7  Se
+  0.0000000000000000  0.0000000000000000  0.3037280662035826   2   8  Se
+  0.6666666666666667  0.3333333333333333  0.3787039395991627   2   9  Se
+  0.0000000000000000  0.0000000000000000  0.4250241177721818   2  10  Se
+  0.6666666666666667  0.3333333333333333  0.5000000000000000   2  11  Se
+  0.3333333333333333  0.6666666666666667  0.5749758645633420   2  12  Se
+  0.6666666666666667  0.3333333333333333  0.6212960521574118   2  13  Se
+  0.3333333333333333  0.6666666666666667  0.6962719343852299   2  14  Se
+  0.0000000000000000  0.0000000000000000  0.7712478166130481   2  15  Se
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+# Basis 
+PAO.EnergyShift  100 meV   ! Extension of AOs (smaller => more extended)
+PAO.SplitNorm      0.15
+PAO.SplitNormH      0.50
+PAO.SoftDefault      T
+%block PAO.Basis
+ Se   3
+   n=4   0   2   ! n,l,#zetas
+    0.0   0.0    ! default
+   n=4   1   2
+    0.0   0.0
+   n=3   2   1
+    0.0
+ Bi   3  
+   n=6   0   2
+    0.0   0.0
+   n=6   1   2
+    0.0   0.0
+   n=5   2   1
+    0.0
+%endblock PAO.Basis
+
+xc.functional         GGA          # Exchange-correlation functional
+xc.authors            PBE          # Exchange-correlation version
+
+%block DM.InitSpin
+ 1  0.  0. 0.
+ 2  0.  0. 0.
+ 3  0.  0. 0.
+ 4  0.  0. 0.
+ 5  0.  0. 0.
+ 6  0.  0. 0.
+ 7  0.  0. 0.
+ 8  0.  0. 0.
+ 9  0.  0. 0.
+10  0.  0. 0.
+11  0.  0. 0.
+12  0.  0. 0.
+13  0.  0. 0.
+14  0.  0. 0.
+15  0.  0. 0.
+%endblock DM.InitSpin
+
+SolutionMethod  diagon
+
+%block kgrid_Monkhorst_Pack
+   3   0   0   0.0
+   0   3   0   0.0
+   0   0   1   0.0
+%endblock kgrid_Monkhorst_Pack
+
+ ################ SCF options #############
+SolutionMethod  diagon
+ElectronicTemperature   0.0010 eV
+MaxSCFIterations   300
+DM.NumberPulay   3
+DM.MixingWeight     0.0100
+MeshCutoff        300.0000 Ry
+DM.Tolerance        1.E-4
+SCF.H.Converge  F
+UseSaveData     F
+DM.UseSaveDM    T
+SpinOrbit       T
+SaveElectrostaticPotential true
+
+# ---- SPIN TEXTURE
+#
+# We use a 'BandPoints' block to generate the states at chosen k-points
+# These keywords are necessary to generate the HSX and WFSX files
+#
+SaveHS T
+WFS.Write.For.Bands T      
+#
+# These keywords could be used to limit the number of bands stored in the WFSX file
+
+Wfs.band.min 190
+Wfs.band.max 230
+
+BandLinesScale      pi/a
+%block BandPoints
+  0.0 0.0 0.0 1
+  0.2 0.1 0.0 1.
+%endblock BandPoints
+
+
+#  *************************************************************************
+#                     Input variables for Denchar
+#  (besides SystemLabel, NumberOfSpecies and ChemicalSpecies, defined above)
+#  *************************************************************************
+
+Write.Denchar true
+Denchar.TypeOfRun      3D
+
+Denchar.PlotCharge          F
+Denchar.PlotWaveFunctions   T
+
+Denchar.CoorUnits       bohr           # Format for coordinate of the points
+                                       # Bohr
+                                       # Ang
+
+Denchar.DensityUnits   Ele/bohr**3     # Units of Charge Density
+                                       # Ele/bohr**3
+                                       # Ele/Ang**3
+                                       # Ele/UnitCell
+
+Denchar.MinX           -0.10  bohr      # Minimum coord of the window in X-dir
+Denchar.MaxX            4.0  bohr      # Maximum coord of the window in X-dir
+
+Denchar.MinY           -0.1  bohr      # Minimum coord of the window in Y-dir
+Denchar.MaxY            5.0  bohr      # Maximum coord of the window in Y-dir
+
+Denchar.MinZ            23.0  bohr      # Maximum coord of the window in X-dir
+Denchar.MaxZ            71.0  bohr      # Maximum coord of the window in X-dir
+
+Denchar.NumberPointsX    30            # Number of points in X-axis
+Denchar.NumberPointsY    30            # Number of points in Y-axis
+Denchar.NumberPointsZ    30            # Number of points in Y-axis
+
+
+#  *************************************************************************
+#                     Input variables for STM
+#  (besides SystemLabel, NumberOfSpecies and ChemicalSpecies, defined above)
+#  *************************************************************************
+STS.NumberOfPoints 30
+STM.DensityUnits   Ele/bohr**3   # Units of Charge Density
+                                     # Ele/bohr**3
+                                     # Ele/Ang**3
+                                     # Ele/UnitCell
+
+STM.MinZ         37.0 Ang     # Minimum coordinate of the window in Z-dir
+STM.MaxZ         37.0 Ang     # Maximum coordinate of the window in Z-dir
+STM.RefZ         37.0 Ang     # Reference plane to propagate wf into vacuum
+STM.VacZ         40.0 Ang     # Reference plane to compute vacuum potential
+
+STM.Emin         -6.3 eV
+STM.Emax         -1.0 eV
+
+STM.NumberPointsX       48         # Number of points in X-axis
+STM.NumberPointsY       48         # Number of points in Y-axis
+STM.NumberPointsZ        1         # Number of points in Z-axis
 
=== added file 'Tests/bi2se3-sts-denchar/bi2se3-sts-denchar.pseudos'
--- Tests/bi2se3-sts-denchar/bi2se3-sts-denchar.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/bi2se3-sts-denchar/bi2se3-sts-denchar.pseudos	2019-09-26 07:10:24 +0000
@@ -0,0 +1,1 @@
+Bi Se
 
=== modified file 'Util/COOP/dm_creator.F90'
--- Util/COOP/dm_creator.F90	2019-02-26 09:08:49 +0000
+++ Util/COOP/dm_creator.F90	2019-09-26 07:10:24 +0000
@@ -78,13 +78,18 @@
   allocate (wk(nkp), pk(3,nkp))
 
   read(wfs_u) nsp
-  non_coll = (nsp == 4)
   read(wfs_u) nao
   read(wfs_u)        !! Symbols, etc
   if (debug) print *, "WFSX read: nkp, nsp, nnao: ", nkp, nsp, nao
 
-  if (non_coll) then
-     nspin_blocks = 1
+  if ( nsp == 8 ) then
+    if (debug) print *, "WFSX from spin-orbit calculation"
+    nspin_blocks = 1
+    non_coll = .true.
+  else if ( nsp == 4 ) then
+    if (debug) print *, "WFSX from non-collinear calculation"
+    nspin_blocks = 1
+    non_coll = .true.
   else
      nspin_blocks = nsp
   endif
 
=== modified file 'Util/COOP/fat.f90'
--- Util/COOP/fat.f90	2019-02-15 15:50:41 +0000
+++ Util/COOP/fat.f90	2019-09-26 07:10:24 +0000
@@ -120,7 +120,7 @@
   allocate (wk(nkp), pk(3,nkp))
 
   read(wfs_u) nsp
-  non_coll = (nsp == 4)
+  non_coll = (nsp >= 4)
   read(wfs_u) nao
   read(wfs_u)        !! Symbols, etc
   if (debug) print *, "WFSX read: nkp, nsp, nnao: ", nkp, nsp, nao
@@ -273,10 +273,12 @@
      enddo
   enddo
 
-  if (non_coll) then
-     write(stt_u,"(/'SPIN (non-coll): ',i2)") nspin_blocks
+  if ( nsp == 8 ) then
+    write(stt_u,"(/'SPIN (spin-orbit): ',i2)") nspin_blocks
+  else if ( nsp == 4 ) then
+    write(stt_u,"(/'SPIN (non-coll): ',i2)") nspin_blocks
   else
-     write(stt_u,"(/'SPIN: ',i2)") nspin_blocks
+    write(stt_u,"(/'SPIN: ',i2)") nspin_blocks
   endif
   
   write(stt_u,"(/'AO LIST:')")
 
=== modified file 'Util/COOP/mprop.f90'
--- Util/COOP/mprop.f90	2019-02-27 15:31:08 +0000
+++ Util/COOP/mprop.f90	2019-09-26 07:10:24 +0000
@@ -141,7 +141,7 @@
   allocate (wk(nkp), pk(3,nkp))
 
   read(wfs_u) wfs_spin_flag   !  1, 2, or 4 
-  non_coll = (wfs_spin_flag == 4)
+  non_coll = (wfs_spin_flag >= 4)
   read(wfs_u) nao
   read(wfs_u)        !! Symbols, etc
   if (debug) print *, "WFSX read: nkp, nspin_flag, nnao: ", nkp, wfs_spin_flag, nao
 
=== modified file 'Util/COOP/spin_texture.f90'
--- Util/COOP/spin_texture.f90	2019-06-07 16:27:46 +0000
+++ Util/COOP/spin_texture.f90	2019-09-26 07:10:24 +0000
@@ -102,7 +102,7 @@
   allocate (wk(nkp), pk(3,nkp))
 
   read(wfs_u) nsp
-  non_coll = (nsp == 4)
+  non_coll = (nsp >= 4)
   if (.not. (non_coll)) then
      STOP "Spin texture not available for collinear-spin"
   endif
 
=== modified file 'Util/Denchar/Src/denchar.f'
--- Util/Denchar/Src/denchar.f	2019-06-08 15:41:39 +0000
+++ Util/Denchar/Src/denchar.f	2019-09-26 07:10:24 +0000
@@ -264,7 +264,7 @@
 
       ! Non-collinear or SOC files have a single "spin" block as opposed to collinear-spin
       ! files, which contain two spin blocks per k section.
-      non_coll = (nspin_wfsx == 4)
+      non_coll = (nspin_wfsx >= 4)
       nspin_blocks = nspin_wfsx
       if (non_coll) nspin_blocks = 1
 
 
=== modified file 'Util/Denchar/Src/wavofr.f'
--- Util/Denchar/Src/wavofr.f	2019-06-08 15:41:39 +0000
+++ Util/Denchar/Src/wavofr.f	2019-09-26 07:10:24 +0000
@@ -49,7 +49,7 @@
 C INTEGER NO               : Total number of orbitals in Supercell
 C INTEGER NO_U             : Total number of orbitals in Unit Cell
 C INTEGER MAXNA            : Maximum number of neighbours of any atom
-C INTEGER NSPIN            : Number of different spin polarizations: 1, 2, 4
+C INTEGER NSPIN            : Number of different spin polarizations: 1, 2, 4, 8
 ! integer nspin_blocks     : blocks in WFSX file       
 ! logical non_coll         : NC/SOC wfs data?
 C INTEGER ISA(NA)          : Species index of each atom
@@ -174,6 +174,17 @@
 
       PI = 4.0D0 * ATAN(1.0D0)
 
+      select case ( NSPIN )
+      case ( 1, 2, 4, 8 )
+! ok
+      case default
+        WRITE(6,*)'BAD NUMBER NSPIN IN WAVOFR.F'
+        WRITE(6,*)'NSPIN = ',NSPIN
+        WRITE(6,*)'IT MUST BE 1, 2, 4 or 8'
+        STOP
+      end select
+
+
       NPLAMAX = NPX * NPY * NPZ
 
       ALLOCATE(PLAPO(NPLAMAX,3))
@@ -343,7 +354,7 @@
      .         FORM = 'FORMATTED')
           REWIND(UNITPH1)
 
-       ELSEIF ((NSPIN .EQ. 2) .or. (NSPIN == 4))  THEN
+       ELSE ! nspin >= 2
           ! We will reuse 'up' and 'down' for the spinor components
            IF (IDIMEN .EQ. 2) THEN
             FN_URE = TRIM(SNAME)//'.CON.K' //
@@ -405,11 +416,6 @@
           OPEN(UNIT = UNITPH2, FILE = FN_DPH, STATUS = 'UNKNOWN',
      .         FORM = 'FORMATTED')
           REWIND(UNITPH2)
-        ELSE
-          WRITE(6,*)'BAD NUMBER NSPIN IN WAVOFR.F'
-          WRITE(6,*)'NSPIN = ',NSPIN
-          WRITE(6,*)'IT MUST BE 1, 2, or 4'
-          STOP
         ENDIF
 
         IF (IDIMEN .EQ. 2) THEN
@@ -465,16 +471,16 @@
                 call write_cube_header(unitmo2,fn_dmo)
                 call write_cube_header(unitph2,fn_dph)
              endif
-          ELSE IF (NSPIN .EQ. 4) THEN
-                call write_cube_header(unitre1,fn_ure)
-                call write_cube_header(unitim1,fn_uim)
-                call write_cube_header(unitmo1,fn_umo)
-                call write_cube_header(unitph1,fn_uph)
-
-                call write_cube_header(unitre2,fn_dre)
-                call write_cube_header(unitim2,fn_dim)
-                call write_cube_header(unitmo2,fn_dmo)
-                call write_cube_header(unitph2,fn_dph)
+          ELSE ! 4 or 8
+            call write_cube_header(unitre1,fn_ure)
+            call write_cube_header(unitim1,fn_uim)
+            call write_cube_header(unitmo1,fn_umo)
+            call write_cube_header(unitph1,fn_uph)
+            
+            call write_cube_header(unitre2,fn_dre)
+            call write_cube_header(unitim2,fn_dim)
+            call write_cube_header(unitmo2,fn_dmo)
+            call write_cube_header(unitph2,fn_dph)
 
           ENDIF
         ENDIF
@@ -576,7 +582,7 @@
                  call mod_and_phase(rwavedn,iwavedn,mwavedn,pwavedn)
               endif
 
-           ELSEIF (NSPIN .EQ. 4) THEN
+           ELSE ! 4 or 8
 
               rwaveup = real(cwave(1), dp)
               iwaveup = aimag(cwave(1))
@@ -605,19 +611,19 @@
                  call write_plapo(unitmo2,mwavedn)
                  call write_plapo(unitph2,pwavedn)
                 endif
-             ELSEIF ( NSPIN .EQ. 4 ) THEN
-
-                 call write_plapo(unitre1,rwaveup)
-                 call write_plapo(unitim1,iwaveup)
-                 call write_plapo(unitmo1,mwaveup)
-                 call write_plapo(unitph1,pwaveup)
-
-                 call write_plapo(unitre2,rwavedn)
-                 call write_plapo(unitim2,iwavedn)
-                 call write_plapo(unitmo2,mwavedn)
-                 call write_plapo(unitph2,pwavedn)
-
-              ENDIF
+             ELSE ! 4 or 8
+
+               call write_plapo(unitre1,rwaveup)
+               call write_plapo(unitim1,iwaveup)
+               call write_plapo(unitmo1,mwaveup)
+               call write_plapo(unitph1,pwaveup)
+               
+               call write_plapo(unitre2,rwavedn)
+               call write_plapo(unitim2,iwavedn)
+               call write_plapo(unitmo2,mwavedn)
+               call write_plapo(unitph2,pwavedn)
+               
+             ENDIF
               
           ELSE IF (IDIMEN .EQ. 3) THEN
              IF (NSPIN .EQ. 1) THEN
@@ -637,17 +643,17 @@
                   MWF(NPO,2) = MWAVEDN
                   PWF(NPO,2) = PWAVEDN
                endif
-             ELSE IF (NSPIN .EQ. 4) THEN
-
-                RWF(NPO,1) = RWAVEUP
-                IMWF(NPO,1) = IWAVEUP
-                MWF(NPO,1) = MWAVEUP
-                PWF(NPO,1) = PWAVEUP
-
-                RWF(NPO,2) = RWAVEDN
-                IMWF(NPO,2) = IWAVEDN
-                MWF(NPO,2) = MWAVEDN
-                PWF(NPO,2) = PWAVEDN
+             ELSE ! 4 or 8
+
+               RWF(NPO,1) = RWAVEUP
+               IMWF(NPO,1) = IWAVEUP
+               MWF(NPO,1) = MWAVEUP
+               PWF(NPO,1) = PWAVEUP
+
+               RWF(NPO,2) = RWAVEDN
+               IMWF(NPO,2) = IWAVEDN
+               MWF(NPO,2) = MWAVEDN
+               PWF(NPO,2) = PWAVEDN
 
             ENDIF
          ENDIF
@@ -668,32 +674,30 @@
            DO NX=1,NPX
               DO NY=1,NPY
 
-                 if ( (NSPIN .EQ. 1)  .or.
-     $                ((NSPIN .EQ.2).and.(ispin.eq.1)) .or. 
-     $                (NSPIN .EQ.4)  ) then
-
-                    WRITE(UNITRE1,'(6e13.5)')
-     $                   (RWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ)
-                    WRITE(UNITIM1,'(6e13.5)')
-     $                  (IMWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ)
-                    WRITE(UNITMO1,'(6e13.5)')
-     $                   (MWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ)
-                    WRITE(UNITPH1,'(6e13.5)')
-     $                   (PWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ)
-                 endif
+                if ( NSPIN == 1 .or.
+     &              (NSPIN == 2 .and. ispin == 1 ) .or.
+     &              NSPIN >= 4 ) then
+                  WRITE(UNITRE1,'(6e13.5)')
+     $                (RWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ)
+                  WRITE(UNITIM1,'(6e13.5)')
+     $                (IMWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ)
+                  WRITE(UNITMO1,'(6e13.5)')
+     $                (MWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ)
+                  WRITE(UNITPH1,'(6e13.5)')
+     $                (PWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ)
+                end if
                  
-                 if ( ((NSPIN .EQ.2).and.(ispin.eq.2)) .or. 
-     $                   (NSPIN .EQ.4)  ) then
-
-                    WRITE(UNITRE2,'(6e13.5)')
-     $                   (RWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ)
-                    WRITE(UNITIM2,'(6e13.5)')
-     $                  (IMWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ)
-                    WRITE(UNITMO2,'(6e13.5)')
-     $                   (MWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ)
-                    WRITE(UNITPH2,'(6e13.5)')
-     $                   (PWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ)
-                 endif
+                if ( (NSPIN == 2 .and. ispin == 2) .or.
+     &              NSPIN >= 4 ) then
+                  WRITE(UNITRE2,'(6e13.5)')
+     $                (RWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ)
+                  WRITE(UNITIM2,'(6e13.5)')
+     $                (IMWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ)
+                  WRITE(UNITMO2,'(6e13.5)')
+     $                (MWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ)
+                  WRITE(UNITPH2,'(6e13.5)')
+     $                (PWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ)
+                end if
               ENDDO
            ENDDO
 
@@ -705,7 +709,7 @@
             WRITE(6,'(A,A)') '   ',FN_IM
             WRITE(6,'(A,A)') '   ',FN_MO
             WRITE(6,'(A,A)') '   ',FN_PH
-          ELSE IF ((NSPIN .EQ. 2) .or. (NSPIN .EQ. 4)) THEN
+          ELSE ! 2, 4 or 8
             WRITE(6,'(A,A)') '   ',FN_URE
             WRITE(6,'(A,A)') '   ',FN_UIM
             WRITE(6,'(A,A)') '   ',FN_DRE
@@ -714,7 +718,7 @@
             WRITE(6,'(A,A)') '   ',FN_UPH
             WRITE(6,'(A,A)') '   ',FN_DMO
             WRITE(6,'(A,A)') '   ',FN_DPH
-            if (nspin == 4) write(6,"(a)")
+            if (nspin >= 4) write(6,"(a)")
      $           "... up and down for spinor components"
           ENDIF
         ENDIF ! 3D
 
=== modified file 'Util/Eig2DOS/Eig2DOS.f90'
--- Util/Eig2DOS/Eig2DOS.f90	2019-03-02 23:57:58 +0000
+++ Util/Eig2DOS/Eig2DOS.f90	2019-09-26 07:10:24 +0000
@@ -56,12 +56,12 @@
   real(dp) :: emin        = -huge(1.0_dp)
   real(dp) :: emax        =  huge(1.0_dp)
   real(dp) :: smear       = 0.2_dp
+  real(dp) :: smear2
   logical  :: loren       = .false.
   logical  :: emin_given  = .false.
   logical  :: emax_given  = .false.
   logical  :: using_weights = .false.
   logical  :: shift_efermi = .false.
-  logical  :: non_coll
   integer  :: min_band = 0
   integer  :: max_band = 0
 
@@ -155,13 +155,18 @@
        action="read")
   read(1,*) Ef
   read(1,*) nband, nspin, nk
-  non_coll = (nspin == 4)
 
-  if (non_coll) then
-     nspin_blocks = 1
-  else
-     nspin_blocks = nspin
-  endif
+  nspin_blocks = 1
+  if ( nspin == 8 ) then
+    write(*,"(a)") "# Eigenvalues calculated from a spin-orbit calculation"
+  else if ( nspin == 4 ) then
+    write(*,"(a)") "# Eigenvalues calculated from a non-collinear calculation"
+  else if ( nspin == 2 ) then
+    write(*,"(a)") "# Eigenvalues calculated from a spin-polarized calculation"
+    nspin_blocks = nspin
+  else if ( nspin == 1 ) then
+    write(*,"(a)") "# Eigenvalues calculated from a non-polarized calculation"
+  end if
   
   write(*,"(a)") "# Eigenvalues read from " // trim(eig_file)
   if ( debug ) print *, "Ef, nband, nspin, nk:", Ef, nband, nspin, nk
@@ -197,6 +202,7 @@
   endif
 
 
+  smear2 = smear * smear
   write(*,"(a,f8.4)") "# Using smearing parameter: ", smear
   write(*,"(a,i0,a)") "# Using ", npts_energy, " points in the energy range"
   write(*,"(2(a,i0))") "# Selected bands: ", min_band ," to: ", max_band
@@ -253,9 +259,9 @@
            do ie = 1, npts_energy
               x = emin + (ie-1)*eincr - e
               if (loren) then
-                 DOS(ie,is) = DOS(ie,is) + weight * smear / (smear*smear + x*x) 
+                 DOS(ie,is) = DOS(ie,is) + weight * smear / (smear2 + x*x) 
               else
-                 DOS(ie,is) = DOS(ie,is) + weight * exp( - x*x/(smear*smear) )
+                 DOS(ie,is) = DOS(ie,is) + weight * exp( - x*x/smear2 )
               end if
            end do
         end do
 
=== modified file 'Util/STM/ol-stm/Src/mainstm.F'
--- Util/STM/ol-stm/Src/mainstm.F	2019-06-06 09:34:48 +0000
+++ Util/STM/ol-stm/Src/mainstm.F	2019-09-26 07:10:24 +0000
@@ -158,7 +158,7 @@
 
       ! Non-collinear or SOC files have a single "spin" block as opposed to collinear-spin
       ! files, which contain two spin blocks per k section.
-      non_coll = (nspin == 4)
+      non_coll = (nspin >= 4)
       nspin_blocks = nspin
       if (non_coll) nspin_blocks = 1
 
 
=== modified file 'Util/STM/ol-stm/Src/stm.f'
--- Util/STM/ol-stm/Src/stm.f	2019-06-05 08:47:53 +0000
+++ Util/STM/ol-stm/Src/stm.f	2019-09-26 07:10:24 +0000
@@ -104,6 +104,7 @@
  
       LOGICAL FIRST
       integer :: idummy, number_of_wfns, spinor_comps
+      integer :: nspin_max
 
       CHARACTER ::  SNAME*40, FNAME*256, stm_label*60
 
@@ -168,10 +169,12 @@
       ALLOCATE(R2IJ(MAXNA))
       ALLOCATE(XIJ(3,MAXNA))
 
+      nspin_max = min(4,nspin)
+
       allocate(CWAVE(spinor_comps))
       ALLOCATE(CW(0:NPX-1,0:NPY-1,spinor_comps))
       ALLOCATE(CWE(0:NPX-1,0:NPY-1,0:NPZ-1,spinor_comps))
-      ALLOCATE(RHO(0:NPX-1,0:NPY-1,0:NPZ-1,nspin))
+      ALLOCATE(RHO(0:NPX-1,0:NPY-1,0:NPZ-1,nspin_max))
 
       FIRST = .TRUE.
       DO I = 1,3
@@ -423,9 +426,9 @@
      $              [0.0_dp, 0.0_dp, ZMIN], ! Extra info for origin
      $              [.true.,.true.,.false.] ! Periodic ?
 
-      WRITE(grid_u) NPX, NPY, NPZ, nspin 
+      WRITE(grid_u) NPX, NPY, NPZ, nspin_max 
 
-      do ispin = 1, nspin
+      do ispin = 1, nspin_max
          DO IZ=0,NPZ-1
             DO IY=0,NPY-1
                WRITE(grid_u) (REAL(RHO(IX,IY,IZ,ispin),sp),IX=0,NPX-1)
 
=== modified file 'Util/WFS/readwfx.f90'
--- Util/WFS/readwfx.f90	2019-03-18 10:51:45 +0000
+++ Util/WFS/readwfx.f90	2019-09-26 07:10:24 +0000
@@ -121,7 +121,7 @@
   read(wfs_u) nk, gamma
 
   read(wfs_u) wfs_spin_flag   !  1, 2, or 4 
-  non_coll = (wfs_spin_flag == 4)
+  non_coll = (wfs_spin_flag >= 4)
   read(wfs_u) nuotot
   read(wfs_u)        !! Symbols, etc
 
 
=== modified file 'Util/WFS/wfs2wfsx.f'
--- Util/WFS/wfs2wfsx.f	2016-01-25 16:00:16 +0000
+++ Util/WFS/wfs2wfsx.f	2019-09-26 07:10:24 +0000
@@ -75,7 +75,7 @@
         orb_info_written = .false.
 
         do iik = 1,nk
-          do iispin = 1,nspin
+          do iispin = 1,min(4,nspin)
 
           read(iu) ik,k(1),k(2),k(3)
           if (ik .ne. iik) stop 'error in index of k-point'
 
=== modified file 'Util/WFS/wfsnc2wfsx.F90'
--- Util/WFS/wfsnc2wfsx.F90	2016-01-25 16:00:16 +0000
+++ Util/WFS/wfsnc2wfsx.F90	2019-09-26 07:10:24 +0000
@@ -103,7 +103,7 @@
                  symfio(j), j=1,nuotot)
 
         do iik = 1,nk
-          do iispin = 1,nspin
+          do iispin = 1,min(4,nspin)
 
              write(io) iik,kp(1:3,iik), wk(iik)
              write(io) iispin
 
=== modified file 'Util/WFS/wfsx2wfs.f'
--- Util/WFS/wfsx2wfs.f	2016-09-12 13:30:18 +0000
+++ Util/WFS/wfsx2wfs.f	2019-09-26 07:10:24 +0000
@@ -61,7 +61,7 @@
      .            symfio(j), j=1,nuotot)
 
         do iik = 1,nk
-          do iispin = 1,nspin
+          do iispin = 1,min(4,nspin)
 
           read(iu) ik,k(1),k(2),k(3)
           write(io) ik,k(1),k(2),k(3)
 
=== modified file 'version.info'
--- version.info	2019-09-25 09:49:02 +0000
+++ version.info	2019-09-26 07:10:24 +0000
@@ -1,1 +1,5 @@
+<<<<<<< TREE
 siesta-4.1-1126
+=======
+siesta-4.1-1119--streamline-file-dim-nc-soc-2
+>>>>>>> MERGE-SOURCE