.gitignore

@@ -10,3 +10,5 @@ build.*
 *.x
 *.swp
 tags
+.DS_Store
+.vscode

.gitlab-ci.yml

@@ -8,9 +8,10 @@ stages:
 build-gfortran-hdf5:
   image: iffregistry.fz-juelich.de/fleur/fleur:latest
   stage: build
-  cache: 
+  artifacts:
      paths:
        - build
+     expire_in: 1h
   script:
     - cd /builds/fleur/fleur; ./configure.sh GITLAB; cd build; make -j 4
 #  only:
@@ -21,9 +22,6 @@ build-gfortran-hdf5:
 test-gfortran-hdf5:
   image: iffregistry.fz-juelich.de/fleur/fleur:latest
   stage: test
-  cache:
-     paths:
-       - build
   script:
     - ulimit -s unlimited ;export juDFT_MPI="mpirun -n 2 --allow-run-as-root ";export OMP_NUM_THREADS=4;cd /builds/fleur/fleur/build;ctest
   artifacts:
@@ -39,10 +37,6 @@ test-gfortran-hdf5:
 pages:
   image: iffregistry.fz-juelich.de/fleur/fleur:latest
   stage: deploy
-  cache:
-     paths:
-        - build
-        - public
   script:
     - echo "HTML should be ready from cache..."
     - mv /builds/fleur/fleur/docs/Docu_main.html /builds/fleur/fleur/public/index.html
@@ -60,10 +54,6 @@ pages:
 doxygen:
   image: iffregistry.fz-juelich.de/fleur/fleur:latest
   stage: html
-  cache:
-     paths:
-       - build
-       - public
   script:
     - cd /builds/fleur/fleur/build ; make doc
     - mkdir ../public
@@ -83,9 +73,10 @@ doxygen:
 build-pgi:
   image: iffregistry.fz-juelich.de/fleur/fleur:pgi
   stage: build
-  cache: 
+  artifacts: 
      paths:
        - build.pgi
+     expire_in: 1h
   script:
     - cd /builds/fleur/fleur; ./configure.sh -l pgi ; cd build.pgi; make 
   allow_failure: true
@@ -100,9 +91,6 @@ test-pgi:
   stage: test
   dependencies:
     - build-pgi
-  cache:
-     paths:
-       - build.pgi
   script:
     - cd /builds/fleur/fleur/build.pgi;ctest
   allow_failure: true
@@ -115,9 +103,10 @@ test-pgi:
 build-intel-static:
   image: iffregistry.fz-juelich.de/fleur/fleur:intel-static
   stage: build
-  cache: 
+  artifacts: 
      paths:
        - build.intel-static
+     expire_in: 1h
   script:
     - set +e && source compilervars.sh intel64 && set -e ; ulimit -s unlimited
     - cd /builds/fleur/fleur; ./configure.sh -l intel-static INTEL_DOCKER_STATIC ; cd build.intel-static; make -j 4
@@ -137,9 +126,10 @@ build-intel-static:
 build-intel:
   image: iffregistry.fz-juelich.de/fleur/fleur:intel-static
   stage: build
-  cache: 
+  artifacts: 
      paths:
        - build.intel.debug
+     expire_in: 1h
   script:
     - set +e && source compilervars.sh intel64 && set -e ; ulimit -s unlimited
     - cd /builds/fleur/fleur; CC=gcc FC=mpiifort  FLEUR_LIBRARIES="-lmkl_scalapack_lp64;-lmkl_blacs_intelmpi_lp64" ./configure.sh -t -d -l intel INTEL_MPI ; cd build.intel.debug; make -j 4
@@ -155,9 +145,6 @@ test-intel:
   stage: test
   dependencies:
      - build-intel
-  cache:
-     paths:
-       - build.intel.debug
   script:
     - set +e && source compilervars.sh intel64 && set -e; ulimit -s unlimited
     - cd /builds/fleur/fleur/build.intel.debug;ctest
@@ -175,9 +162,6 @@ test-intel:
 gfortran-coverage:
   image: iffregistry.fz-juelich.de/fleur/fleur:latest
   stage: html
-  cache: 
-     paths:
-        - build
   script:
       - cd /builds/fleur/fleur; ./configure.sh -l coverage -flags --coverage GITLAB; cd build.coverage; make -j 4
       - lcov --capture --initial -d CMakeFiles -o baseline.info

cdn/cdntot.f90

@@ -74,6 +74,54 @@ CONTAINS
       END DO ! loop over spins
    END SUBROUTINE integrate_cdn
 
+   SUBROUTINE integrate_realspace(xcpot, atoms, sym, sphhar, input, &
+                                  stars, cell, oneD, vacuum, noco, mt, is, hint)
+      use m_types
+      use m_mt_tofrom_grid
+      use m_pw_tofrom_grid
+      use m_constants
+      implicit none
+      CLASS(t_xcpot), INTENT(inout)   :: xcpot
+      TYPE(t_atoms),INTENT(IN)      :: atoms
+      TYPE(t_sym), INTENT(in)       :: sym
+      TYPE(t_sphhar), INTENT(IN)    :: sphhar
+      TYPE(t_input), INTENT(IN)     :: input
+      TYPE(t_stars), INTENT(IN)     :: stars
+      TYPE(t_cell), INTENT(IN)      :: cell
+      TYPE(t_oneD), INTENT(in)      :: oneD
+      TYPE(t_vacuum), INTENT(in)    :: vacuum
+      TYPE(t_noco), INTENT(in)      :: noco
+      real, intent(inout)           :: mt(:,:,:), is(:,:)
+      character(len=*), intent(in), optional :: hint
+      integer                       :: n_atm, i
+
+      TYPE(t_potden)                :: tmp_potden
+      REAL                          :: q(input%jspins), qis(input%jspins), &
+                                       qmt(atoms%ntype,input%jspins), qvac(2,input%jspins),&
+                                       qtot, qistot
+
+      call tmp_potden%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN)
+      call init_mt_grid(input%jspins, atoms, sphhar, xcpot, sym)
+      do n_atm =1,atoms%ntype
+         call mt_from_grid(atoms, sphhar, n_atm, input%jspins, mt(:,:,n_atm), &
+                           tmp_potden%mt(:,0:,n_atm,:))
+
+         do i=1,atoms%jri(n_atm)
+            tmp_potden%mt(i,:,n_atm,:) = tmp_potden%mt(i,:,n_atm,:) * atoms%rmsh(i,n_atm)**2
+         enddo
+      enddo
+      call finish_mt_grid()
+
+      call init_pw_grid(xcpot, stars, sym, cell)
+      call pw_from_grid(xcpot, stars, .False., is, tmp_potden%pw)
+      call finish_pw_grid()
+
+      call integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, tmp_potden, &
+                                   q, qis, qmt, qvac, qtot, qistot)
+
+      call print_cdn_inte(q, qis, qmt, qvac, qtot, qistot, hint)
+   END SUBROUTINE integrate_realspace
+
    SUBROUTINE cdntot(stars,atoms,sym,vacuum,input,cell,oneD,&
                      den,l_printData,qtot,qistot)
 
@@ -154,4 +202,29 @@ CONTAINS
 
       CALL timestop("cdntot")
    END SUBROUTINE cdntot
