eparas.f90 8.26 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
MODULE m_eparas
  !***********************************************************************
  ! Calculates qlo, enerlo and sqlo, which are needed to determine the 
  ! new energy parameters.
  ! Philipp Kurz 99/04
  !***********************************************************************
  ! also the 'normal' energy parameters are now included...
  !
  ! if (l_mcd) then mcd contains mcd spectrum: first index = polarization
  ! second = core level ; third = band index                  gb.2001
  ! corrected to work also for multiple LO's of same l at the same atom
  !                                                           gb.2005
  !*************** ABBREVIATIONS *****************************************
  ! qlo     : charge density of one local orbital at the current k-point
  ! sqlo    : qlo integrated over the Brillouin zone
  ! enerlo  : qlo*energy integrated over the Brillouin zone
  !***********************************************************************
  !
CONTAINS
26
  SUBROUTINE eparas(jsp,atoms,noccbd, mpi,ikpt,ne,we,eig,skip_t,l_evp,eigVecCoeffs,&
27
                    usdus,regCharges,mcd,l_mcd)
28 29
    USE m_types
    IMPLICIT NONE
30 31 32 33 34 35
    TYPE(t_usdus),         INTENT(IN)    :: usdus
    TYPE(t_mpi),           INTENT(IN)    :: mpi
    TYPE(t_atoms),         INTENT(IN)    :: atoms
    TYPE(t_eigVecCoeffs),  INTENT(IN)    :: eigVecCoeffs
    TYPE(t_regionCharges), INTENT(INOUT) :: regCharges
    TYPE(t_mcd),           INTENT(INOUT) :: mcd
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: noccbd,jsp     
    INTEGER, INTENT (IN) :: ne,ikpt  ,skip_t
    LOGICAL, INTENT (IN) :: l_mcd,l_evp
    !     ..
    !     .. Array Arguments ..
    REAL,    INTENT (IN)  :: eig(:)!(dimension%neigd),
    REAL,    INTENT (IN)  :: we(noccbd) 

    !     ..
    !     .. Local Scalars ..
    INTEGER i,l,lo,lop ,natom,nn,ntyp,m
    INTEGER nt1,nt2,lm,n,ll1,ipol,icore,index
    REAL fac
    COMPLEX suma,sumb,sumab,sumba
    !     ..
    !     .. Local Arrays ..
Daniel Wortmann's avatar
Daniel Wortmann committed
54 55
    REAL qlo(noccbd,atoms%nlod,atoms%nlod,atoms%ntype)
    REAL qaclo(noccbd,atoms%nlod,atoms%ntype),qbclo(noccbd,atoms%nlod,atoms%ntype)
56 57 58 59 60 61 62
    !     ..
    !
    !---> initialize ener, sqal, enerlo and sqlo on first call
    !

    IF ((ikpt.LE.mpi%isize).AND..NOT.l_evp) THEN
       IF (l_mcd) THEN
63
          mcd%mcd(:,:,:) = 0.0
64
       ENDIF
65 66 67 68 69
       regCharges%ener(:,:,jsp) = 0.0
       regCharges%sqal(:,:,jsp) = 0.0
       regCharges%qal(:,:,:,jsp) = 0.0
       regCharges%enerlo(:,:,jsp) = 0.0
       regCharges%sqlo(:,:,jsp) = 0.0
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
    END IF
    !
    !--->    l-decomposed density for each occupied state
    !
    !         DO 140 i = (skip_t+1),ne    ! this I need for all states
    DO i = 1,ne              ! skip in next loop
       nt1 = 1
       DO n = 1,atoms%ntype
          fac = 1./atoms%neq(n)
          nt2 = nt1 + atoms%neq(n) - 1
          DO l = 0,3
             suma = CMPLX(0.,0.)
             sumb = CMPLX(0.,0.)
             ll1 = l* (l+1)
             DO m = -l,l
                lm = ll1 + m
                IF ( .NOT.l_mcd ) THEN
                   DO natom = nt1,nt2
88 89
                      suma = suma + eigVecCoeffs%acof(i,lm,natom,jsp)*CONJG(eigVecCoeffs%acof(i,lm,natom,jsp))
                      sumb = sumb + eigVecCoeffs%bcof(i,lm,natom,jsp)*CONJG(eigVecCoeffs%bcof(i,lm,natom,jsp))
90 91 92 93 94
                   ENDDO
                ELSE
                   suma = CMPLX(0.,0.) ; sumab = CMPLX(0.,0.) 
                   sumb = CMPLX(0.,0.) ; sumba = CMPLX(0.,0.)
                   DO natom = nt1,nt2
95 96 97 98
                      suma = suma + eigVecCoeffs%acof(i,lm,natom,jsp)*CONJG(eigVecCoeffs%acof(i,lm,natom,jsp))
                      sumb = sumb + eigVecCoeffs%bcof(i,lm,natom,jsp)*CONJG(eigVecCoeffs%bcof(i,lm,natom,jsp))
                      sumab= sumab + eigVecCoeffs%acof(i,lm,natom,jsp) *CONJG(eigVecCoeffs%bcof(i,lm,natom,jsp))
                      sumba= sumba + eigVecCoeffs%bcof(i,lm,natom,jsp) *CONJG(eigVecCoeffs%acof(i,lm,natom,jsp))
99
                   ENDDO
100
                   DO icore = 1, mcd%ncore(n)
101 102
                      DO ipol = 1, 3
                         index = 3*(n-1) + ipol
