Finished merge

cdn/cdnval.F90

@@ -84,6 +84,9 @@ CONTAINS
     USE m_doswrite
     USE m_cylpts
     USE m_cdnread, ONLY : cdn_read0, cdn_read
+    USE m_corespec, only : l_cs    ! calculation of core spectra (EELS)
+    USE m_corespec_io, only : corespec_init
+    USE m_corespec_eval, only : corespec_gaunt,corespec_rme,corespec_dos,corespec_ddscs
 #ifdef CPP_MPI
     USE m_mpi_col_den ! collect density data from parallel nodes
 #endif
@@ -324,6 +327,11 @@ CONTAINS
        ALLOCATE ( m_mcd(1,1,1,1),mcd(1,1,1) )
     ENDIF
 
+! calculation of core spectra (EELS) initializations -start-
+    CALL corespec_init(atoms)
+    IF(l_cs.AND.jspin.EQ.1) CALL corespec_gaunt()
+! calculation of core spectra (EELS) initializations -end-
+
     ALLOCATE ( kveclo(atoms%nlotot) )
 
     IF (mpi%irank==0) THEN
@@ -361,7 +369,7 @@ CONTAINS
     ncored = 0
 
     ALLOCATE ( flo(atoms%jmtd,2,atoms%nlod,dimension%jspd) )
-    DO  n = 1,atoms%ntype
+    DO n = 1,atoms%ntype
        IF (input%cdinf.AND.mpi%irank==0) WRITE (6,FMT=8001) n
        DO  l = 0,atoms%lmax(n)
           DO ispin = jsp_start,jsp_end
@@ -388,6 +396,11 @@ CONTAINS
                ncore,e_mcd,m_mcd)
           ncored = max(ncore(n),ncored)
        END IF
+
+       IF(l_cs) CALL corespec_rme(atoms,input,n,dimension%nstd,&
+                                  input%jspins,jspin,results%ef,&
+                                  dimension%msh,vr(:,0,:,:),f,g)
+
        !
        !--->   generate the extra wavefunctions for the local orbitals,
        !--->   if there are any.
@@ -802,6 +815,13 @@ CONTAINS
                 DEALLOCATE (acoflo,bcoflo,cveccof)
                 CALL timestop("cdnval: force_a12/21")
              END IF
+
+             IF(l_cs) THEN
+                CALL corespec_dos(atoms,usdus,ispin,dimension%lmd,kpts%nkpt,ikpt,&
+                                  dimension%neigd,noccbd,results%ef,banddos%sig_dos,&
+                                  eig,we,acof(1,0,1,ispin),bcof(1,0,1,ispin),&
+                                  ccof(-atoms%llod,1,1,1,ispin))
+             END IF
           END DO !--->    end loop over ispin
 
           IF (noco%l_mperp) THEN
@@ -872,6 +892,9 @@ CONTAINS
     END DO
     CALL timestop("cdnval: mpi_col_den")
 #endif
+
+    IF(l_cs) CALL corespec_ddscs(jspin,input%jspins)
+
     IF (((jspin.eq.input%jspins).OR.noco%l_mperp) .AND. (banddos%dos.or.banddos%vacdos.or.input%cdinf) ) THEN
        CALL timestart("cdnval: dos")
        IF (mpi%irank==0) THEN

cdn/n_mat.f90

@@ -11,6 +11,7 @@ MODULE m_nmat
   !     all atoms are stored in lda_u(), if lda_u()<0, no +U is used.
   !     For details see Eq.(12) of Shick et al. PRB 60, 10765 (1999)
   !     Part of the LDA+U package                   G.B., Oct. 2000
+  !     Extension to multiple U per atom type by G.M. 2017
   !     ************************************************************
 CONTAINS
   SUBROUTINE n_mat(atoms,sym, ne,usdus,jspin,we, acof,bcof,ccof, n_mmp)
@@ -34,7 +35,7 @@ CONTAINS
     !     ..
     !     .. Local Scalars ..
     COMPLEX c_0