+
+   SUBROUTINE print_cdn_inte(q, qis, qmt, qvac, qtot, qistot, hint)
+      use  ieee_arithmetic
+      implicit none
+      REAL, INTENT(in)                       :: q(:), qis(:), qmt(:,:), qvac(:,:), qtot, qistot
+      character(len=*), intent(in), optional :: hint
+      integer                                :: n_mt
+      
+
+      if(present(hint)) write (*,*) "DEN of ", hint
+      write (*,*) "q   = ", q
+      write (*,*) "qis = ", qis
+
+      write (*,*) "qmt"
+      do n_mt = 1,size(qmt, dim=1)
+         write (*,*) "mt = ", n_mt, qmt(n_mt,:)
+      enddo
+
+      if(.not. any(ieee_is_nan(qvac))) then
+         write (*, *) "qvac",    qvac
+      endif
+      write (*, *) "qtot",    qtot
+      write (*, *) "qis_tot", qistot
+      write (*, *) "-------------------------"
+   END SUBROUTINE print_cdn_inte
 END MODULE m_cdntot

cdn/cdnval.F90

@@ -88,14 +88,14 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
 #endif
 
    ! Local Scalars
-   INTEGER :: ikpt,jsp_start,jsp_end,ispin,jsp
+   INTEGER :: ikpt,ikpt_i,jsp_start,jsp_end,ispin,jsp
    INTEGER :: iErr,nbands,noccbd,iType
-   INTEGER :: skip_t,skip_tt,nStart,nEnd,nbasfcn
+   INTEGER :: skip_t,skip_tt,nbasfcn
    LOGICAL :: l_orbcomprot, l_real, l_dosNdir, l_corespec
 
    ! Local Arrays
-   REAL,    ALLOCATABLE :: we(:)
-   REAL,    ALLOCATABLE :: eig(:)
+   REAL,ALLOCATABLE :: we(:),eig(:)
+   INTEGER,ALLOCATABLE :: ev_list(:)
    REAL,    ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:),flo(:,:,:,:) ! radial functions
 
    TYPE (t_lapw)             :: lapw
@@ -171,36 +171,28 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
 
    skip_tt = dot_product(enpara%skiplo(:atoms%ntype,jspin),atoms%neq(:atoms%ntype))
    IF (noco%l_soc.OR.noco%l_noco) skip_tt = 2 * skip_tt
-   ALLOCATE (we(MAXVAL(cdnvalJob%noccbd(:))))
-   ALLOCATE (eig(MAXVAL(cdnvalJob%noccbd(:))))
-   jsp = MERGE(1,jspin,noco%l_noco)
 
-   DO ikpt = cdnvalJob%ikptStart, cdnvalJob%nkptExtended, cdnvalJob%ikptIncrement
+   jsp = MERGE(1,jspin,noco%l_noco)
 
-      IF (ikpt.GT.kpts%nkpt) THEN
-#ifdef CPP_MPI
-         CALL MPI_BARRIER(mpi%mpi_comm,iErr) ! Synchronizes the RMA operations
-#endif
-         EXIT
-      END IF
+   DO ikpt_i = 1,size(cdnvalJob%k_list)
+      ikpt=cdnvalJob%k_list(ikpt_i)
 
       CALL lapw%init(input,noco, kpts,atoms,sym,ikpt,cell,.false., mpi)
       skip_t = skip_tt
-      noccbd = cdnvalJob%noccbd(ikpt)
-      nStart = cdnvalJob%nStart(ikpt)
-      nEnd = cdnvalJob%nEnd(ikpt)
-      we(1:noccbd) = cdnvalJob%weights(1:noccbd,ikpt)
-      eig(1:noccbd) = results%eig(nStart:nEnd,ikpt,jsp)
+      ev_list=cdnvaljob%compact_ev_list(ikpt_i,banddos%dos)
+      noccbd = SIZE(ev_list)
+      we  = cdnvalJob%weights(ev_list,ikpt)
+      eig = results%eig(ev_list,ikpt,jsp)
 
       IF (cdnvalJob%l_evp) THEN
-         IF (nStart > skip_tt) skip_t = 0
-         IF (nEnd <= skip_tt) skip_t = noccbd
-         IF ((nStart <= skip_tt).AND.(nEnd > skip_tt)) skip_t = mod(skip_tt,noccbd)
+         IF (minval(ev_list) > skip_tt) skip_t = 0
+         IF (maxval(ev_list) <= skip_tt) skip_t = noccbd
+         IF ((minval(ev_list) <= skip_tt).AND.(maxval(ev_list) > skip_tt)) skip_t = mod(skip_tt,noccbd)
       END IF
 
       nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
       CALL zMat%init(l_real,nbasfcn,noccbd)
-      CALL read_eig(eig_id,ikpt,jsp,n_start=nStart,n_end=nEnd,neig=nbands,zmat=zMat)
+      CALL read_eig(eig_id,ikpt,jsp,list=ev_list,neig=nbands,zmat=zMat)
 #ifdef CPP_MPI
       CALL MPI_BARRIER(mpi%mpi_comm,iErr) ! Synchronizes the RMA operations
 #endif

cdn_mt/abcof.F90

@@ -144,11 +144,11 @@ CONTAINS
                    ENDIF
 
                    IF (noco%l_noco) THEN
-                      !--->            generate the complex conjgates of the spinors (chi)
-                      ccchi(1,1) = CONJG( EXP(-ImagUnit*noco%alph(n)/2)*COS(noco%beta(n)/2))
-                      ccchi(1,2) = CONJG(-EXP(-ImagUnit*noco%alph(n)/2)*SIN(noco%beta(n)/2))
-                      ccchi(2,1) = CONJG( EXP( ImagUnit*noco%alph(n)/2)*SIN(noco%beta(n)/2))
-                      ccchi(2,2) = CONJG( EXP( ImagUnit*noco%alph(n)/2)*COS(noco%beta(n)/2))
+                      !--->            generate the spinors (chi)
+                      ccchi(1,1) =  EXP(ImagUnit*noco%alph(n)/2)*COS(noco%beta(n)/2)
+                      ccchi(1,2) = -EXP(ImagUnit*noco%alph(n)/2)*SIN(noco%beta(n)/2)
+                      ccchi(2,1) =  EXP(-ImagUnit*noco%alph(n)/2)*SIN(noco%beta(n)/2)
+                      ccchi(2,2) =  EXP(-ImagUnit*noco%alph(n)/2)*COS(noco%beta(n)/2)
                       IF (noco%l_ss) THEN
                          !--->              the coefficients of the spin-down basis functions are
                          !--->              stored in the second half of the eigenvector
@@ -188,7 +188,7 @@ CONTAINS
                       wronk = usdus%uds(l,n,jspin)*usdus%dus(l,n,jspin) - usdus%us(l,n,jspin)*usdus%duds(l,n,jspin)
                       IF (apw(l,n)) THEN
                          fj(l) = 1.0*const * fj(l)/usdus%us(l,n,jspin)
-                         dfj(l) = 0.0d0
+                         dfj(l) = 0.0
                       ELSE
                          dfj(l) = const* (usdus%dus(l,n,jspin)*fj(l)-df*usdus%us(l,n,jspin))/wronk
                          fj(l) = const* (df*usdus%uds(l,n,jspin)-fj(l)*usdus%duds(l,n,jspin))/wronk

cdn_mt/abcof3.F90

@@ -107,7 +107,7 @@ CONTAINS
              wronk = usdus%uds(l,n,jspin)*usdus%dus(l,n,jspin)-usdus%us(l,n,jspin)*usdus%duds(l,n,jspin) !Wronski determinante
              IF (apw(l,n)) THEN
                 fj(l) = 1.0*const * fj(l)/usdus%us(l,n,jspin)
