tlo.f90 10.6 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
MODULE m_tlo
      USE m_juDFT
!***********************************************************************
!     sets up the extra t-matrix elements due to the local orbitals.
!     only non=zero elements are calculated
!
!     p.kurz jul. 1996
!***********************************************************************
      CONTAINS
16
        SUBROUTINE tlo(atoms,sphhar,jspin,jsp,ntyp,enpara,lh0,input,vr,&
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
             flo,f,g,usdus,uuilon,duilon,ulouilopn, tlmplm )
          !
          !*************** ABBREVIATIONS *****************************************
          ! tuulo      : t-matrix element of the lo and the apw radial fuction
          ! tdulo      : t-matrix element of the lo and the energy derivativ of 
          !              the apw radial fuction
          ! tuloulo    : t-matrix element of two los
          !c***********************************************************************
          !
          USE m_intgr, ONLY : intgr3  
          USE m_gaunt, ONLY: gaunt1
          USE m_types
          IMPLICIT NONE
          TYPE(t_input),INTENT(IN)    :: input
          TYPE(t_sphhar),INTENT(IN)   :: sphhar
          TYPE(t_atoms),INTENT(IN)    :: atoms
          TYPE(t_usdus),INTENT(IN)    :: usdus
          TYPE(t_tlmplm),INTENT(INOUT):: tlmplm
          TYPE(t_enpara),INTENT(IN)   :: enpara
          !     ..
          !     .. Scalar Arguments ..
38
          INTEGER, INTENT (IN) :: jspin,jsp,ntyp ,lh0
39 40 41 42 43 44 45 46 47 48
          !     ..
          !     .. Array Arguments ..
          REAL,    INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd)
          REAL,    INTENT (IN) :: f(atoms%jmtd,2,0:atoms%lmaxd),g(atoms%jmtd,2,0:atoms%lmaxd)
          REAL,    INTENT (IN) :: flo(atoms%jmtd,2,atoms%nlod)
          REAL,    INTENT (IN) :: uuilon(atoms%nlod,atoms%ntypd),duilon(atoms%nlod,atoms%ntypd)
          REAL,    INTENT (IN) :: ulouilopn(atoms%nlod,atoms%nlod,atoms%ntypd)
          !     ..
          !     .. Local Scalars ..
          COMPLEX ci,cil
49
          INTEGER i,l,lh,lm ,lmin,lmp,lo,lop,loplo,lp,lpmax,lpmax0,lpmin,lpmin0,lpp ,mem,mp,mpp,m,lmx,mlo,mlolo
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
          !     ..
          !     .. Local Arrays ..
          REAL x(atoms%jmtd),ulovulo(atoms%nlod*(atoms%nlod+1)/2,lh0:sphhar%nlhd)
          REAL uvulo(atoms%nlod,0:atoms%lmaxd,lh0:sphhar%nlhd),dvulo(atoms%nlod,0:atoms%lmaxd,lh0:sphhar%nlhd)
          !     ..

          ci = CMPLX(0.0,1.0)
          DO lo = 1,atoms%nlo(ntyp)
             l = atoms%llo(lo,ntyp)
             DO lp = 0,atoms%lmax(ntyp)
                lmin = ABS(lp-l)
                !               lmin = lp - l
                lmx = lp + l
                DO lh = lh0,sphhar%nlh(atoms%ntypsy(atoms%nat))
                   lpp = sphhar%llh(lh,atoms%ntypsy(atoms%nat))
                   IF ((MOD(l+lp+lpp,2).EQ.1) .OR. (lpp.LT.lmin) .OR.&
                        (lpp.GT.lmx)) THEN
                      uvulo(lo,lp,lh) = 0.0
                      dvulo(lo,lp,lh) = 0.0
                   ELSE
                      DO i = 1,atoms%jri(ntyp)
                         x(i) = (f(i,1,lp)*flo(i,1,lo)+ f(i,2,lp)*flo(i,2,lo))*vr(i,lh)
                      END DO
                      CALL intgr3(x,atoms%rmsh(:,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),uvulo(lo,lp,lh))
                      DO i = 1,atoms%jri(ntyp)
                         x(i) = (g(i,1,lp)*flo(i,1,lo)+ g(i,2,lp)*flo(i,2,lo))*vr(i,lh)
                      END DO
                      CALL intgr3(x,atoms%rmsh(:,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),dvulo(lo,lp,lh))
                   END IF
                END DO
             END DO
          END DO
          loplo = 0
          DO lop = 1,atoms%nlo(ntyp)
             lp = atoms%llo(lop,ntyp)
             DO lo = 1,lop
                l = atoms%llo(lo,ntyp)
                loplo = loplo + 1
                IF (loplo>size(ulovulo,1))  CALL juDFT_error("loplo too large!!!" ,calledby ="tlo")
                DO lh = lh0,sphhar%nlh(atoms%ntypsy(atoms%nat))
                   lpp = sphhar%llh(lh,atoms%ntypsy(atoms%nat))
                   lmin = ABS(lp - l)
                   lmx = lp + l
                   IF ((MOD(l+lp+lpp,2).EQ.1).OR.(lpp.LT.lmin).OR.(lpp.GT.lmx)) THEN
                      ulovulo(loplo,lh) = 0.0
                   ELSE
                      DO i = 1,atoms%jri(ntyp)
                         x(i) = (flo(i,1,lop)*flo(i,1,lo)+flo(i,2,lop)*flo(i,2,lo))*vr(i,lh)
                      END DO
                      CALL intgr3(x,atoms%rmsh(:,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),ulovulo(loplo,lh))
                   END IF
                END DO
             END DO
          END DO
          !---> generate the different t matrices
          !---> but first initialize them ( done in eigen )
          !     
          !---> generate the t-matrices. for optimal performance consider only
          !---> those combinations of l,l',l'',m,m',m'' that satisfy the three
          !---> conditions for non-zero gaunt-coeff. i.e.
          !---> |l - l''| <= l' <= l + l'' (triangular condition)
          !---> m' = m + m'' and l + l' + l'' even
          !---> loop over the local orbitals
