wann_mmk0_updown_sph.f 9.8 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
      MODULE m_wann_mmk0_updown_sph
      CONTAINS
      SUBROUTINE wann_mmk0_updown_sph(
     >               l_noco,alph,beta,
11
     >               llod,noccbd,nlod,natd,ntypd,lmaxd,lmax,lmd,
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
     >               ntype,neq,nlo,llo,
     >               radial1_ff,radial1_gg,
     >               radial1_fg,radial1_gf,
     >               radial1_flo,radial1_glo,
     >               radial1_lof,radial1_log, 
     >               radial1_lolo,
     >               acof,bcof,ccof,
     >               ddn,uulon,dulon,uloulopn,
     =               mmn)
c************************************************************
c     Overlaps of the spin-down parts of the Bloch functions
c     with the spin-up parts.
c                           Frank Freimuth
c************************************************************
      implicit none
      logical, intent (in)  :: l_noco
      integer, intent (in)  :: llod,nlod,natd,ntypd,lmaxd,lmd
29
      integer, intent (in)  :: lmax(:) !(ntypd)
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
      integer, intent (in)  :: ntype,noccbd
      REAL,    INTENT (IN)  :: alph(ntypd),beta(ntypd)
      integer, intent (in)  :: neq(ntypd)
      integer, intent (in)  :: nlo(ntypd),llo(nlod,ntypd)
      real,    intent (in)  :: radial1_ff(:,:,0:,0:,:)
      real,    intent (in)  :: radial1_gg(:,:,0:,0:,:)      
      real,    intent (in)  :: radial1_fg(:,:,0:,0:,:)
      real,    intent (in)  :: radial1_gf(:,:,0:,0:,:)
      real,intent(in)       :: radial1_flo(:,:,0:,:,:)
      real,intent(in)       :: radial1_glo(:,:,0:,:,:)
      real,intent(in)       :: radial1_lof(:,:,:,0:,:)
      real,intent(in)       :: radial1_log(:,:,:,0:,:)
      real,intent(in)       :: radial1_lolo(:,:,:,:,:)
      real,    intent (in)  :: ddn(0:lmaxd,ntypd,2)
      real,    intent (in)  :: uloulopn(nlod,nlod,ntypd,2)
      real,    intent (in)  :: uulon(nlod,ntypd,2),dulon(nlod,ntypd,2)
      complex, intent (in)  :: ccof(-llod:llod,noccbd,nlod,natd,2)
      complex, intent (in)  :: acof(noccbd,0:lmd,natd,2)
      complex, intent (in)  :: bcof(noccbd,0:lmd,natd,2)
      complex, intent (inout) :: mmn(noccbd,noccbd)

      integer           :: i,j,l,lo,lop,m,natom,nn,ntyp
      integer           :: nt1,nt2,lm,n,ll1,i1spin,i2spin
      complex           :: suma,sumb,sumc,sumd
      complex           :: suma12(2,2),sumb12(2,2)
      complex           :: sumc12(2,2),sumd12(2,2)
56 57 58
      complex, allocatable :: qlo(:,:,:,:,:)
      complex, allocatable :: qaclo(:,:,:,:),qbclo(:,:,:,:)
      complex, allocatable :: qcloa(:,:,:,:),qclob(:,:,:,:)
59 60 61 62 63
      COMPLEX           :: ccchi(2,2),ci

      ci = cmplx(0.0,1.0)
      allocate (qlo(noccbd,noccbd,nlod,nlod,ntypd), 
     +          qaclo(noccbd,noccbd,nlod,ntypd),
64 65 66 67
     +          qbclo(noccbd,noccbd,nlod,ntypd),
     +          qcloa(noccbd,noccbd,nlod,ntypd),
     +          qclob(noccbd,noccbd,nlod,ntypd))
 