-                dfj(l) = 0.0d0
+                dfj(l) = 0.0
              ELSE
                 dfj(l) = const* (usdus%dus(l,n,jspin)*fj(l)-df*usdus%us(l,n,jspin))/wronk
                 fj(l) = const* (df*usdus%uds(l,n,jspin)-fj(l)*usdus%duds(l,n,jspin))/wronk

cdn_mt/orbMagMoms.f90

@@ -31,16 +31,21 @@ SUBROUTINE orbMagMoms(input,atoms,noco,clmom)
    WRITE (6,FMT=9020)
    CALL openXMLElement('orbitalMagneticMomentsInMTSpheres',(/'units'/),(/'muBohr'/))
    DO iType = 1, atoms%ntype
-      IF (noco%l_noco) THEN
-         thetai = noco%beta(iType)
-         phii   = noco%alph(iType)
-      END IF
-
       ! magn. moment(-)
       slxmom = clmom(1,iType,1)+clmom(1,iType,2)
       slymom = clmom(2,iType,1)+clmom(2,iType,2)
       slmom =  clmom(3,iType,1)+clmom(3,iType,2)
 
+      IF (noco%l_noco) THEN
+         thetai = noco%beta(iType)
+         phii   = noco%alph(iType)
+
+         !Fix of sign of moment in first variation calculations. Perhaps it would be better to understand this :-(
+         !slxmom=-1*slxmom
+         slymom=-1*slymom
+         !slmom=-1*slmom
+      END IF
+
       ! rotation: orbital moment || spin moment (extended to incude phi - hopefully)
       slmom   = cos(thetai)*slmom + sin(thetai)*(cos(phii)*slxmom + sin(phii)*slymom)
       clmom(3,iType,1) = cos(thetai)*clmom(3,iType,1) + &

cmake/compileenv.txt

@@ -16,9 +16,20 @@ elseif (EXISTS ${CMAKE_SOURCE_DIR}/version)
   file(READ ${CMAKE_SOURCE_DIR}/version git_describe)
 endif()
 
-#normalize the strings
-string(STRIP  ${git_hash} git_hash)
-string(STRIP ${git_describe} git_describe)
-string(STRIP ${git_branch} git_branch)
-
+#normalize the strings, fix for problems in git commands above
+if (git_hash)
+  string(STRIP  ${git_hash} git_hash)
+else()
+   set(git_hash unkown)
+endif()
+if (git_describe)
+  string(STRIP ${git_describe} git_describe)
+else()
+   set(git_describe unkown)
+endif()
+if (git_branch)   
+  string(STRIP ${git_branch} git_branch)
+else()  
+   set(git_branch unkown)
+endif()
 file(GENERATE OUTPUT ${CMAKE_SOURCE_DIR}/init/compileinfo.h CONTENT "gitdesc=\"${git_describe}\"\ncompile_date=\"${compile_time}\"\ncompile_user=\"${compile_user}\"\ncompile_host=\"${compile_host}\"\ngitbranch=\"${git_branch}\"\ngithash=\"${git_hash}\"\ncompile_flags=\"${CMAKE_Fortran_FLAGS}\"\nlink_flags=\"${FLEUR_LIBRARIES}\"\n")

cmake/compilerflags.cmake

@@ -17,8 +17,9 @@ if (${CMAKE_Fortran_COMPILER_ID} MATCHES "Intel")
       set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -mkl -qopenmp -assume byterecl")
    endif()     
    set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -xHost -O2 -g")
+   #set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -xMIC-AVX512 -O2")
    if (${CMAKE_Fortran_COMPILER_VERSION} VERSION_LESS "19.0.0.0")
-       set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -C -traceback -O0 -g -ftrapuv -check uninit -check pointers -DCPP_DEBUG -warn=all")
+       set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -C -traceback -O0 -g -ftrapuv -check uninit -check pointers -DCPP_DEBUG -warn all")
    else()
        set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -CB -traceback -O0 -g -ftrapuv -check uninit -check pointers -DCPP_DEBUG")
    endif()

cmake/machines/HAZELHEN.sh

+export CRAYPE_LINK_TYPE=dynamic

cmake/tests/test_MPI.cmake

@@ -16,9 +16,13 @@ if (DEFINED CLI_FLEUR_USE_MPI)
 	   set(FLEUR_USE_MPI FALSE)
        endif()
    endif()	   
-endif()  
+endif()
+
+if ( "$ENV{I_MPI_ROOT}" MATCHES ".*_2019.*")
+   message("It looks like you are using the IntelMPI2019 library. This is buggy. \n You might want to do 'export MPIR_CVAR_CH4_OFI_ENABLE_RMA=0' to enable running FLEUR")
+endif()
 
-if (FLEUR_USE_MPI)
+If (FLEUR_USE_MPI)
    set(FLEUR_MPI_DEFINITIONS ${FLEUR_MPI_DEFINITIONS} "CPP_MPI")
    set(FLEUR_USE_SERIAL FALSE)
 else()
@@ -27,4 +31,4 @@ endif()
 
 if (DEFINED CLI_FLEUR_USE_SERIAL)
   set(FLEUR_USE_SERIAL CLI_FLEUR_USE_SERIAL)
-endif()
\ No newline at end of file
+endif()

configure.sh

@@ -23,7 +23,8 @@ echo "------------ Welcome to the FLEUR configuration script -------------"
 
 
 #Check if we are using the git version and update if pull was used as an argument
-if test -d $DIR/.git 
+if test -d $DIR/.git
+then
     #Check if hook is installed and install it if needed
     if test -h $DIR/.git/hooks/pre-commit
     then
@@ -33,13 +34,12 @@ if test -d $DIR/.git
         ln -s $DIR/tests/git-hooks/pre-commit $DIR/.git/hooks
         echo "Git version found, hook installed"
     fi
-then
-   if [ $gitupdate -gt 0 ] 
-   then
+    if [ $gitupdate -gt 0 ] 
+    then
        cd $DIR 
        git pull
        cd -
-   fi
+    fi
 fi
 
 

diagonalization/available_solvers.F90

@@ -48,8 +48,11 @@ MODULE m_available_solvers
 #else
   INTEGER,PARAMETER:: diag_elpa_1node=-14
 #endif 
+#ifdef CPP_SCALAPACK
   INTEGER,PARAMETER:: diag_debugout=201
-  
+#else
+  INTEGER,PARAMETER:: diag_debugout=20
+#endif  
   INTEGER,PARAMETER::diag_all_solver(9)=(/diag_elpa,diag_elemental,diag_scalapack,diag_magma,diag_chase,diag_cusolver,diag_lapack,diag_elpa_1node,diag_debugout/)
   
 CONTAINS

diagonalization/elpa.F90

@@ -192,19 +192,19 @@ CONTAINS
        n_col = indxl2g(i,     nb, hmat%blacsdata%mycol, 0, hmat%blacsdata%npcol)
        n_row = numroc (n_col, nb, hmat%blacsdata%myrow, 0, hmat%blacsdata%nprow)
        IF (hmat%l_real) THEN
-          hmat%data_r(n_row+1:hmat%matsize1,i) = 0.d0 
+          hmat%data_r(n_row+1:hmat%matsize1,i) = 0.0 
        ELSE
-          hmat%data_c(n_row+1:hmat%matsize1,i) = 0.d0
+          hmat%data_c(n_row+1:hmat%matsize1,i) = 0.0
        ENDIF
     ENDDO
  
     ! Use the ev_dist array to store the calculated values for the lower part.
     IF (hmat%l_real) THEN
-       CALL pdtran(hmat%global_size1,hmat%global_size1,1.d0,hmat%data_r,1,1,&
-                         hmat%blacsdata%blacs_desc,0.d0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc)
+       CALL pdtran(hmat%global_size1,hmat%global_size1,1.0,hmat%data_r,1,1,&
+                         hmat%blacsdata%blacs_desc,0.0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc)
     ELSE
-       CALL pztranc(hmat%global_size1,hmat%global_size2,cmplx(1.d0,0.d0),hmat%data_c,1,1,&
-                         hmat%blacsdata%blacs_desc,cmplx(0.d0,0.d0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc)
+       CALL pztranc(hmat%global_size1,hmat%global_size2,cmplx(1.0,0.0),hmat%data_c,1,1,&
+                         hmat%blacsdata%blacs_desc,cmplx(0.0,0.0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc)
     ENDIF
 
     ! Copy the calculated values to the lower part of the H matrix
@@ -259,11 +259,11 @@ CONTAINS
 
     ! 2b. tmp2 = ev_dist**T
     IF (hmat%l_real) THEN
-       CALL pdtran(ev_dist%global_size1,ev_dist%global_size1,1.d0,ev_dist%data_r,1,1,&
-                          ev_dist%blacsdata%blacs_desc,0.d0,tmp2_r,1,1,ev_dist%blacsdata%blacs_desc)
+       CALL pdtran(ev_dist%global_size1,ev_dist%global_size1,1.0,ev_dist%data_r,1,1,&
+                          ev_dist%blacsdata%blacs_desc,0.0,tmp2_r,1,1,ev_dist%blacsdata%blacs_desc)
     ELSE
        CALL pztranc(ev_dist%global_size1,ev_dist%global_size1,cmplx(1.0,0.0),ev_dist%data_c,1,1,&
-                          ev_dist%blacsdata%blacs_desc,cmplx(0.d0,0.d0),tmp2_c,1,1,ev_dist%blacsdata%blacs_desc)
+                          ev_dist%blacsdata%blacs_desc,cmplx(0.0,0.0),tmp2_c,1,1,ev_dist%blacsdata%blacs_desc)
     ENDIF
 
     ! 2c. A =  U**-T * tmp2 ( = U**-T * Aorig * U**-1 )
@@ -307,11 +307,11 @@ CONTAINS
     ! Set lower half from upper half
 
     IF (hmat%l_real) THEN
-       CALL pdtran(hmat%global_size1,hmat%global_size1,1.d0,hmat%data_r,1,1,&
-                          hmat%blacsdata%blacs_desc,0.d0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc)
+       CALL pdtran(hmat%global_size1,hmat%global_size1,1.0,hmat%data_r,1,1,&
+                          hmat%blacsdata%blacs_desc,0.0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc)
     ELSE
        CALL pztranc(hmat%global_size1,hmat%global_size1,cmplx(1.0,0.0),hmat%data_c,1,1,&
-                          hmat%blacsdata%blacs_desc,cmplx(0.d0,0.d0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc)
+                          hmat%blacsdata%blacs_desc,cmplx(0.0,0.0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc)
     ENDIF
 
 
@@ -396,11 +396,11 @@ CONTAINS
 
     ! mult_ah_b_complex needs the transpose of U**-1, thus tmp2 = (U**-1)**T
     IF (hmat%l_real) THEN
-       CALL pdtran(smat%global_size1,smat%global_size1,1.d0,smat%data_r,1,1,&
-                         smat%blacsdata%blacs_desc,0.d0,tmp2_r,1,1,smat%blacsdata%blacs_desc)
+       CALL pdtran(smat%global_size1,smat%global_size1,1.0,smat%data_r,1,1,&
+                         smat%blacsdata%blacs_desc,0.0,tmp2_r,1,1,smat%blacsdata%blacs_desc)
     ELSE
-       CALL pztranc(smat%global_size1,smat%global_size1,cmplx(1.d0,0.d0),smat%data_c,1,1,&
-                         smat%blacsdata%blacs_desc,cmplx(0.d0,0.d0),tmp2_c,1,1,smat%blacsdata%blacs_desc)
+       CALL pztranc(smat%global_size1,smat%global_size1,cmplx(1.0,0.0),smat%data_c,1,1,&
+                         smat%blacsdata%blacs_desc,cmplx(0.0,0.0),tmp2_c,1,1,smat%blacsdata%blacs_desc)
     ENDIF
 
 #if defined (CPP_ELPA_201705003)

diagonalization/lapack_diag.F90

@@ -23,6 +23,7 @@ IMPLICIT NONE
     REAL               :: dumrwork(1),abstol
     COMPLEX            :: dumwork(1)
     REAL,external      :: dlamch
+    REAL               :: eigTemp(hmat%matsize1)
 
     
     ALLOCATE(t_mat::zmat)
@@ -31,24 +32,25 @@ IMPLICIT NONE
     IF (hmat%l_real) THEN
        ALLOCATE(iwork(5*hmat%matsize1),ifail(hmat%matsize1))
        CALL dsygvx(1,'V','I','U', hmat%matsize1,hmat%data_r,SIZE(hmat%data_r,1),smat%data_r,SIZE(smat%data_r,1),&
-            0.0,0.0,1,ne,abstol,m,eig,zmat%data_r,SIZE(zmat%data_r,1),dumrwork,-1, iwork, ifail, info)
+            0.0,0.0,1,ne,abstol,m,eigTemp,zmat%data_r,SIZE(zmat%data_r,1),dumrwork,-1, iwork, ifail, info)
        lwork=dumrwork(1)
        ALLOCATE(rwork(lwork))
        IF (info.NE.0) CALL judft_error("Diagonalization via LAPACK failed (Workspace)",no=info)
        CALL dsygvx(1,'V','I','U', hmat%matsize1,hmat%data_r,SIZE(hmat%data_r,1),smat%data_r,SIZE(smat%data_r,1),&
-            0.0,0.0,1,ne,abstol,m,eig,zmat%data_r,SIZE(zmat%data_r,1),rwork, lwork, iwork, ifail, info)
+            0.0,0.0,1,ne,abstol,m,eigTemp,zmat%data_r,SIZE(zmat%data_r,1),rwork, lwork, iwork, ifail, info)
     ELSE
        ALLOCATE(rwork(7*hmat%matsize1),iwork(5*hmat%matsize1),ifail(hmat%matsize1))
        !Do a workspace query
        CALL zhegvx(1,'V','I','U',hmat%matsize1,hmat%data_c,SIZE(hmat%data_c,1),smat%data_c,SIZE(smat%data_c,1),&
-            0.0,0.0,1,ne,abstol,m,eig,zmat%data_c,SIZE(zmat%data_c,1),dumwork,-1,rwork,iwork,ifail,info)
+            0.0,0.0,1,ne,abstol,m,eigTemp,zmat%data_c,SIZE(zmat%data_c,1),dumwork,-1,rwork,iwork,ifail,info)
        lwork=dumwork(1)
        ALLOCATE(work(lwork))
        IF (info.NE.0) CALL judft_error("Diagonalization via LAPACK failed (Workspace)",no=info)
        !Perform diagonalization
        CALL zhegvx(1,'V','I','U',hmat%matsize1,hmat%data_c,SIZE(hmat%data_c,1),smat%data_c,SIZE(smat%data_c,1),&
-            0.0,0.0,1,ne,abstol,m,eig,zmat%data_c,SIZE(zmat%data_c,1),work,lwork,rwork,iwork,ifail,info)
+            0.0,0.0,1,ne,abstol,m,eigTemp,zmat%data_c,SIZE(zmat%data_c,1),work,lwork,rwork,iwork,ifail,info)
     ENDIF
