Initial commit for the integration of EELS as provided bei Roman Kovacik

I adapted Roman's code slightly such that it fits into the new fleur version.
Nothing is tested yet, so testing and polishing still needs to be done.

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

cmake/Files_and_Targets.txt

@@ -31,6 +31,7 @@ include(inpgen/CMakeLists.txt)
 include(docs/CMakeLists.txt)
 include(tests/CMakeLists.txt)
 include(mpi/CMakeLists.txt)
+include(eels/CMakeLists.txt)
 
 
 #include(wannier/CMakeLists.txt)

eels/CMakeLists.txt

+set(fleur_F90 ${fleur_F90}
+eels/corespec.f90
+eels/corespec_io.f90
+eels/corespec_eval.f90
+)

eels/corespec.f90

+!--------------------------------------------------------------------------------
+! Copyright (c) 2017 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.
+!--------------------------------------------------------------------------------
+
+module m_corespec
+
+  implicit none
+
+! 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
+
+  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/)
+  integer, parameter :: sign(4) = (/1,-1,0,0/)
+
+  character(len=2), parameter :: ssep="> "
+  character(len=3), parameter :: fsos1="t55"
+  character(len=3), parameter :: fsos2="t90"
+  character(len=15), parameter :: fsb="(a,('"//ssep//"'),a,"//fsos1
+  character(len=9), parameter :: fse=fsos2//",2x,a)"
+  character(len=11), parameter :: csmsgerr=" ... STOP !"
+  character(len=14), parameter :: csmsgwar=" ... WARNING !"
+  character(len=32), parameter :: csmsgs = fsb//")"
+  character(len=64), parameter :: csmsgsss = fsb//",a,a5,"//fse
+  character(len=64), parameter :: csmsgsis = fsb//",a,i5,"//fse
+  character(len=64), parameter :: csmsgsisis = fsb//",a,i5,a,i5,"//fse
+  character(len=64), parameter :: csmsgsfs = fsb//",a,f8.3,"//fse
+  character(len=64), parameter :: csmsgses = fsb//",a,es12.3,"//fse
+
+! INPUT PARAMETERS
+
+  type csitype
+     sequence
+     integer :: verb  ! output verbosity
+     integer :: type  ! atomic type used for calculation of core spectra
+     character(len=1) :: edge  ! edge character (K,L,M,N,O,P)
+     integer :: edgeidx(11)  ! l-j edges
+     integer :: lx  ! maximum lmax considered in spectra calculation
+     real :: ek0  ! kinetic energy of incoming electrons
+     real :: emn  ! energy spectrum lower bound
+     real :: emx  ! energy spectrum upper bound
+     real :: ein  ! energy spectrum increment
+  end type csitype
+
+! VARIABLES
+
+  logical :: l_cs
+
+  integer :: l1,l2,la1,la2,li
+  integer :: m1,m2,mu1,mu2,mi
+  integer :: lx,ln,lax,lan,lix,lin
+
+  character(len=32) :: smeno
+
+  type (csitype) :: csi
+
+  type csvtype
+     sequence
+     integer :: nc  ! main quantum no. of the core level of atomic type
+     integer :: nljc  ! number of l-j edge lines
+     integer, allocatable :: lc(:)  ! edge angular quantum nos.; nljc elements
+     real, allocatable :: eedge(:)  ! lc-dep. edge energy; nljc elements
+     real, allocatable :: occ(:)  ! lc-dep. occupation; nljc elements
+     integer :: nex  ! no. of energy sampling points
+     real, allocatable :: egrid(:)  ! energy grid; 0:nex elements
+     real, allocatable :: eos(:)  ! energy grid / sigma
+     real, allocatable :: eloss(:,:)  ! efermi-eedge+egrid
+     integer :: nen  ! minimum index for which egrid >=0
+     integer :: nqv  ! no. of q vectors
+     real :: qv0  ! |q| of incoming electrons
+     real, allocatable :: qv1(:,:,:)  ! |q| of outgoing electrons
+     real, allocatable :: qv(:,:,:,:)  ! delta q vectors
+     real :: gamma  ! gamma = 1+ek0/mc2
+     real :: beta  ! beta = v/c = 1/sqrt(1-1/gamma^2)
+     real, allocatable :: gaunt(:,:,:,:,:,:)  ! gaunt coefficients
+     real, allocatable :: fc(:,:,:,:)  ! core radial function
+     real, allocatable :: fv(:,:,:,:)  ! valence radial function
+     real, allocatable :: fb(:,:,:,:,:)  ! bessel function
+     real, allocatable :: rme(:,:,:,:,:,:,:,:)  ! matrix elements
+     real, allocatable :: dose(:,:,:,:,:)  ! dos (bands)
+     real, allocatable :: dosb(:,:,:,:,:)  ! dos (bands)
+     complex, allocatable :: ddscs(:,:,:,:,:)  ! dos (bands)
+     
+  end type csvtype
+
+  type (csvtype) :: csv
+
+end module m_corespec

