Line data    Source code

       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2020 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : MODULE m_types_mcd
       7             :   use m_judft
       8             :   use m_types_eigdos
       9             :   implicit none
      10             :   PRIVATE
      11             :   public t_mcd
      12             :   TYPE,extends(t_eigdos):: t_mcd
      13             :     REAL                 :: emcd_lo, emcd_up, maxE_mcd
      14             :     INTEGER, ALLOCATABLE :: ncore(:)
      15             :     INTEGER,ALLOCATABLE  :: n_dos_to_type(:)
      16             :     REAL,    ALLOCATABLE :: e_mcd(:,:,:)
      17             :     REAL,    ALLOCATABLE :: mcd(:,:,:,:,:)
      18             :     COMPLEX, ALLOCATABLE :: m_mcd(:,:,:,:)
      19             :     REAL,    ALLOCATABLE :: mcd_grid(:)
      20             :   CONTAINS
      21             :     PROCEDURE,PASS :: init => mcd_init
      22             :     procedure      :: get_dos_grid
      23             :     PROCEDURE      :: make_dos
      24             :     PROCEDURE      :: get_weight_eig
      25             :     PROCEDURE      :: get_num_weights
      26             :     PROCEDURE      :: get_weight_name
      27             :     END TYPE t_mcd
      28             : contains
      29             : 
      30           0 : function get_dos_grid(this)
      31             :   class(t_mcd),intent(in):: this
      32             :   real,allocatable:: get_dos_grid(:)
      33             : 
      34             :   INTEGER :: ind,ntype,nc
      35             :   REAL:: e_core
      36           0 :   get_dos_grid=this%mcd_grid
      37             : end function
      38             : 
      39           0 : subroutine make_dos(eigdos,kpts,input,banddos,efermi)
      40             :     use m_types_banddos
      41             :     use m_types_input
      42             :     use m_types_kpts
      43             : 
      44             :     class(t_mcd),intent(inout)   :: eigdos
      45             :     type(t_banddos),intent(in)   :: banddos
      46             :     type(t_input),intent(in)     :: input
      47             :     type(t_kpts),intent(in)      :: kpts
      48             :     real,intent(in)              :: efermi
      49             : 
      50             :     integer ::n,i,ind,ntype,nc,k,l,jspin
      51             :     real    :: e1,e2,e_lo,e_up,fac
      52           0 :     real,allocatable:: dos(:,:,:)
      53           0 :     if (allocated(eigdos%dos)) return
      54             :     !Call the routine of the parent-class
      55           0 :     call t_eigdos_make_dos(eigdos,kpts,input,banddos,efermi)
      56             : 
      57             :     !Only unoccupied states
      58           0 :     DO n=1,size(eigdos%dos_grid)
      59           0 :       if (eigdos%dos_grid(n)<0.0)  eigdos%dos(n,:,:)=0.0
      60             :     enddo
      61             : 
      62             :     !Map the values to MCD grid
      63             : 
      64           0 :     e_lo =  minval(eigdos%e_mcd)-efermi-maxval(eigdos%dos_grid) - 3.0*banddos%sig_dos
      65           0 :     e_up =  eigdos%maxE_mcd-efermi + 3.0*banddos%sig_dos
      66           0 :     ALLOCATE(eigdos%mcd_grid(size(eigdos%dos_grid)))
      67           0 :     DO i=1,size(eigdos%dos_grid)
      68           0 :       eigdos%mcd_grid(i)=e_lo+(i-1)*(e_up-e_lo)/(size(eigdos%mcd_grid)-1)
      69             :     ENDDO
      70             : 
      71           0 :     allocate(dos,mold=eigdos%dos)
      72           0 :     dos=0.0
      73           0 :     ind=0
      74           0 :     DO ntype=1,size(eigdos%ncore)
      75           0 :       DO nc=1,eigdos%ncore(ntype)
      76           0 :         DO k = 1,3
      77           0 :           ind=ind+1
      78           0 :           DO jspin=1,size(eigdos%e_mcd,2)
      79           0 :             DO i=1,size(eigdos%dos_grid)-1
      80           0 :               e1=-1*eigdos%dos_grid(i)-efermi+eigdos%e_mcd(ntype,jspin,nc)
      81           0 :               e2=-1*eigdos%dos_grid(i+1)-efermi+eigdos%e_mcd(ntype,jspin,nc)
      82           0 :               DO l=1,size(eigdos%mcd_grid)
      83           0 :                 IF ((e2.LE.eigdos%mcd_grid(l)).AND. (e1.GT.eigdos%mcd_grid(l))) THEN
      84           0 :                   fac = (eigdos%mcd_grid(l)-e1)/(e2-e1)
      85           0 :                   dos(l,jspin,ind) = dos(l,jspin,ind)+ eigdos%dos(i,jspin,ind)*(1.-fac) + fac * eigdos%dos(i+1,jspin,ind)
      86             :                 ENDIF
      87             :               ENDDO
      88             :             ENDDO
      89             :           ENDDO
      90             :         ENDDO
      91             :       ENDDO
      92             :     ENDDO
      93           0 :     eigdos%dos=dos
      94           0 :   end subroutine
      95             : 
      96           0 : function get_weight_eig(this,id)
      97             :   class(t_mcd),intent(in):: this
      98             :   INTEGER,intent(in)         :: id
      99             :   real,allocatable:: get_weight_eig(:,:,:)
     100             : 
     101             :   INTEGER :: ind,ntype,nc
     102             : 
     103           0 :   ind=0
     104           0 :   DO ntype=1,size(this%ncore)
     105           0 :     DO nc=1,this%ncore(ntype)
     106           0 :       ind=ind+1
     107           0 :       if (ind==id) get_weight_eig=this%mcd(ntype,nc,:,:,:)
     108           0 :       ind=ind+1
     109           0 :       if (ind==id) get_weight_eig=this%mcd(ntype+1,nc,:,:,:)
     110           0 :       ind=ind+1
     111           0 :       if (ind==id) get_weight_eig=this%mcd(ntype+2,nc,:,:,:)
     112           0 :       IF(ind>id) return
     113             :     ENDDO
     114             :   ENDDO
     115           0 :   IF(ind>id)CALL judft_error("Types_mcd: data not found")
     116             : 
     117             : END function
     118             : 
     119           0 : integer function get_num_weights(this)
     120             :   class(t_mcd),intent(in):: this
     121           0 :   get_num_weights=3*sum(this%ncore)
     122           0 : end function
     123             : 
     124           0 :   character(len=20) function get_weight_name(this,id)
     125             :     class(t_mcd),intent(in):: this
     126             :     INTEGER,intent(in)         :: id
     127             : 
     128             :     character(len=3):: c
     129             :     INTEGER :: ind,n_dos,nc,n
     130           0 :     ind=0
     131           0 :     DO n=1,size(this%mcd,1)
     132           0 :       n_dos=(n-1)/3+1
     133           0 :       select case(mod(n,3))
     134             :       case(1)
     135           0 :         c="pos"
     136             :       case(2)
     137           0 :         c="neg"
     138             :       case(0)
     139           0 :         c="cir"
     140             :       end select
     141           0 :       DO nc=1,this%ncore(n_dos)
     142           0 :         ind=ind+1
     143           0 :         if (ind==id) THEN
     144           0 :           write(get_weight_name,"(a,i0,a,i0,a)") "At",this%n_dos_to_type(n_dos),"NC",nc,c
     145           0 :           RETURN
     146           0 :         ELSE IF(ind>id) then
     147           0 :           CALL judft_error("Types_mcd: data not found")
     148             :         ENDIF
     149             :       ENDDO
     150             :     ENDDO
     151             :   end function
     152             : 
     153             : 
     154             : 
     155         672 : SUBROUTINE mcd_init(thisMCD,banddos,input,atoms,kpts,eig)
     156             :   USE m_types_setup
     157             :   USE m_types_kpts
     158             : 
     159             :    IMPLICIT NONE
     160             : 
     161             :    CLASS(t_mcd),          INTENT(INOUT) :: thisMCD
     162             :    TYPE(t_banddos),       INTENT(IN)    :: banddos
     163             : 
     164             :    TYPE(t_input),         INTENT(IN)    :: input
     165             :    TYPE(t_atoms),         INTENT(IN)    :: atoms
     166             :    TYPE(t_kpts),          INTENT(IN)    :: kpts
     167             :    real,INTENT(IN)                      :: eig(:,:,:)
     168             : 
     169             :    integer :: ntype !no of types for which MCD is calculated
     170             : 
     171        1374 :    thisMCD%n_dos_to_type=banddos%dos_typelist
     172         672 :    ntype=size(banddos%dos_typelist)
     173         672 :    thisMCD%name_of_dos="MCD"
     174        2016 :    ALLOCATE (thisMCD%ncore(ntype))
     175        3360 :    ALLOCATE (thisMCD%e_mcd(ntype,input%jspins,29))
     176         672 :    IF (banddos%l_mcd) THEN
     177           0 :       thisMCD%emcd_lo = banddos%e_mcd_lo
     178           0 :       thisMCD%emcd_up = banddos%e_mcd_up
     179           0 :       ALLOCATE (thisMCD%m_mcd(29,(3+1)**2,3*ntype,2))
     180           0 :       ALLOCATE (thisMCD%mcd(3*ntype,29,input%neig,kpts%nkpt,input%jspins) )
     181           0 :       IF (.NOT.banddos%dos) WRITE (*,*) 'For mcd-spectra set banddos%dos=T!'
     182             :    ELSE
     183         672 :       ALLOCATE(thisMCD%dos(0,0,0)) !indicated no DOS should be calculated
     184         672 :       ALLOCATE (thisMCD%m_mcd(1,1,1,1))
     185        2016 :       ALLOCATE (thisMCD%mcd(1,1,1,1,input%jspins))
     186             :    ENDIF
     187             : 
     188         672 :    thisMCD%maxE_mcd = -1000000.0
     189         692 :    thisMCD%ncore = 0
     190       51944 :    thisMCD%e_mcd = 0.0
     191        6012 :    thisMCD%mcd = 0.0
     192        3360 :    thisMCD%m_mcd = CMPLX(0.0,0.0)
     193             : 
     194      406524 :    thisMCD%eig=eig
     195             : 
     196             : 
     197             : 
     198             : 
     199         672 : END SUBROUTINE mcd_init
     200           0 : end module m_types_mcd