+    eig(:SIZE(eig)) = eigTemp(:SIZE(eig))
     IF (info.NE.0) CALL judft_error("Diagonalization via LAPACK failed(zhegvx/dsygvx)",no=info)
     IF (m.NE.ne) CALL judft_error("Diagonalization via LAPACK failed failed without explicit errorcode.")
   END SUBROUTINE lapack_diag

dos/evaldos.f90

@@ -119,8 +119,8 @@
       ENDDO
  
       IF ( l_mcd ) THEN ! create an energy grid for mcd-spectra
-        e_lo =  9.9d+9 
-        e_up = -9.9d+9     
+        e_lo =  9.9*10.0**9 
+        e_up = -9.9*10.0**9     
         DO jspin = 1,input%jspins
           DO n = 1,atoms%ntype
             DO icore = 1 , mcd%ncore(n)

eels/corespec.f90

@@ -12,11 +12,11 @@ module m_corespec
 
 ! PARAMETERS
 
-  complex, parameter :: cone = cmplx(1.d0,0.d0)
-  complex, parameter :: cimu = cmplx(0.d0,1.d0)
-  real, parameter :: alpha = 7.29735257d-3
-  real, parameter :: mec2 = 0.51099891d6
-  real, parameter :: ecoredeep = 0.5d0
+  complex, parameter :: cone = cmplx(1.0,0.0)
+  complex, parameter :: cimu = cmplx(0.0,1.0)
+  real, parameter :: alpha = 7.29735257e-3
+  real, parameter :: mec2 = 0.51099891e6
+  real, parameter :: ecoredeep = 0.5
 
   integer, parameter :: edgel(11) = (/0,1,1,2,2,3,3,4,4,5,5/)
   integer, parameter :: edgej(11) = (/1,1,3,3,5,5,7,7,9,9,11/)