113
          mlo=sum(atoms%nlo(:ntyp-1))
114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
          DO lo = 1,atoms%nlo(ntyp)
             l = atoms%llo(lo,ntyp)
             DO m = -l,l
                !--->       loop over the lattice harmonics
                DO lh = lh0,sphhar%nlh(atoms%ntypsy(atoms%nat))
                   lpp = sphhar%llh(lh,atoms%ntypsy(atoms%nat))
                   lpmin0 = ABS(l-lpp)
                   lpmax0 = l + lpp
                   !--->          check that lpmax is smaller than the max l of the
                   !--->          wavefunction expansion at this atom
                   lpmax = MIN(lpmax0,atoms%lmax(ntyp))
                   !--->          make sure that l + l'' + lpmax is even
                   lpmax = lpmax - MOD(l+lpp+lpmax,2)
                   DO mem = 1,sphhar%nmem(lh,atoms%ntypsy(atoms%nat))
                      mpp = sphhar%mlh(mem,lh,atoms%ntypsy(atoms%nat))
                      mp = m + mpp
                      lpmin = MAX(lpmin0,ABS(mp))
                      !--->             make sure that l + l'' + lpmin is even
                      lpmin = lpmin + MOD(ABS(lpmax-lpmin),2)
                      !--->             loop over l'
                      DO lp = lpmin,lpmax,2
                         lmp = lp* (lp+1) + mp
                         cil = ((ci** (l-lp))*sphhar%clnu(mem,lh,atoms%ntypsy(atoms%nat)))* gaunt1(lp,lpp,l,mp,mpp,m,atoms%lmaxd)
137 138 139 140
                         tlmplm%tuulo(lmp,m,lo+mlo,jsp) = &
                              tlmplm%tuulo(lmp,m,lo+mlo,jsp) + cil*uvulo(lo,lp,lh)
                         tlmplm%tdulo(lmp,m,lo+mlo,jsp) = &
                              tlmplm%tdulo(lmp,m,lo+mlo,jsp) + cil*dvulo(lo,lp,lh)
141 142 143 144 145 146 147
                      END DO
                   END DO
                END DO
             END DO
          END DO
          !---> generate the t-matrix including two local orbitals for lo' >= lo
          !---> loop over lo'
148
          mlolo=dot_product(atoms%nlo(:ntyp-1),atoms%nlo(:ntyp-1)+1)/2