68 69 70 71 72 73 74 75 76 77 78 79
c---> performs summations of the overlaps of the wavefunctions
      do i = 1,noccbd            
       do j = 1,noccbd
         nt1 = 1
         do n = 1,ntype
            if(l_noco)then
               ccchi(1,1) = conjg( exp( ci*alph(n)/2)*cos(beta(n)/2))
               ccchi(1,2) = conjg(-exp( ci*alph(n)/2)*sin(beta(n)/2))
               ccchi(2,1) = conjg( exp(-ci*alph(n)/2)*sin(beta(n)/2))
               ccchi(2,2) = conjg( exp(-ci*alph(n)/2)*cos(beta(n)/2))
            endif
            nt2 = nt1 + neq(n) - 1
80
            do l = 0,lmax(n)
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 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158
             if(.not.l_noco)then  
               suma = cmplx(0.,0.)
               sumb = cmplx(0.,0.)
               sumc = cmplx(0.,0.)
               sumd = cmplx(0.,0.)
               ll1 = l* (l+1)
               do m = -l,l
                  lm = ll1 + m
                  do natom = nt1,nt2
                    suma = suma + acof(i,lm,natom,1)*
     +                      conjg(acof(j,lm,natom,2))
                    sumb = sumb + bcof(i,lm,natom,1)*
     +                      conjg(bcof(j,lm,natom,2))
                    sumc = sumc + acof(i,lm,natom,1)*
     +                      conjg(bcof(j,lm,natom,2))
                    sumd = sumd + bcof(i,lm,natom,1)*
     +                      conjg(acof(j,lm,natom,2))
                  enddo !natom
               enddo !m      
               mmn(i,j) = mmn(i,j) + ( suma*radial1_ff(1,2,l,l,n)+
     +                                 sumb*radial1_gg(1,2,l,l,n)+
     +                                 sumc*radial1_fg(1,2,l,l,n)+
     +                                 sumd*radial1_gf(1,2,l,l,n)  )      
             else
               suma12 = cmplx(0.,0.)
               sumb12 = cmplx(0.,0.)
               sumc12 = cmplx(0.,0.)
               sumd12 = cmplx(0.,0.)
               ll1 = l* (l+1)
               do i1spin=1,2
                do i2spin=1,2
                 do m = -l,l
                  lm = ll1 + m
                  do natom = nt1,nt2
                    suma12(i1spin,i2spin) = suma12(i1spin,i2spin) 
     +                      + acof(i,lm,natom,i1spin)*
     +                      conjg(acof(j,lm,natom,i2spin))
                    sumb12(i1spin,i2spin) = sumb12(i1spin,i2spin) 
     +                      + bcof(i,lm,natom,i1spin)*
     +                      conjg(bcof(j,lm,natom,i2spin))
                    sumc12(i1spin,i2spin) = sumc12(i1spin,i2spin) 
     +                      + acof(i,lm,natom,i1spin)*
     +                      conjg(bcof(j,lm,natom,i2spin))
                    sumd12(i1spin,i2spin) = sumd12(i1spin,i2spin) 
     +                      + bcof(i,lm,natom,i1spin)*
     +                      conjg(acof(j,lm,natom,i2spin))
                  enddo !natom
                 enddo !m
                 mmn(i,j) =    mmn(i,j)
     &                           +
     &             suma12(i1spin,i2spin)*radial1_ff(i1spin,i2spin,l,l,n)
     &                    *ccchi(1,i2spin)*conjg(ccchi(2,i1spin))

     &                           +
     &             sumb12(i1spin,i2spin)*radial1_gg(i1spin,i2spin,l,l,n)
     &                    *ccchi(1,i2spin)*conjg(ccchi(2,i1spin))

     &                           +
     &             sumc12(i1spin,i2spin)*radial1_fg(i1spin,i2spin,l,l,n)
     &                    *ccchi(1,i2spin)*conjg(ccchi(2,i1spin))

     &                           +
     &             sumd12(i1spin,i2spin)*radial1_gf(i1spin,i2spin,l,l,n)
     &                    *ccchi(1,i2spin)*conjg(ccchi(2,i1spin))
                enddo !i2spin
               enddo !i1spin 
             endif   

            enddo !l
            nt1 = nt1 + neq(n)
         enddo !n   
       enddo   ! cycle by j-band
      enddo  !  cycle by i-band

c---> initialize qlo arrays
      qlo(:,:,:,:,:) = 0.0
      qaclo(:,:,:,:) = 0.0
      qbclo(:,:,:,:) = 0.0