eels/corespec_io.f90

@@ -42,10 +42,10 @@ MODULE m_corespec_io
     csi%edge = ""
     csi%edgeidx = 0
     csi%lx = -1
-    csi%ek0 = 0.d0
-    csi%emn = -2.d0
-    csi%emx = 20.d0
-    csi%ein = 0.1d0
+    csi%ek0 = 0.0
+    csi%emn = -2.0
+    csi%emx = 20.0
+    csi%ein = 0.1
 
     csi%nqphi = 12
     csi%nqr = 20
@@ -182,12 +182,12 @@ MODULE m_corespec_io
            &"maximum l: ","csi%lx = ",csi%lx,"will be used"
 
 ! csi%ek0
-    if(csi%ek0.le.0.d0) then
+    if(csi%ek0.le.0.0) then
       write(*,csmsgs)  trim(smeno),"found csi%ek0 <= 0.0 !"//csmsgerr ; stop
     endif
-    csi%ek0 = csi%ek0*1000.d0 ! conversion from keV to eV
-    csv%gamma = 1.d0+csi%ek0/mec2
-    csv%beta = sqrt(1.d0-1.d0/(csv%gamma**2))
+    csi%ek0 = csi%ek0*1000.0 ! conversion from keV to eV
+    csv%gamma = 1.0+csi%ek0/mec2
+    csv%beta = sqrt(1.0-1.0/(csv%gamma**2))
     if(csi%verb.eq.1) then
       write(*,csmsgses)  trim(smeno),&
            &"kinetic energy of incoming electrons: ","csi%ek0 = ",csi%ek0,&
@@ -204,8 +204,8 @@ MODULE m_corespec_io
     if(csi%emn.gt.csi%emx) then
       write(*,csmsgs)  trim(smeno),"found csi%emn > csi%emx !"//csmsgerr ; stop
     endif
-    if(csi%ein.le.0.d0) then
-      write(*,csmsgs)  trim(smeno),"found csi%ein <= 0.d0 !"//csmsgerr ; stop
+    if(csi%ein.le.0.0) then
+      write(*,csmsgs)  trim(smeno),"found csi%ein <= 0.0 !"//csmsgerr ; stop
     endif
     if(((csi%emx-csi%emn)/csi%ein)-int((csi%emx-csi%emn)/csi%ein).ne.0) then
       write(*,csmsgs)  trim(smeno),&
@@ -216,7 +216,7 @@ MODULE m_corespec_io
     csv%egrid = (/(csi%emn+csi%ein*dble(i), i = 0,csv%nex)/)
     csv%nen = 0
 !!$    do i = 0,csv%nex
-!!$      if(csv%egrid(i).ge.0.d0) then
+!!$      if(csv%egrid(i).ge.0.0) then
 !!$        csv%nen = i
 !!$        exit
 !!$      endif
@@ -240,9 +240,9 @@ MODULE m_corespec_io
            &csv%nen,"will be used"
 
     if(.not.allocated(csv%eedge)) allocate(csv%eedge(csv%nljc))
-    csv%eedge = 0.d0
+    csv%eedge = 0.0
     if(.not.allocated(csv%occ)) allocate(csv%occ(csv%nljc))
-    csv%occ = 0.d0
+    csv%occ = 0.0
 
     l_cs = .true.
 

eigen/eigen.F90

@@ -18,11 +18,10 @@ CONTAINS
    !>
    !> The matrices generated and diagonalized here are of type m_mat as defined in m_types_mat.
    !>@author D. Wortmann
-   SUBROUTINE eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,sym,kpts,DIMENSION,vacuum,input,&
+   SUBROUTINE eigen(mpi,stars,sphhar,atoms,xcpot,sym,kpts,DIMENSION,vacuum,input,&
                     cell,enpara,banddos,noco,oneD,hybrid,iter,eig_id,results,inden,v,vx)
 
 #include"cpp_double.h"
-      USE m_constants, ONLY : pi_const,sfp_const
       USE m_types
       USE m_eigen_hssetup
       USE m_pot_io
@@ -51,7 +50,6 @@ CONTAINS
       TYPE(t_oneD),INTENT(IN)      :: oneD
       TYPE(t_hybrid),INTENT(INOUT) :: hybrid
       TYPE(t_enpara),INTENT(INOUT) :: enpara
-      TYPE(t_obsolete),INTENT(IN)  :: obsolete
       TYPE(t_input),INTENT(IN)     :: input
       TYPE(t_vacuum),INTENT(IN)    :: vacuum
       TYPE(t_noco),INTENT(IN)      :: noco
@@ -77,20 +75,18 @@ CONTAINS
 
       ! Local Scalars
       INTEGER jsp,nk,nred,ne_all,ne_found
-      INTEGER ne,lh0
+      INTEGER ne, nk_i
       INTEGER isp,i,j,err
       LOGICAL l_wu,l_file,l_real,l_zref
       INTEGER :: solver=0
       ! Local Arrays
-      INTEGER              :: ierr(3)
+      INTEGER              :: ierr
       INTEGER              :: neigBuffer(kpts%nkpt,input%jspins)
 
       COMPLEX              :: unfoldingBuffer(SIZE(results%unfolding_weights,1),kpts%nkpt,input%jspins) ! needed for unfolding bandstructure mpi case
 
-      INTEGER, PARAMETER   :: lmaxb = 3
       REAL,    ALLOCATABLE :: bkpt(:)
-      REAL,    ALLOCATABLE :: eig(:)
-      COMPLEX, ALLOCATABLE :: vs_mmp(:,:,:,:)
+      REAL,    ALLOCATABLE :: eig(:), eigBuffer(:,:,:)
 
       INTEGER                   :: jsp_m, i_kpt_m, i_m
 
@@ -106,7 +102,9 @@ CONTAINS
       character(len=300)        :: errmsg
 
       call ud%init(atoms,input%jspins)