-    INTEGER i,j,k,l ,mp,n,it,is,isi,natom,n_ldau,lp,m
+    INTEGER i,j,k,l ,mp,n,it,is,isi,natom,natomTemp,n_ldau,lp,m,i_u
     INTEGER ilo,ilop,ll1,nn,lmp,lm
     REAL fac
     !     ..
@@ -45,16 +46,17 @@ CONTAINS
     !
     ! calculate n_mat:
     !
-    n_ldau = 0
     natom = 0
+    i_u = 1
     DO n = 1,atoms%ntype
-       IF (atoms%lda_u(n)%l.GE.0) THEN
-          n_ldau = n_ldau + 1
-          n_tmp(:,:) =cmplx(0.0,0.0)
-          l = atoms%lda_u(n)%l
+       DO WHILE (i_u.LE.atoms%n_u)
+          IF (atoms%lda_u(i_u)%atomType.GT.n) EXIT
+          natomTemp = natom
+          n_tmp(:,:) = cmplx(0.0,0.0)
+          l = atoms%lda_u(i_u)%l
           ll1 = (l+1)*l 
           DO nn = 1, atoms%neq(n)
-             natom = natom + 1
+             natomTemp = natomTemp + 1
              !
              !  prepare n_mat in local frame (in noco-calculations this depends 
              !                                also on alpha(n) and beta(n) )
@@ -66,8 +68,8 @@ CONTAINS
                    c_0 = cmplx(0.0,0.0)
                    DO i = 1,ne
                       c_0 = c_0 +  we(i) * ( usdus%ddn(l,n,jspin) *&
-                           conjg(bcof(i,lmp,natom))*bcof(i,lm,natom) +&
-                           conjg(acof(i,lmp,natom))*acof(i,lm,natom) )
+                           conjg(bcof(i,lmp,natomTemp))*bcof(i,lm,natomTemp) +&
+                           conjg(acof(i,lmp,natomTemp))*acof(i,lm,natomTemp) )
                    ENDDO
                    n_tmp(m,mp) = c_0 
                 ENDDO
@@ -85,17 +87,17 @@ CONTAINS
                          c_0 = cmplx(0.0,0.0)
                          DO i = 1,ne
                             c_0 = c_0 +  we(i) * (  usdus%uulon(ilo,n,jspin) * (&
-                                 conjg(acof(i,lmp,natom))*ccof(m,i,ilo,natom) +&
-                                 conjg(ccof(mp,i,ilo,natom))*acof(i,lm,natom) )&
+                                 conjg(acof(i,lmp,natomTemp))*ccof(m,i,ilo,natomTemp) +&
+                                 conjg(ccof(mp,i,ilo,natomTemp))*acof(i,lm,natomTemp) )&
                                  + usdus%dulon(ilo,n,jspin) * (&
-                                 conjg(bcof(i,lmp,natom))*ccof(m,i,ilo,natom) +&
-                                 conjg(ccof(mp,i,ilo,natom))*bcof(i,lm,natom)))
+                                 conjg(bcof(i,lmp,natomTemp))*ccof(m,i,ilo,natomTemp) +&
+                                 conjg(ccof(mp,i,ilo,natomTemp))*bcof(i,lm,natomTemp)))
                          ENDDO
                          DO ilop = 1, atoms%nlo(n)
                             IF (atoms%llo(ilop,n).EQ.l) THEN
                                DO i = 1,ne
                                   c_0 = c_0 +  we(i) * usdus%uloulopn(ilo,ilop,n,jspin) *&
-                                       conjg(ccof(mp,i,ilop,natom)) *ccof(m ,i,ilo ,natom)
+                                       conjg(ccof(mp,i,ilop,natomTemp)) *ccof(m ,i,ilo ,natomTemp)
                                ENDDO
                             ENDIF
                          ENDDO
@@ -108,10 +110,10 @@ CONTAINS
              !
              !  n_mmp should be rotated by D_mm' ; compare force_a21
              !
-             DO it = 1, sym%invarind(natom)
+             DO it = 1, sym%invarind(natomTemp)
 
