tlo.f90 10.5 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
                       na,flo,f,g,usdus,uuilon,duilon,ulouilopn, tlmplm )
18 19 20 21 22 23 24 25 26 27 28
          !
          !*************** 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
29
          use m_constants
30 31 32 33 34 35 36 37 38
          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 ..
39
          INTEGER, INTENT (IN) :: jspin,jsp,ntyp ,lh0,na
40 41 42 43 44
          !     ..
          !     .. 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)
Daniel Wortmann's avatar
Daniel Wortmann committed
45 46
          REAL,    INTENT (IN) :: uuilon(atoms%nlod,atoms%ntype),duilon(atoms%nlod,atoms%ntype)
          REAL,    INTENT (IN) :: ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype)
47 48
          !     ..
          !     .. Local Scalars ..
49
          COMPLEX cil
50
          INTEGER i,l,lh,lm ,lmin,lmp,lo,lop,loplo,lp,lpmax,lpmax0,lpmin,lpmin0,lpp ,mem,mp,mpp,m,lmx,mlo,mlolo
51 52 53 54 55 56 57 58 59 60 61 62
          !     ..
          !     .. 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)
          !     ..

          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
63 64
                DO lh = lh0,sphhar%nlh(atoms%ntypsy(na))
                   lpp = sphhar%llh(lh,atoms%ntypsy(na))
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
                   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")
89 90
                DO lh = lh0,sphhar%nlh(atoms%ntypsy(na))
                   lpp = sphhar%llh(lh,atoms%ntypsy(na))
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
                   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
          DO lo = 1,atoms%nlo(ntyp)
             l = atoms%llo(lo,ntyp)
             DO m = -l,l
                !--->       loop over the lattice harmonics
118 119
                DO lh = lh0,sphhar%nlh(atoms%ntypsy(na))
                   lpp = sphhar%llh(lh,atoms%ntypsy(na))
120 121 122 123 124 125 126
                   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)
127 128
                   DO mem = 1,sphhar%nmem(lh,atoms%ntypsy(na))
                      mpp = sphhar%mlh(mem,lh,atoms%ntypsy(na))
129 130 131 132 133 134 135
                      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
136
                         cil = ((ImagUnit** (l-lp))*sphhar%clnu(mem,lh,atoms%ntypsy(na)))* 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
          DO lop = 1,atoms%nlo(ntyp)
             lp = atoms%llo(lop,ntyp)
             DO mp = -lp,lp
                !--->       loop over the lattice harmonics
153 154 155 156
                DO lh = lh0,sphhar%nlh(atoms%ntypsy(na))
                   lpp = sphhar%llh(lh,atoms%ntypsy(na))
                   DO mem = 1,sphhar%nmem(lh,atoms%ntypsy(na))
                      mpp = sphhar%mlh(mem,lh,atoms%ntypsy(na))
157 158 159 160 161 162 163
                      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
164
                            cil = ((ImagUnit** (l-lp))*sphhar%clnu(mem,lh,atoms%ntypsy(na)))* 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