-      ALLOCATE (eig(DIMENSION%neigd),bkpt(3))
+      ALLOCATE(eig(DIMENSION%neigd))
+      ALLOCATE(bkpt(3))
+      ALLOCATE(eigBuffer(DIMENSION%neigd,kpts%nkpt,input%jspins))
 
       l_real=sym%invs.AND..NOT.noco%l_noco
 
@@ -129,11 +127,12 @@ CONTAINS
       neigBuffer = 0
       results%neig = 0
       results%eig = 1.0e300
+      eigBuffer = 1.0e300
       unfoldingBuffer = CMPLX(0.0,0.0)
 
       DO jsp = 1,MERGE(1,input%jspins,noco%l_noco)
-         k_loop:DO nk = mpi%n_start,kpts%nkpt,mpi%n_stride
-
+         k_loop:DO nk_i = 1,size(mpi%k_list)
+            nk=mpi%k_list(nk_i)
             ! Set up lapw list
             CALL lapw%init(input,noco, kpts,atoms,sym,nk,cell,l_zref, mpi)
             call timestart("Setup of H&S matrices")
@@ -234,6 +233,7 @@ CONTAINS
                 !        Mai 2019                 U. Alekseeva      
                 CALL write_eig(eig_id, nk,jsp,ne_found,ne_all,&
                            eig(:ne_all),n_start=mpi%n_size,n_end=mpi%n_rank,zMat=zMat)
+                eigBuffer(:ne_all,nk,jsp) = eig(:ne_all)
             ELSE
                 CALL write_eig(eig_id, nk,jsp,ne_found,&
                            n_start=mpi%n_size,n_end=mpi%n_rank,zMat=zMat)
@@ -246,6 +246,8 @@ CONTAINS
             CALL timestop("EV output")
 
             IF (banddos%unfoldband) THEN
+               IF(modulo (kpts%nkpt,mpi%n_size).NE.0) call juDFT_error("number kpts needs to be multiple of number mpi threads",& 
+                   hint=errmsg, calledby="eigen.F90")
                CALL calculate_plot_w_n(banddos,cell,kpts,smat_unfold,zMat,lapw,nk,jsp,eig,results,input,atoms,unfoldingBuffer,mpi)
                CALL smat_unfold%free()
                DEALLOCATE(smat_unfold, stat=dealloc_stat, errmsg=errmsg)
@@ -264,24 +266,15 @@ CONTAINS
        CALL MPI_ALLREDUCE(unfoldingBuffer,results%unfolding_weights,SIZE(results%unfolding_weights,1)*SIZE(results%unfolding_weights,2)*SIZE(results%unfolding_weights,3),CPP_MPI_COMPLEX,MPI_SUM,mpi%mpi_comm,ierr)
       END IF
       CALL MPI_ALLREDUCE(neigBuffer,results%neig,kpts%nkpt*input%jspins,MPI_INTEGER,MPI_SUM,mpi%mpi_comm,ierr)
+      CALL MPI_ALLREDUCE(eigBuffer(:dimension%neigd,:,:),results%eig(:dimension%neigd,:,:),dimension%neigd*kpts%nkpt*input%jspins,MPI_DOUBLE_PRECISION,MPI_MIN,mpi%mpi_comm,ierr)
       CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
+
 #else
       results%neig(:,:) = neigBuffer(:,:)
+      results%eig(:dimension%neigd,:,:) = eigBuffer(:dimension%neigd,:,:)
       results%unfolding_weights(:,:,:) = unfoldingBuffer(:,:,:)
 #endif
 
-      ! Sorry for the following strange workaround to fill the results%eig array.
-      ! At some point someone should have a closer look at how the eigenvalues are
-      ! distributed and fill the array without using the eigenvalue-IO.
-      DO jsp = 1,MERGE(1,input%jspins,noco%l_noco)
-         DO nk = 1,kpts%nkpt
-            CALL read_eig(eig_id,nk,jsp,results%neig(nk,jsp),results%eig(:,nk,jsp))
-#ifdef CPP_MPI
-            CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
-#endif
-         END DO
-      END DO
-
       !IF (hybrid%l_hybrid.OR.hybrid%l_calhf) CALL close_eig(eig_id)
 
       IF( input%jspins .EQ. 1 .AND. hybrid%l_hybrid ) THEN

eigen/hs_int.F90

@@ -48,7 +48,7 @@ CONTAINS
                    ii=-1*ii
                    in = stars%ig(ii(1),ii(2),ii(3))
                    IF (in.EQ.0) CYCLE
-                   th = stars%rgphs(ii(1),ii(2),ii(3))*CONJG(vpw(in,3))           
+                   th = stars%rgphs(ii(1),ii(2),ii(3))*conjg(vpw(in,3))           
                    ts=0.0
                 ELSEIF(ispin==2.and.jspin==1) THEN
                 !   ii = -1*ii

eigen/hsmt.F90

@@ -112,13 +112,13 @@ CONTAINS
              CALL hsmt_spinor(ispin,n,noco,chi)
              DO iintsp=1,2
                 DO jintsp=1,2