-                fac = 1.0  /  ( sym%invarind(natom) * atoms%neq(n) )
-                is = sym%invarop(natom,it)
+                fac = 1.0  /  ( sym%invarind(natomTemp) * atoms%neq(n) )
+                is = sym%invarop(natomTemp,it)
                 isi = sym%invtab(is)
                 d_tmp(:,:) = cmplx(0.0,0.0)
                 DO m = -l,l
@@ -123,16 +125,16 @@ CONTAINS
                 n1_tmp =  matmul( nr_tmp, d_tmp )
                 DO m = -l,l
                    DO mp = -l,l
-                      n_mmp(m,mp,n_ldau) = n_mmp(m,mp,n_ldau) +conjg(n1_tmp(m,mp)) * fac
+                      n_mmp(m,mp,i_u) = n_mmp(m,mp,i_u) + conjg(n1_tmp(m,mp)) * fac
                    ENDDO
                 ENDDO
 
              ENDDO
 
           ENDDO ! sum  over equivalent atoms
-       ELSE
-          natom = natom + atoms%neq(n)
-       ENDIF
+          i_u = i_u + 1
+       END DO
+       natom = natom + atoms%neq(n)
     ENDDO     ! loop over atom types
 
     !     do m=-l,l

cdn_mt/cdnmt.f90

-!--------------------------------------------------------------------------------
 ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
 ! This file is part of FLEUR and available as free software under the conditions
 ! of the MIT license as expressed in the LICENSE file in more detail.
@@ -67,12 +66,12 @@ CONTAINS
     !     ..
     !     .. Local Scalars ..
     INTEGER itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu
-    INTEGER ilo,ilop
+    INTEGER ilo,ilop,i
     REAL s,wronk,sumlm,qmtt
     COMPLEX cs
     !     ..
     !     .. Local Arrays ..
-    REAL qmtl(0:atoms%lmaxd),qmtllo(0:atoms%lmaxd)
+    REAL qmtl(0:atoms%lmaxd,jspd,atoms%ntype),qmtllo(0:atoms%lmaxd)
     CHARACTER(LEN=20) :: attributes(6)
 
     !     ..
@@ -80,8 +79,10 @@ CONTAINS
     REAL,    ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:)
     COMPLEX, ALLOCATABLE :: rho21(:,:,:)
     !
+
+    CALL timestart("cdnmt")
+
     IF (noco%l_mperp) THEN
-       ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jspd),g(atoms%jmtd,2,0:atoms%lmaxd,jspd) )
        ALLOCATE ( usdus%us(0:atoms%lmaxd,atoms%ntype,jspd),usdus%uds(0:atoms%lmaxd,atoms%ntype,jspd) )
        ALLOCATE ( usdus%dus(0:atoms%lmaxd,atoms%ntype,jspd),usdus%duds(0:atoms%lmaxd,atoms%ntype,jspd) )
        ALLOCATE ( usdus%ddn(0:atoms%lmaxd,atoms%ntype,jspd) )
@@ -90,22 +91,34 @@ CONTAINS
           rho21(:,:,:) = cmplx(0.0,0.0)
        ENDIF
     ELSE