149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
          DO lop = 1,atoms%nlo(ntyp)
             lp = atoms%llo(lop,ntyp)
             DO mp = -lp,lp
                !--->       loop over the lattice harmonics
                DO lh = lh0,sphhar%nlh(atoms%ntypsy(atoms%nat))
                   lpp = sphhar%llh(lh,atoms%ntypsy(atoms%nat))
                   DO mem = 1,sphhar%nmem(lh,atoms%ntypsy(atoms%nat))
                      mpp = sphhar%mlh(mem,lh,atoms%ntypsy(atoms%nat))
                      m = mp - mpp
                      !--->             loop over lo
                      DO lo = 1,lop
                         l = atoms%llo(lo,ntyp)
                         loplo = ((lop-1)*lop)/2 + lo
                         IF ((ABS(l-lpp).LE.lp) .AND. (lp.LE. (l+lpp)) .AND.&
                              (MOD(l+lp+lpp,2).EQ.0) .AND. (ABS(m).LE.l)) THEN
                            cil = ((ci** (l-lp))*sphhar%clnu(mem,lh,atoms%ntypsy(atoms%nat)))* gaunt1(lp,lpp,l,mp,mpp,m,atoms%lmaxd)
165
                            tlmplm%tuloulo(mp,m,loplo+mlolo,jsp) = tlmplm%tuloulo(mp,m,loplo+mlolo,jsp) + cil*ulovulo(loplo,lh)
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
                         END IF
                      END DO
                   END DO
                END DO
             END DO
          END DO
          !---> add the diagonal terms from the muffin-tin hamiltonian. these
          !---> terms have to be made hermitian. if second variation is switched
          !---> on, the t-matrices contain only the contributions from the
          !---> non-spherical hamiltonian.
          IF (.NOT.input%secvar) THEN
             DO lo = 1,atoms%nlo(ntyp)
                l = atoms%llo(lo,ntyp)
                DO m = -l,l
                   lm = l* (l+1) + m
181
                   tlmplm%tuulo(lm,m,lo+mlo,jsp) = tlmplm%tuulo(lm,m,lo+mlo,jsp) + 0.5 * usdus%uulon(lo,ntyp,jspin) *&
182
                        ( enpara%el0(l,ntyp,jspin)+enpara%ello0(lo,ntyp,jspin) )
183 184
                   tlmplm%tdulo(lm,m,lo+mlo,jsp) = tlmplm%tdulo(lm,m,lo+mlo,jsp) + 0.5 * usdus%dulon(lo,ntyp,jspin) *&
                        ( enpara%el0(l,ntyp,jspin)+enpara%ello0(lo,ntyp,jspin) ) + 0.5 * usdus%uulon(lo,ntyp,jspin)
185
                   IF (atoms%ulo_der(lo,ntyp).GE.1) THEN
186 187
                      tlmplm%tuulo(lm,m,lo+mlo,jsp) = tlmplm%tuulo(lm,m,lo+mlo,jsp) + 0.5 * uuilon(lo,ntyp)
                      tlmplm%tdulo(lm,m,lo+mlo,jsp) = tlmplm%tdulo(lm,m,lo+mlo,jsp) + 0.5 * duilon(lo,ntyp)
188 189 190
                   ENDIF
                   !+apw_lo
                   IF (atoms%l_dulo(lo,ntyp)) THEN         
191 192
                      tlmplm%tuulo(lm,m,lo+mlo,jsp) = tlmplm%tuulo(lm,m,lo+mlo,jsp) + 0.5
                      tlmplm%tdulo(lm,m,lo+mlo,jsp) = 0.0
193 194 195 196 197 198 199 200 201
                   ENDIF
                   !+apw_lo
                 END DO
             END DO
             DO lop = 1,atoms%nlo(ntyp)
                lp = atoms%llo(lop,ntyp)
                DO lo = atoms%lo1l(lp,ntyp),lop
                   loplo = ((lop-1)*lop)/2 + lo
                   DO m = -lp,lp
202
                      tlmplm%tuloulo(m,m,loplo+mlolo,jsp) = tlmplm%tuloulo(m,m,loplo+mlolo,jsp) + 0.5* (enpara%ello0(lop,ntyp,jspin)+&
203
                           enpara%ello0(lo,ntyp,jspin))* usdus%uloulopn(lop,lo,ntyp,jspin) + 0.5* (ulouilopn(lop,lo,ntyp) +&
204 205 206 207 208 209 210 211
                           ulouilopn(lo,lop,ntyp))
                   END DO
                END DO
             END DO
          END IF

        END SUBROUTINE tlo
      END MODULE m_tlo