-                   CALL hsmt_sph(n,atoms,mpi,ispin,input,noco,iintsp,jintsp,chi(jintsp,iintsp),&
+                   CALL hsmt_sph(n,atoms,mpi,ispin,input,noco,iintsp,jintsp,chi(iintsp,jintsp),&
                                  lapw,enpara%el0,td%e_shift(n,ispin),usdus,fj(:,0:,ispin,:),gj(:,0:,ispin,:),&
                                  smat(iintsp,jintsp),hmat(iintsp,jintsp))
-                   CALL hsmt_nonsph(n,mpi,sym,atoms,ispin,iintsp,jintsp,chi(jintsp,iintsp),noco,cell,&
+                   CALL hsmt_nonsph(n,mpi,sym,atoms,ispin,iintsp,jintsp,chi(iintsp,jintsp),noco,cell,&
                                     lapw,td,fj(:,0:,ispin,:),gj(:,0:,ispin,:),hmat(iintsp,jintsp))
                    CALL hsmt_lo(input,atoms,sym,cell,mpi,noco,lapw,usdus,td,fj(:,0:,ispin,:),gj(:,0:,ispin,:),&
-                          n,chi(jintsp,iintsp),ispin,iintsp,jintsp,hmat(iintsp,jintsp),smat(iintsp,jintsp))
+                          n,chi(iintsp,jintsp),ispin,iintsp,jintsp,hmat(iintsp,jintsp),smat(iintsp,jintsp))
                 ENDDO
              ENDDO
           ENDIF

eigen/hsmt_distspin.F90

@@ -17,7 +17,7 @@ CONTAINS
 
     DO iintsp=1,2
        DO jintsp=1,2
-          mat(jintsp,iintsp)%data_c(:,:)=chi(iintsp,jintsp)*mat_tmp%data_c(:,:)+mat(jintsp,iintsp)%data_c(:,:)
+          mat(jintsp,iintsp)%data_c(:,:)=chi(jintsp,iintsp)*mat_tmp%data_c(:,:)+mat(jintsp,iintsp)%data_c(:,:)
        ENDDO
     ENDDO
   END SUBROUTINE hsmt_distspins

eigen/hsmt_fjgj.F90

@@ -44,7 +44,7 @@ CONTAINS
              gg = rk(k)*gb(l)
              IF ( apw(l) ) THEN
                 fj(k,l,ispin) = 1.0*con1 * ff / us(l,ispin)
-                gj(k,l,ispin) = 0.0d0
+                gj(k,l,ispin) = 0.0
              ELSE
                 IF (l_flag) THEN
                    DO jspin = 1, jspins
@@ -167,7 +167,7 @@ CONTAINS
              gg = lapw%rk(k,intspin)*gb(l)
              IF ( apw(l) ) THEN
                 fj(k,l,ispin,intspin) = 1.0*con1 * ff / usdus%us(l,n,ispin)
-                gj(k,l,ispin,intspin) = 0.0d0
+                gj(k,l,ispin,intspin) = 0.0
              ELSE
                 IF (noco%l_constr.or.l_socfirst) THEN
                    DO jspin = 1, input%jspins

eigen/hsmt_soc_offdiag.f90

@@ -37,8 +37,7 @@ CONTAINS
     COMPLEX:: chi(2,2,2,2),angso(lapw%nv(1),2,2)
     REAL, ALLOCATABLE :: plegend(:,:),dplegend(:,:)
     COMPLEX, ALLOCATABLE :: cph(:)
-  
-    
+   
     CALL timestart("offdiagonal soc-setup")
     
     DO l = 0,atoms%lmaxd
@@ -87,6 +86,7 @@ CONTAINS
        DO  l = 1,atoms%lmax(n)
           DO j1=1,2
              DO j2=1,2
+             !DO j2=j1,j1
                 DO kj = 1,ki
                    fct  =cph(kj) * dplegend(kj,l)*fl2p1(l)*(&
                         fj(ki,l,j1)*fj(kj,l,j2) *td%rsoc%rsopp(n,l,j1,j2) + &

eigen/hsmt_sph.F90

@@ -398,8 +398,8 @@ SUBROUTINE hsmt_sph_cpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_s
 
       !--->          update overlap and l-diagonal hamiltonian matrix
       kj_end = MIN(ki,lapw%nv(iintsp))
-      VecHelpS = 0.d0
-      VecHelpH = 0.d0
+      VecHelpS = 0.0
+      VecHelpH = 0.0
       DO  l = 0,atoms%lmax(n)
          fjkiln = fj(ki,l,jintsp)
          gjkiln = gj(ki,l,jintsp)

eigen/hsmt_spinor.F90

@@ -24,21 +24,21 @@ CONTAINS
 
     !--->       set up the spinors of this atom within global
     !--->       spin-coordinateframe
-    chi(1,1) =  exp(-ImagUnit*noco%alph(n)/2)*cos(noco%beta(n)/2)
-    chi(1,2) = -exp(-ImagUnit*noco%alph(n)/2)*sin(noco%beta(n)/2)
-    chi(2,1) =  exp(ImagUnit*noco%alph(n)/2)*sin(noco%beta(n)/2)
-    chi(2,2) =  exp(ImagUnit*noco%alph(n)/2)*cos(noco%beta(n)/2)
+    chi(1,1) =  exp(ImagUnit*noco%alph(n)/2)*cos(noco%beta(n)/2)
+    chi(1,2) = -EXP(ImagUnit*noco%alph(n)/2)*SIN(noco%beta(n)/2)
+    chi(2,1) =  EXP(-ImagUnit*noco%alph(n)/2)*SIN(noco%beta(n)/2)
+    chi(2,2) =  EXP(-ImagUnit*noco%alph(n)/2)*COS(noco%beta(n)/2)
     !--->       and determine the prefactors for the Hamitonian- and
     !--->       overlapp-matrix elements
     IF (isp<3) THEN
        isp1=isp
        isp2=isp
     ELSEIF(isp==3) THEN
-       isp1=1
-       isp2=2
-    ELSE
        isp1=2
        isp2=1
+    ELSE
+       isp1=1
+       isp2=2
     ENDIF
     
     chi_mat(1,1) = chi(1,isp1)*CONJG(chi(1,isp2))
@@ -77,10 +77,10 @@ CONTAINS
     
     !--->       set up the spinors of this atom within global
     !--->       spin-coordinateframe
-    chi(1,1) =  exp(-ImagUnit*noco%alph(n)/2)*cos(noco%beta(n)/2)
-    chi(1,2) = -exp(-ImagUnit*noco%alph(n)/2)*sin(noco%beta(n)/2)
-    chi(2,1) =  exp(ImagUnit*noco%alph(n)/2)*sin(noco%beta(n)/2)
-    chi(2,2) =  EXP(ImagUnit*noco%alph(n)/2)*COS(noco%beta(n)/2)
+    chi(1,1) =  exp(ImagUnit*noco%alph(n)/2)*cos(noco%beta(n)/2)
+    chi(1,2) = -EXP(ImagUnit*noco%alph(n)/2)*SIN(noco%beta(n)/2)
+    chi(2,1) =  EXP(-ImagUnit*noco%alph(n)/2)*SIN(noco%beta(n)/2)
+    chi(2,2) =  EXP(-ImagUnit*noco%alph(n)/2)*COS(noco%beta(n)/2)
 
     isigma_x=MATMUL(CONJG(TRANSPOSE(chi)), MATMUL(isigma(:,:,1),chi))
     isigma_y=MATMUL(CONJG(TRANSPOSE(chi)), MATMUL(isigma(:,:,2),chi))

eigen/old/hsmt_sph_new.F90

@@ -88,7 +88,7 @@ CONTAINS
                 gg = lapw%rk(k,iintsp)*gb(l)
         !        IF ( apw(l) ) THEN
         !           fj(k,l,n,iintsp) = 1.0*con1 * ff / usdus%us(l,n,isp)
-        !           gj(k,l,n,iintsp) = 0.0d0
+        !           gj(k,l,n,iintsp) = 0.0
         !        ELSE
                    !--->                 in a spin-spiral calculation fj and gj are needed
                    !--->                 both interstitial spin directions at the same time

eigen/vacfun.f90

@@ -53,7 +53,7 @@ CONTAINS
     REAL u(vacuum%nmzd,size(duz,1),input%jspins),ud(vacuum%nmzd,size(duz,1),input%jspins)
     REAL v(3),x(vacuum%nmzd), qssbti(2,2)
     !     ..
-    fac=MERGE(1.0,-1.0,jspin1<=jspin2)
+    fac=MERGE(1.0,-1.0,jspin1>=jspin2)
     ipot=MERGE(jspin1,3,jspin1==jspin2)
 
     tuuv=0.0;tudv=0.0;tddv=0.0;tduv=0.0

eigen_soc/abcof_soc.F90

@@ -139,7 +139,7 @@ CONTAINS
                       wronk = usdus%uds(l,n,jspin)*usdus%dus(l,n,jspin) - usdus%us(l,n,jspin)*usdus%duds(l,n,jspin)
                       IF (apw(l,n)) THEN
                          fj(l) = 1.0*const * fj(l)/usdus%us(l,n,jspin)
-                         dfj(l) = 0.0d0
+                         dfj(l) = 0.0
                       ELSE
                          dfj(l) = const* (usdus%dus(l,n,jspin)*fj(l)-df*usdus%us(l,n,jspin))/wronk
                          fj(l) = const* (df*usdus%uds(l,n,jspin)-fj(l)*usdus%duds(l,n,jspin))/wronk

eigen_soc/alineso.F90

@@ -97,7 +97,7 @@ CONTAINS
           WRITE(6,*) "Non-SOC ev for nk,jsp:",nk,jsp
           WRITE(6,"(6(f10.6,1x))") eig(:ne,jsp)
        ENDIF
-       CALL read_eig(eig_id,nk,jsp,n_start=1,n_end=ne,zmat=zmat(jsp))
+       CALL read_eig(eig_id,nk,jsp,list=[(i,i=1,ne)],zmat=zmat(jsp))
 
        ! write(*,*) 'process',irank,' reads ',nk
 
@@ -297,7 +297,7 @@ else
           IF (i1.EQ.1) nn = 0
           IF (i1.EQ.2) nn = nsz(1)
 
-          zhelp2(:,:) = 0.d0
+          zhelp2(:,:) = 0.0
           DO j = 1,nsize
              DO i = 1,nsz(jsp)
                 zhelp2(i,j) =  CONJG(hso(i+nn,j))
@@ -305,11 +305,11 @@ else
           ENDDO  ! j
 
           if (l_real) THEN
-             CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd,CMPLX(1.d0,0.d0),CMPLX(zmat(jsp)%data_r(:,:)),&
-                  zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.d0,0.d0), zso(1,1,jsp2),zmat(1)%matsize1)
+             CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd,CMPLX(1.0,0.0),CMPLX(zmat(jsp)%data_r(:,:)),&
+                  zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.0,0.0), zso(1,1,jsp2),zmat(1)%matsize1)
           else