-       ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
-       ALLOCATE ( g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
        ALLOCATE (   usdus%us(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) )
        ALLOCATE (  usdus%uds(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) )
        ALLOCATE (  usdus%dus(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) )
        ALLOCATE ( usdus%duds(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) )
        ALLOCATE (  usdus%ddn(0:atoms%lmaxd,atoms%ntype,jsp_start:jsp_end) )
     ENDIF
-    WRITE (6,FMT=8000)
-    WRITE (16,FMT=8000)
-8000 FORMAT (/,5x,'l-like charge',/,t6,'atom',t15,'s',t24,'p',&
-         &     t33,'d',t42,'f',t51,'total')
-    CALL timestart("cdnmt")
 
-    na = 1
+    !$OMP PARALLEL DEFAULT(none) &
+    !$OMP SHARED(usdus,rho,chmom,clmom,qa21,rho21) &
+    !$OMP SHARED(atoms,jsp_start,jsp_end,epar,vr,uu,dd,du,sphhar,uloulopn,ello,aclo,bclo,cclo) &
+    !$OMP SHARED(acnmt,bcnmt,ccnmt,orb,orbl,orblo,ddnmt,udnmt,dunmt,uunmt,mt21,lo21,uloulop21)&
+    !$OMP SHARED(uloulopn21,noco,l_fmpl,uunmt21,ddnmt21,dunmt21,udnmt21,jspd)&
+    !$OMP PRIVATE(itype,na,ispin,l,f,g,nodeu,noded,wronk,i,j,s,qmtllo,qmtt,qmtl,nd,lh,lp,llp,cs)
+    IF (noco%l_mperp) THEN
+       ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jspd),g(atoms%jmtd,2,0:atoms%lmaxd,jspd) )
+    ELSE
+       ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
+       ALLOCATE ( g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
+    ENDIF
+
+    qmtl = 0
+    
+    !$OMP DO
     DO itype = 1,atoms%ntype
+       na = 1
+       DO i = 1, itype - 1
+          na = na + atoms%neq(i)
+       ENDDO
        !--->    spherical component
        DO ispin = jsp_start,jsp_end
           DO l = 0,atoms%lmax(itype)
@@ -126,7 +139,6 @@ CONTAINS
              qmtllo(l) = 0.0
           END DO
 
-
           CALL rhosphnlo(itype,atoms,sphhar,&
                uloulopn(1,1,itype,ispin),usdus%dulon(1,itype,ispin),&
                usdus%uulon(1,itype,ispin),ello(1,itype,ispin),&
@@ -138,26 +150,13 @@ CONTAINS
 
 
           !--->       l-decomposed density for each atom type
-          qmtt = 0.
+          qmtt = 0.0
           DO l = 0,atoms%lmax(itype)
-             qmtl(l) = ( uu(l,itype,ispin)+dd(l,itype,ispin)&
+             qmtl(l,ispin,itype) = ( uu(l,itype,ispin)+dd(l,itype,ispin)&
                   &              *usdus%ddn(l,itype,ispin) )/atoms%neq(itype) + qmtllo(l)
-             qmtt = qmtt + qmtl(l)
+             qmtt = qmtt + qmtl(l,ispin,itype)
           END DO
           chmom(itype,ispin) = qmtt
-          WRITE (6,FMT=8100) itype, (qmtl(l),l=0,3),qmtt
-          WRITE (16,FMT=8100) itype, (qmtl(l),l=0,3),qmtt
-8100      FORMAT (' -->',i3,2x,4f9.5,2x,f9.5)
-
-          attributes = ''
-          WRITE(attributes(1),'(i0)') itype
-          WRITE(attributes(2),'(f12.7)') qmtt
-          WRITE(attributes(3),'(f12.7)') qmtl(0)
-          WRITE(attributes(4),'(f12.7)') qmtl(1)
-          WRITE(attributes(5),'(f12.7)') qmtl(2)
-          WRITE(attributes(6),'(f12.7)') qmtl(3)
-          CALL writeXMLElementForm('mtCharge',(/'atomType','total   ','s       ','p       ','d       ','f       '/),attributes,&
-                                   reshape((/8,5,1,1,1,1,6,12,12,12,12,12/),(/6,2/)))
 
           !+soc
           !--->       spherical angular component
@@ -254,9 +253,35 @@ CONTAINS
           ENDIF ! l_fmpl
        ENDIF ! noco%l_mperp
 
-       na = na + atoms%neq(itype)
     ENDDO ! end of loop over atom types
+    !$OMP END DO
+    DEALLOCATE ( f,g)
+    !$OMP END PARALLEL
+
+    WRITE (6,FMT=8000)
+    WRITE (16,FMT=8000)
+8000 FORMAT (/,5x,'l-like charge',/,t6,'atom',t15,'s',t24,'p',&
+         &     t33,'d',t42,'f',t51,'total')
+
+    DO itype = 1,atoms%ntype
+       DO ispin = jsp_start,jsp_end
+          WRITE ( 6,FMT=8100) itype, (qmtl(l,ispin,itype),l=0,3),chmom(itype,ispin)
+          WRITE (16,FMT=8100) itype, (qmtl(l,ispin,itype),l=0,3),chmom(itype,ispin)
+8100      FORMAT (' -->',i3,2x,4f9.5,2x,f9.5)
+          attributes = ''
+          WRITE(attributes(1),'(i0)') itype
+          WRITE(attributes(2),'(f12.7)') chmom(itype,ispin)
+          WRITE(attributes(3),'(f12.7)') qmtl(0,ispin,itype)
+          WRITE(attributes(4),'(f12.7)') qmtl(1,ispin,itype)
+          WRITE(attributes(5),'(f12.7)') qmtl(2,ispin,itype)
+          WRITE(attributes(6),'(f12.7)') qmtl(3,ispin,itype)
+          CALL writeXMLElementForm('mtCharge',(/'atomType','total   ','s       ','p       ','d       ','f       '/),attributes,&
+                                   reshape((/8,5,1,1,1,1,6,12,12,12,12,12/),(/6,2/)))
+       ENDDO
+    ENDDO
+
     CALL timestop("cdnmt")
+
     !---> for testing: to plot the offdiag. part of the density matrix it
     !---> is written to the file rhomt21. This file can read in pldngen.
     IF (l_fmpl) THEN
@@ -267,7 +292,5 @@ CONTAINS
     ENDIF
     !---> end of test output
 
-    DEALLOCATE ( f,g)
-
   END SUBROUTINE cdnmt
 END MODULE m_cdnmt

cmake/Files_and_Targets.txt

@@ -32,8 +32,10 @@ include(docs/CMakeLists.txt)
 include(tests/CMakeLists.txt)
 include(mpi/CMakeLists.txt)
 include(hybrid/CMakeLists.txt)
+include(eels/CMakeLists.txt)
 
 #include(wannier/CMakeLists.txt)
+#include(wannier/uhu/CMakeLists.txt)
 
 set(inpgen_F77 
 global/constants.f inpgen/element.f inpgen/atom_input.f inpgen/crystal.f inpgen/lattice2.f inpgen/setab.f inpgen/super_check.f

cmake/tests/test_ELPA.cmake

@@ -30,11 +30,14 @@ LINK_LIBRARIES ${FLEUR_LIBRARIES})
 try_compile(FLEUR_USE_ELPA_NEW ${CMAKE_BINARY_DIR} ${CMAKE_SOURCE_DIR}/cmake/tests/test_ELPA_NEW.f90
 LINK_LIBRARIES ${FLEUR_LIBRARIES})
 try_compile(FLEUR_USE_ELPA_201605003 ${CMAKE_BINARY_DIR} ${CMAKE_SOURCE_DIR}/cmake/tests/test_ELPA_201605003.f90
+LINK_LIBRARIES ${FLEUR_LIBRARIES})
+try_compile(FLEUR_USE_ELPA_201605004 ${CMAKE_BINARY_DIR} ${CMAKE_SOURCE_DIR}/cmake/tests/test_ELPA_201605004.f90
 LINK_LIBRARIES ${FLEUR_LIBRARIES})
     message("Version check for ELPA:")
     message("OLD ELPA      : ${FLEUR_USE_ELPA_OLD}")
     message("NEW ELPA      : ${FLEUR_USE_ELPA_NEW}")
     message("201605003 ELPA: ${FLEUR_USE_ELPA_201605003}")
+    message("201605004 ELPA: ${FLEUR_USE_ELPA_201605004}")
 #Set preprocessor switches
    if (FLEUR_USE_ELPA_OLD)
        set(FLEUR_USE_ELPA TRUE)
@@ -48,4 +51,8 @@ LINK_LIBRARIES ${FLEUR_LIBRARIES})
        set(FLEUR_USE_ELPA TRUE)
        set(FLEUR_MPI_DEFINITIONS ${FLEUR_MPI_DEFINITIONS} "CPP_ELPA" "CPP_ELPA2" "CPP_ELPA_201605003")
    endif()
+   if (FLEUR_USE_ELPA_201605004)
+       set(FLEUR_USE_ELPA TRUE)
+       set(FLEUR_MPI_DEFINITIONS ${FLEUR_MPI_DEFINITIONS} "CPP_ELPA" "CPP_ELPA2" "CPP_ELPA_201605004")
+   endif()
 endif() 

cmake/tests/test_ELPA_201605004.f90

+      program test
+      use elpa1
+      integer:: ierr,mpi_subcom, myrowblacs, mycolblacs
+      integer:: mpi_comm_rows,mpi_comm_cols,m,nb,myrowsssca,mycolssca
+      logical :: ok
+      real :: bsca(10,10),asca(10,10),eigvec(10)
+      ok= elpa_mult_at_b_real('U', 'L',m, m,bsca,myrowssca,mycolssca,asca,SIZE(asca,1),SIZE(asca,2),nb,&
+         mpi_comm_rows, mpi_comm_cols,eigvec,myrowssca,mycolssca)
+      !ok=CHOLESKY_real (m,bsca,SIZE(bsca,1),nb,mycolssca,mpi_comm_rows,mpi_comm_cols,.false.)      
+      end

cmake/tests/test_XML.cmake

@@ -21,7 +21,7 @@ message("XML Library found for linking:${FLEUR_USE_XML}")
 
 if (FLEUR_USE_XML)
    try_compile(FLEUR_USE_XML ${CMAKE_BINARY_DIR} ${CMAKE_SOURCE_DIR}/cmake/tests/test_XML.c
-   LINK_LIBRARIES "-lxml2")
+   CMAKE_FLAGS "-DCMAKE_C_LINK_EXECUTABLE='echo no linking'" LINK_LIBRARIES "-lxml2")
    if (NOT FLEUR_USE_XML)
       find_package(LibXml2)
       set(CMAKE_C_FLAGS "-I${LIBXML2_INCLUDE_DIR}")

core/cored.F90

@@ -224,7 +224,7 @@ CONTAINS
 8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f9.6,5x,&
          &       'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x,&
          &       'weight')