103 104 105 106 107
                         mcd%mcd(index,icore,i)=mcd%mcd(index,icore,i) + fac*(&
                              suma * CONJG(mcd%m_mcd(icore,lm+1,index,1))*mcd%m_mcd(icore,lm+1,index,1)  +&
                              sumb * CONJG(mcd%m_mcd(icore,lm+1,index,2))*mcd%m_mcd(icore,lm+1,index,2)  +&
                              sumab* CONJG(mcd%m_mcd(icore,lm+1,index,2))*mcd%m_mcd(icore,lm+1,index,1)  +&
                              sumba* CONJG(mcd%m_mcd(icore,lm+1,index,1))*mcd%m_mcd(icore,lm+1,index,2)  ) 
108 109 110 111
                      ENDDO
                   ENDDO
                ENDIF     ! end MCD
             ENDDO
112
             regCharges%qal(l,n,i,jsp) = (suma+sumb*usdus%ddn(l,n,jsp))/atoms%neq(n)
113 114 115 116 117 118 119 120 121 122 123 124
          ENDDO
          nt1 = nt1 + atoms%neq(n)
       ENDDO
    ENDDO
    !
    !--->    perform Brillouin zone integration and summation over the
    !--->    bands in order to determine the energy parameters for each
    !--->    atom and angular momentum
    !
    DO l = 0,3
       DO n = 1,atoms%ntype
          DO i = (skip_t+1),noccbd
125 126
             regCharges%ener(l,n,jsp) = regCharges%ener(l,n,jsp) + regCharges%qal(l,n,i,jsp)*we(i)*eig(i)
             regCharges%sqal(l,n,jsp) = regCharges%sqal(l,n,jsp) + regCharges%qal(l,n,i,jsp)*we(i)
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
          ENDDO
       ENDDO
    ENDDO

    !---> initialize qlo

    qlo(:,:,:,:) = 0.0
    qaclo(:,:,:) = 0.0
    qbclo(:,:,:) = 0.0

    !---> density for each local orbital and occupied state

    natom = 0
    DO ntyp = 1,atoms%ntype
       DO nn = 1,atoms%neq(ntyp)
          natom = natom + 1
          DO lo = 1,atoms%nlo(ntyp)
             l = atoms%llo(lo,ntyp)
             ll1 = l* (l+1)
             DO m = -l,l
                lm = ll1 + m
                DO i = 1,ne
                   qbclo(i,lo,ntyp) = qbclo(i,lo,ntyp) +REAL(&
150 151
                        eigVecCoeffs%bcof(i,lm,natom,jsp)*CONJG(eigVecCoeffs%ccof(m,i,lo,natom,jsp))+&
                        eigVecCoeffs%ccof(m,i,lo,natom,jsp)*CONJG(eigVecCoeffs%bcof(i,lm,natom,jsp)) )
152
                   qaclo(i,lo,ntyp) = qaclo(i,lo,ntyp) + REAL(&
153 154
                        eigVecCoeffs%acof(i,lm,natom,jsp)*CONJG(eigVecCoeffs%ccof(m,i,lo,natom,jsp))+&
                        eigVecCoeffs%ccof(m,i,lo,natom,jsp)*CONJG(eigVecCoeffs%acof(i,lm,natom,jsp)) )
155 156 157 158 159 160 161
                ENDDO
             ENDDO
             DO lop = 1,atoms%nlo(ntyp)
                IF (atoms%llo(lop,ntyp).EQ.l) THEN
                   DO m = -l,l
                      DO i = 1,ne
                         qlo(i,lop,lo,ntyp) = qlo(i,lop,lo,ntyp) +  REAL(&
162
                              CONJG(eigVecCoeffs%ccof(m,i,lop,natom,jsp))*eigVecCoeffs%ccof(m,i,lo,natom,jsp))
163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178
                      ENDDO
                   ENDDO
                ENDIF
             ENDDO
          ENDDO
       ENDDO
    ENDDO

    !---> perform brillouin zone integration and sum over bands

    DO ntyp = 1,atoms%ntype
       DO lo = 1,atoms%nlo(ntyp)
          l = atoms%llo(lo,ntyp)
          ! llo > 3 used for unoccupied states only
          IF( l .GT. 3 ) CYCLE
          DO i = 1,ne
179
             regCharges%qal(l,ntyp,i,jsp)= regCharges%qal(l,ntyp,i,jsp)  + ( 1.0/atoms%neq(ntyp) )* (&
180 181 182 183 184
                  qaclo(i,lo,ntyp)*usdus%uulon(lo,ntyp,jsp)+qbclo(i,lo,ntyp)*usdus%dulon(lo,ntyp,jsp)     )
          END DO
          DO lop = 1,atoms%nlo(ntyp)
             IF (atoms%llo(lop,ntyp).EQ.l) THEN
                DO i = 1,ne
185 186 187
                   regCharges%enerlo(lo,ntyp,jsp) = regCharges%enerlo(lo,ntyp,jsp) +qlo(i,lop,lo,ntyp)*we(i)*eig(i)
                   regCharges%sqlo(lo,ntyp,jsp) = regCharges%sqlo(lo,ntyp,jsp) + qlo(i,lop,lo,ntyp)*we(i)
                   regCharges%qal(l,ntyp,i,jsp)= regCharges%qal(l,ntyp,i,jsp)  + ( 1.0/atoms%neq(ntyp) ) *&
188 189 190 191 192 193 194 195 196
                        qlo(i,lop,lo,ntyp)*usdus%uloulopn(lop,lo,ntyp,jsp)
                ENDDO
             ENDIF
          ENDDO
       END DO
    END DO

  END SUBROUTINE eparas
END MODULE m_eparas