-             CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd, CMPLX(1.d0,0.d0),zmat(jsp)%data_c(:,:),&
-                  zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.d0,0.d0), zso(:,:,jsp2),zmat(1)%matsize1)
+             CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd, CMPLX(1.0,0.0),zmat(jsp)%data_c(:,:),&
+                  zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.0,0.0), zso(:,:,jsp2),zmat(1)%matsize1)
           endif
 
        ENDDO    !isp

eigen_soc/eigenso.F90

@@ -49,24 +49,30 @@ CONTAINS
     TYPE(t_potden),INTENT(IN)     :: vTot
     TYPE(t_enpara),INTENT(IN)     :: enpara
     TYPE(t_results),INTENT(INOUT) :: results
+
+#ifdef CPP_MPI
+    INCLUDE 'mpif.h'
+#endif
+
     !     ..
     !     .. Scalar Arguments ..
     INTEGER, INTENT (IN) :: eig_id       
     !     ..
     !     ..
     !     .. Local Scalars ..
-    INTEGER i,j,nk,jspin,n ,l
+    INTEGER i,j,nk,nk_i,jspin,n ,l
     ! INTEGER n_loc,n_plus,i_plus,
-    INTEGER n_end,nsz,nmat,n_stride
+    INTEGER nsz,nmat,n_stride
     LOGICAL l_socvec   !,l_all
     INTEGER wannierspin
-    TYPE(t_usdus):: usdus
+    TYPE(t_usdus)        :: usdus
     !     ..
     !     .. Local Arrays..
     CHARACTER*3 chntype
 
     TYPE(t_rsoc) :: rsoc
-    REAL,    ALLOCATABLE :: eig_so(:) 
+    INTEGER, ALLOCATABLE :: neigBuffer(:,:)
+    REAL,    ALLOCATABLE :: eig_so(:), eigBuffer(:,:,:)
     COMPLEX, ALLOCATABLE :: zso(:,:,:)
 
     TYPE(t_mat)::zmat
@@ -112,30 +118,20 @@ CONTAINS
     !
 
 
-    ALLOCATE( eig_so(2*DIMENSION%neigd) )
+    ALLOCATE (eig_so(2*DIMENSION%neigd))
+    ALLOCATE (eigBuffer(2*DIMENSION%neigd,kpts%nkpt,wannierspin))
+    ALLOCATE (neigBuffer(kpts%nkpt,wannierspin))
+    results%eig = 1.0e300
+    eigBuffer = 1.0e300
+    results%neig = 0
+    neigBuffer = 0
     rsoc%soangl(:,:,:,:,:,:) = CONJG(rsoc%soangl(:,:,:,:,:,:))
     CALL timestop("eigenso: spnorb")
     !
     !--->    loop over k-points: each can be a separate task
-    !
-    !n_loc = INT(kpts%nkpt/mpi%isize)
-    !n_plus = kpts%nkpt - mpi%isize*n_loc
-    !i_plus = -1
-    !IF (mpi%irank.LT.n_plus) i_plus = 0
-    !n_end = (mpi%irank+1)+(n_loc+i_plus)*mpi%isize
-    !
-#if defined(CPP_MPI)
-     n_stride = kpts%nkpt/mpi%n_groups
-#else
-     n_stride = 1
-#endif
-     n_end = kpts%nkpt
-    !write(*,'(4i12)') mpi%irank, mpi%n_groups, n_stride, mpi%n_start
-    !
-    !--->  start loop k-pts
-    !
-    ! DO  nk = mpi%irank+1,n_end,mpi%isize
-     DO nk = mpi%n_start,n_end,n_stride
+    DO nk_i=1,SIZE(mpi%k_list)
+        nk=mpi%k_list(nk_i)
+     !DO nk = mpi%n_start,n_end,n_stride
        CALL lapw%init(input,noco, kpts,atoms,sym,nk,cell,.FALSE., mpi)
        ALLOCATE( zso(lapw%nv(1)+atoms%nlotot,2*DIMENSION%neigd,wannierspin))
        zso(:,:,:) = CMPLX(0.0,0.0)
@@ -154,31 +150,33 @@ CONTAINS
        IF (mpi%n_rank==0) THEN
           IF (input%eonly) THEN
              CALL write_eig(eig_id, nk,jspin,neig=nsz,neig_total=nsz, eig=eig_so(:nsz))
-          ELSE          
+             STOP 'jspin is undefined here (eigenso - eonly branch)'
+             eigBuffer(:nsz,nk,jspin) = eig_so(:nsz)
+             neigBuffer(nk,jspin) = nsz
+          ELSE
              CALL zmat%alloc(.FALSE.,SIZE(zso,1),nsz)
              DO jspin = 1,wannierspin
                 CALL timestart("eigenso: write_eig")  
                 zmat%data_c=zso(:,:nsz,jspin)
                 CALL write_eig(eig_id, nk,jspin,neig=nsz,neig_total=nsz, eig=eig_so(:nsz),zmat=zmat)
-   
-                CALL timestop("eigenso: write_eig")  
+                eigBuffer(:nsz,nk,jspin) = eig_so(:nsz)
+                neigBuffer(nk,jspin) = nsz
+                CALL timestop("eigenso: write_eig")
              ENDDO
           ENDIF ! (input%eonly) ELSE
        ENDIF ! n_rank == 0
        DEALLOCATE (zso)
     ENDDO ! DO nk 
 
-    ! Sorry for the following strange workaround to fill the results%neig and results%eig arrays.
-    ! At some point someone should have a closer look at how the eigenvalues are
-    ! distributed and fill the arrays without using the eigenvalue-IO.
-    DO jspin = 1, wannierspin
-       DO nk = 1,kpts%nkpt
-          CALL read_eig(eig_id,nk,jspin,results%neig(nk,jspin),results%eig(:,nk,jspin))
 #ifdef CPP_MPI
-          CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
+    CALL MPI_ALLREDUCE(neigBuffer,results%neig,kpts%nkpt*wannierspin,MPI_INTEGER,MPI_SUM,mpi%mpi_comm,ierr)
+    CALL MPI_ALLREDUCE(eigBuffer(:2*dimension%neigd,:,:),results%eig(:2*dimension%neigd,:,:),&
+                       2*dimension%neigd*kpts%nkpt*wannierspin,MPI_DOUBLE_PRECISION,MPI_MIN,mpi%mpi_comm,ierr)
+    CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
+#else
+    results%neig(:,:) = neigBuffer(:,:)
+    results%eig(:2*dimension%neigd,:,:) = eigBuffer(:2*dimension%neigd,:,:)
 #endif
-<