-8010 FORMAT (12x,2f5.0,f6.1,f10.4,f10.0)
+8010 FORMAT (12x,2f5.0,f6.1,f10.4,f12.4)
 8020 FORMAT (f20.8,'  electrons lost from core.')
 8030 FORMAT (10x,'atom type',i3,'  (spin',i2,') ',/,10x,&
          &       'kinetic energy=',e20.12,5x,'sum of the eigenvalues=',&

core/etabinit.F90

@@ -127,7 +127,7 @@ CONTAINS
     ENDDO
 8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f8.6,5x,&
                 'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x, 'weight')
-8010 FORMAT (12x,2f5.0,f6.1,f10.4,f10.0)
+8010 FORMAT (12x,2f5.0,f6.1,f10.4,f12.4)
 8020 FORMAT (/,/,12x,'core e.v. initialization')
 
   END SUBROUTINE etabinit

core/setcor.f90

@@ -23,7 +23,7 @@ CONTAINS
     REAL,INTENT (OUT)    :: occ(:,:)
     !     ..
     !     .. Local Scalars ..
-    INTEGER iz,jz,jz0,k,n,m,i,jspin,tempInt
+    INTEGER iz,jz,jz0,k,n,m,i,tempInt
     INTEGER k_h(2),n_h(2)
     REAL fj,l,bmu_l,o_h(2), fac(2),tempReal
     LOGICAL l_clf
@@ -53,8 +53,6 @@ CONTAINS
        ENDDO
        CLOSE (61)
        RETURN
-    ELSE
-       jspin=1
     ENDIF
 
     IF (atoms%zatom(itype)>92.01e0)  CALL juDFT_error(" z > 92",calledby ="setcor"&
@@ -214,6 +212,9 @@ CONTAINS
     ! modify default electron configuration according to explicitely provided setting in inp.xml