eels/corespec_io.f90

+!--------------------------------------------------------------------------------
+! Copyright (c) 2017 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.
+!--------------------------------------------------------------------------------
+
+MODULE m_corespec_io
+
+  USE m_corespec
+  USE m_types
+  USE m_juDFT
+
+  IMPLICIT NONE
+
+  CONTAINS
+
+!===============================================================================
+!
+!  S U B R O U T I N E   C O R E S P E C _ I N I T
+!
+!-------------------------------------------------------------------------------
+!
+  SUBROUTINE corespec_init(atoms)
+
+    IMPLICIT NONE
+
+    TYPE(t_atoms),INTENT(IN)   :: atoms
+
+    INTEGER                    :: ui,i
+    LOGICAL                    :: lexist
+
+    namelist /csinp/ csi
+
+    smeno = "corespec_init"
+
+    write(*,'(/,a)') trim(smeno)//ssep
+
+    l_cs = .false.
+
+! initialization of input parameters: type csi
+
+    csi%verb = 1
+    csi%type = 0
+    csi%edge = ""
+    csi%edgeidx = 0
+    csi%lx = -1
+    csi%ek0 = 0.d0
+    csi%emn = -2.d0
+    csi%emx = 20.d0
+    csi%ein = 0.1d0
+
+! reading of input parameters from 'corespec_inp' file
+
+    call iounit(ui)
+    inquire(file="corespec_inp",exist=lexist)
+    if(lexist) then
+      open(ui,file="corespec_inp",status='old')
+      read(ui,nml=csinp)
+      close(ui)
+    else
+      return
+    endif
+
+    IF(ANY(atoms%nlo(:).NE.0)) CALL juDFT_error("EELS + LOs not available at the moment!" ,calledby ="corespec_io")
+
+! sanity check of the input parameters; if they are not correct, program stops
+! unit conversion if necessary
+! if csi%verb = 1, detailed information is written to stdout
+
+! csi%type
+    if(csi%type.le.0) then
+      write(*,csmsgs)  trim(smeno),"found csi%type <= 0 !"//csmsgerr ; stop
+    endif
+    if(csi%type.gt.atoms%ntype) then
+      write(*,csmsgs)  trim(smeno),"found csi%type > atoms%ntype!"//csmsgerr ; stop
+    endif
+    if(csi%verb.eq.1) write(*,csmsgsis)  trim(smeno),&
+           &"atomic type: ","csi%type = ",csi%type,"will be used"
+
+! csi%edge -> csv%nc
+    csv%nc = 0
+    if(len(trim(csi%edge)).ne.0) csv%nc = ichar(csi%edge)-ichar('K')+1
+    if(csi%edge.eq."") then
+      write(*,csmsgs) trim(smeno),"found empty csi%edge !"//csmsgerr ; stop
+    endif
+    if(csv%nc.lt.1.or.csv%nc.gt.6) then
+      write(*,csmsgs) trim(smeno),&
+           &"specify csi%edge as one of {K,L,M,N,O,P} !"//csmsgerr ; stop
+    endif
+    if(csi%verb.eq.1) write(*,csmsgsss)  trim(smeno),&
+           &"edge: ","csi%edge = ",csi%edge,"will be used"
+    if(csi%verb.eq.1) write(*,csmsgsis)  trim(smeno),&
+           &"main quantum no.: ","csv%nc = ",csv%nc,"will be used"
+
+! csi%edgeidx(:) -> csv%nljc
+    if(maxval(csi%edgeidx).gt.(2*csv%nc-1)) then
+      write(*,csmsgs) trim(smeno),&
+         &"found csi%edgeidx > 2*csv%nc-1 !"//csmsgerr ; stop
+    endif
+    if(count(csi%edgeidx.gt.0).gt.(2*csv%nc-1)) then
+      write(*,csmsgs) trim(smeno),&
+         &"found more than 2*csv%nc-1 of csi%edgeidx > 0 !"//csmsgerr ; stop
+    endif
+    if((csv%nc-1)**2+maxval(csi%edgeidx).gt.atoms%ncst(csi%type)) then
+      write(*,csmsgs) trim(smeno),&
+         &"found (csv%nc-1)^2+maxval(csi%edgeidx) > atoms%ncst(csi%type)!"//csmsgerr
+      stop
+    endif
+    csv%nljc = count(csi%edgeidx.gt.0)
+    if(.not.allocated(csv%lc)) allocate(csv%lc(csv%nljc))
+    csv%lc = (/(edgel(csi%edgeidx(i)), i = 1,csv%nljc)/)
+!!$    print*,csv%lc
+    if(csi%verb.eq.1) write(*,csmsgsis)  trim(smeno),&
+           &"nljc edge lines: ","csv%nljc = ",csv%nljc,"will be used"
+
+! csi%lx
+    if(csi%lx.lt.0) then
+      write(*,csmsgs)  trim(smeno),"found csi%lx < 0 !"//csmsgerr ; stop
+    endif
+    if(csi%lx.gt.atoms%lmax(csi%type)) then
+      write(*,csmsgs)  trim(smeno),&
+           &"found csi%lx > atoms%lmax(csi%type)!"//csmsgerr ; stop
+    endif
+    if(csi%verb.eq.1) write(*,csmsgsis)  trim(smeno),&
+           &"maximum l: ","csi%lx = ",csi%lx,"will be used"
+
+! csi%ek0
+    if(csi%ek0.le.0.d0) 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))
+    if(csi%verb.eq.1) then
+      write(*,csmsgses)  trim(smeno),&
+           &"kinetic energy of incoming electrons: ","csi%ek0 = ",csi%ek0,&
+           &"eV will be used"
+      write(*,csmsgses)  trim(smeno),&
+           &"Lorentz factor (gamma): ","csv%gamma = ",csv%gamma,&
+           &"will be used"
+      write(*,csmsgses)  trim(smeno),&
+           &"v/c factor (beta): ","csv%beta = ",csv%beta,&
+           &"will be used"
+    endif
+
+! csi%emn csi%emx csi%ein -> csv%nex csv%egrid(0:csv%nex)
+    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
+    endif
+    if(((csi%emx-csi%emn)/csi%ein)-int((csi%emx-csi%emn)/csi%ein).ne.0) then
+      write(*,csmsgs)  trim(smeno),&
+           &"found non-integer (csi%emx-csi%emn)/csi%ein !"//csmsgerr ; stop
+    endif
+    csv%nex = int((csi%emx-csi%emn)/csi%ein)
+    if(.not.allocated(csv%egrid)) allocate(csv%egrid(0:csv%nex))
+    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
+!!$        csv%nen = i
+!!$        exit
+!!$      endif
+!!$    enddo
+
+!!$    print*,csv%egrid
+!!$    print*,csv%nen
+    if(csi%verb.eq.1) write(*,csmsgsfs)  trim(smeno),&
+           &"energy spectrum lower bound: ","csi%emn = ",csi%emn,&
+           &"eV will be used"
+    if(csi%verb.eq.1) write(*,csmsgsfs)  trim(smeno),&
+           &"energy spectrum upper bound: ","csi%emx = ",csi%emx,&
+           &"eV will be used"
+    if(csi%verb.eq.1) write(*,csmsgsfs)  trim(smeno),&
+           &"energy spectrum increment: ","csi%ein = ",csi%ein,"eV will be used"
+    if(csi%verb.eq.1) write(*,csmsgsis)  trim(smeno),&
+           &"no. of energy spectrum grid points: ","csv%nex = 0 : ",&
+           &csv%nex,"will be used"
+    if(csi%verb.eq.1) write(*,csmsgsis)  trim(smeno),&
+           &"minimum index for which egrid >= 0: ","csv%nen = ",&
+           &csv%nen,"will be used"
+
+    if(.not.allocated(csv%eedge)) allocate(csv%eedge(csv%nljc))
+    csv%eedge = 0.d0
+    if(.not.allocated(csv%occ)) allocate(csv%occ(csv%nljc))
+    csv%occ = 0.d0
+
+    l_cs = .true.
+
+    if(csi%verb.eq.1) write(*,*) ""
+
+  end subroutine corespec_init
+!
+!===============================================================================
+!===============================================================================
+!
+!  S U B R O U T I N E   I O U N I T
+!
+!-------------------------------------------------------------------------------
+!
+  subroutine iounit(unit)
+!
+! Assignes unused integer to I/O unit
+!
+    implicit none
+!
+    integer, intent(out) :: unit
+!
+    integer :: i
+    logical :: fopen
+!
+    i = 9
+    fopen = .true.
+    do while(fopen)
+      i = i+1
+      inquire(unit=i,opened=fopen)
+    enddo
+!
+    unit = i
+!
+  end subroutine iounit
+!
+!===============================================================================
+
+end module m_corespec_io