159 160 161
      qcloa(:,:,:,:) = 0.0
      qclob(:,:,:,:) = 0.0

162 163 164 165 166 167 168 169 170 171
c---> prepare the coefficients
      natom = 0
      do ntyp = 1,ntype
         do nn = 1,neq(ntyp)
            natom = natom + 1
            do lo = 1,nlo(ntyp)
               l = llo(lo,ntyp)
               ll1 = l* (l+1)
               do m = -l,l
                  lm = ll1 + m
172 173 174 175 176 177 178 179 180 181 182 183 184 185 186
                  do j = 1,noccbd
                   do i = 1,noccbd
                     qbclo(i,j,lo,ntyp) = qbclo(i,j,lo,ntyp) +
     +                  bcof(i,lm,natom,1)*conjg(ccof(m,j,lo,natom,2)) 
c     +                +ccof(m,i,lo,natom,1)*conjg(bcof(j,lm,natom,2)) 
                     qaclo(i,j,lo,ntyp) = qaclo(i,j,lo,ntyp) +
     +                  acof(i,lm,natom,1)*conjg(ccof(m,j,lo,natom,2)) 
c     +                 +ccof(m,i,lo,natom,1)*conjg(acof(j,lm,natom,2))

                     qclob(i,j,lo,ntyp) = qclob(i,j,lo,ntyp)
c     +                bcof(i,lm,natom,1)*conjg(ccof(m,j,lo,natom,2)) 
     +                +ccof(m,i,lo,natom,1)*conjg(bcof(j,lm,natom,2)) 
                     qcloa(i,j,lo,ntyp) = qcloa(i,j,lo,ntyp) 
c     +                acof(i,lm,natom,1)*conjg(ccof(m,j,lo,natom,2))+ 
     +                +ccof(m,i,lo,natom,1)*conjg(acof(j,lm,natom,2))
187 188 189 190 191 192 193 194
                   enddo
                  enddo
               enddo
               do lop = 1,nlo(ntyp)
                 if (llo(lop,ntyp).eq.l) then
                   do m = -l,l
                     do i = 1,noccbd
                      do j = 1,noccbd
195 196 197
                       qlo(i,j,lo,lop,ntyp) = qlo(i,j,lo,lop,ntyp) + 
     +                        conjg(ccof(m,j,lop,natom,2))
     *                                  *ccof(m,i,lo,natom,1)
198 199 200 201 202 203 204 205 206 207 208 209 210
                      enddo
                     enddo
                   enddo
                 endif
               enddo
            enddo
         enddo
      enddo
c---> perform summation of the coefficients with the integrals
c---> of the radial basis functions
      do ntyp = 1,ntype
         do lo = 1,nlo(ntyp)
            l = llo(lo,ntyp)
211 212
            do j = 1,noccbd
             do i = 1,noccbd
213
               mmn(i,j)= mmn(i,j)  + 
214 215 216 217 218 219 220
     +    qaclo(i,j,lo,ntyp)*radial1_flo(1,2,l,lo,ntyp)     
     +                      +
     +    qbclo(i,j,lo,ntyp)*radial1_glo(1,2,l,lo,ntyp)
     +                      +
     +    qcloa(i,j,lo,ntyp)*radial1_lof(1,2,lo,l,ntyp)     
     +                      +
     +    qclob(i,j,lo,ntyp)*radial1_log(1,2,lo,l,ntyp)  
221 222 223 224
             enddo
            enddo 
            do lop = 1,nlo(ntyp)
               if (llo(lop,ntyp).eq.l) then
225 226
               do j = 1,noccbd
                do i = 1,noccbd
227
                 mmn(i,j) = mmn(i,j)  + 
228
     +   qlo(i,j,lop,lo,ntyp)*radial1_lolo(1,2,lop,lo,ntyp)
229 230 231 232 233 234
                enddo
               enddo
               endif
            enddo
         enddo 
      enddo 
235
      deallocate ( qlo,qaclo,qbclo,qcloa,qclob )
236 237 238

      END SUBROUTINE wann_mmk0_updown_sph
      END MODULE m_wann_mmk0_updown_sph
239