kp_perturbation.F90 41.3 KB
Newer Older
1 2 3 4 5 6 7
      MODULE  m_kp_perturbation

      CONTAINS

      SUBROUTINE ibs_correction(&
                        nk,atoms,&
                        dimension,input,jsp,&
Daniel Wortmann's avatar
Daniel Wortmann committed
8
                        hybdat,hybrid,&
9
                        lapw,kpts,nkpti,&
Daniel Wortmann's avatar
Daniel Wortmann committed
10 11
                        cell,mnobd,&
                        sym,&
12 13 14 15 16 17 18 19 20 21
                        proj_ibsc,olap_ibsc)

      USE m_sphbes
      USE m_dsphbs
      USE m_constants
      USE m_ylm
      USE m_gaunt
      USE m_apws
      USE m_util
      USE m_types
Daniel Wortmann's avatar
Daniel Wortmann committed
22
      USE m_io_hybrid
23
      IMPLICIT NONE
Daniel Wortmann's avatar
Daniel Wortmann committed
24 25

      TYPE(t_hybdat),INTENT(IN)   :: hybdat
26
      TYPE(t_dimension),INTENT(IN)   :: dimension
Daniel Wortmann's avatar
Daniel Wortmann committed
27
      TYPE(t_hybrid),INTENT(INOUT)   :: hybrid
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
      TYPE(t_input),INTENT(IN)   :: input
      TYPE(t_sym),INTENT(IN)   :: sym
      TYPE(t_cell),INTENT(IN)   :: cell
      TYPE(t_kpts),INTENT(IN)   :: kpts
      TYPE(t_atoms),INTENT(IN)   :: atoms
      TYPE(t_lapw),INTENT(IN)   :: lapw
      
      ! - scalars -
      INTEGER, INTENT(IN)   :: jsp 
      INTEGER, INTENT(IN)   ::  mnobd
      INTEGER, INTENT(IN)   ::  nk  ,nkpti
      
            
      
!       REAL   , INTENT(INOUT)::  ibs_corr(3,3,nbands)
      
      ! - arrays -
      

      
      COMPLEX, INTENT(INOUT)::  olap_ibsc(3,3,mnobd,mnobd)
Daniel Wortmann's avatar
Daniel Wortmann committed
49
      COMPLEX, INTENT(INOUT)::  proj_ibsc(:,:,:)!(3,mnobd,hybdat%nbands(nk))
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
      ! - local scalars -
      INTEGER               ::  i,itype,ieq,iatom,iatom1 ,iband,iband1
      INTEGER               ::  iband2,ilo,ibas,ic,ikpt,ikvec,invsfct
      INTEGER               ::  irecl_cmt,irecl_z
      INTEGER               ::  j,m
      INTEGER               ::  l1,m1,p1,l2,m2,p2,l ,p,lm,&
                                lmp,lmp1,lmp2,lm1,lm2
      INTEGER               ::  ok,ig
      INTEGER               ::  idum
      REAL                  ::  const
      REAL                  ::  ka,kb
      REAL                  ::  kvecn
      REAL                  ::  olap_udot,olap_uulo,olap_udotulo 
      REAL                  ::  rdum
      REAL                  ::  ws
      COMPLEX               ::  phase
      COMPLEX               ::  cj,cdj
      COMPLEX               ::  denom,enum
      COMPLEX               ::  cdum,cdum1,cdum2
      COMPLEX, PARAMETER    ::  img = (0d0,1d0)
      ! - local arrays -
      INTEGER               ::  lmp_start(atoms%ntype)
      REAL                  ::  alo(atoms%nlod,atoms%ntype),blo(atoms%nlod,atoms%ntype),&
                                clo(atoms%nlod,atoms%ntype)
      REAL                  ::  u1_lo(atoms%jmtd,atoms%nlod,atoms%ntype),&
                                u2_lo(atoms%jmtd,atoms%nlod,atoms%ntype)
      REAL                  ::  kvec(3),qvec(3)
      REAL                  ::  sbes(0:atoms%lmaxd+1),dsbes(0:atoms%lmaxd+1)
      REAL                  ::  bas1_tmp(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd+1,atoms%ntype),&
                                bas2_tmp(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd+1,atoms%ntype)
      REAL                  ::  bas1_MT_tmp(hybrid%maxindx,0:atoms%lmaxd+1,atoms%ntype),&
                                drbas1_MT_tmp(hybrid%maxindx,0:atoms%lmaxd+1,atoms%ntype)
      REAL                  ::  ru1(atoms%jmtd,3,mnobd),ru2(atoms%jmtd,3,mnobd)
      REAL                  ::  iu1(atoms%jmtd,3,mnobd),iu2(atoms%jmtd,3,mnobd)
      REAL                  ::  rintegrand(atoms%jmtd),iintegrand(atoms%jmtd),&
                                integrand(atoms%jmtd)
  
      COMPLEX               ::  f(atoms%jmtd,mnobd)
Daniel Wortmann's avatar
Daniel Wortmann committed
88
      COMPLEX               ::  carr(3),carr2(3,hybdat%nbands(nk))
89 90 91 92
      COMPLEX               ::  ylm( (atoms%lmaxd+2)**2 )
      COMPLEX,ALLOCATABLE   ::  u1(:,:,:,:,:),u2(:,:,:,:,:)
      COMPLEX,ALLOCATABLE   ::  cmt_lo(:,:,:,:)
      COMPLEX,ALLOCATABLE   ::  cmt_apw(:,:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
93 94 95
      TYPE(t_mat)           ::  z
      REAL                  ::  work_r(dimension%neigd)
      COMPLEX               ::  work_c(dimension%neigd)
96 97


Daniel Wortmann's avatar
Daniel Wortmann committed
98
      !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
99 100 101

      

Daniel Wortmann's avatar
Daniel Wortmann committed
102 103
      bas1_tmp(:,:,0:atoms%lmaxd,:)    = hybdat%bas1(:,:,0:atoms%lmaxd,:)
      bas2_tmp(:,:,0:atoms%lmaxd,:)    = hybdat%bas2(:,:,0:atoms%lmaxd,:)
104

Daniel Wortmann's avatar
Daniel Wortmann committed
105 106
      bas1_MT_tmp(:,0:atoms%lmaxd,:) =   hybdat%bas1_MT(:,0:atoms%lmaxd,:)
      drbas1_MT_tmp(:,0:atoms%lmaxd,:) = hybdat%drbas1_MT(:,0:atoms%lmaxd,:)
107
       
Daniel Wortmann's avatar
Daniel Wortmann committed
108 109
      bas1_tmp(:,:,atoms%lmaxd+1,:)    = hybdat%bas1(:,:,atoms%lmaxd,:)
      bas2_tmp(:,:,atoms%lmaxd+1,:)    = hybdat%bas2(:,:,atoms%lmaxd,:)
110
        
Daniel Wortmann's avatar
Daniel Wortmann committed
111 112
      bas1_MT_tmp(:,atoms%lmaxd+1,:) =   hybdat%bas1_MT(:,atoms%lmaxd,:)
      drbas1_MT_tmp(:,atoms%lmaxd+1,:) = hybdat%drbas1_MT(:,atoms%lmaxd,:)
113 114 115 116


      
      ! read in z coefficient from direct access file z at k-point nk
Daniel Wortmann's avatar
Daniel Wortmann committed
117 118

      call read_z(z,nk)
119 120 121 122 123 124

      ! construct local orbital consisting of radial function times spherical harmonic
      ! where the radial function vanishes on the MT sphere boundary
      ! with this the local orbitals have a trivial k-dependence

      ! compute radial lo matching coefficients
Daniel Wortmann's avatar
Daniel Wortmann committed
125
      hybrid%nindx = 2
126 127 128
      DO itype = 1,atoms%ntype
        DO ilo = 1,atoms%nlo(itype)
          l              = atoms%llo(ilo,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
129 130
          hybrid%nindx(l,itype) = hybrid%nindx(l,itype) + 1          
          p              = hybrid%nindx(l,itype)
131

Daniel Wortmann's avatar
Daniel Wortmann committed
132
          ws = -wronskian(hybdat%bas1_MT(1,l,itype),hybdat%drbas1_MT(1,l,itype), hybdat%bas1_MT(2,l,itype),hybdat%drbas1_MT(2,l,itype) )
133
 
Daniel Wortmann's avatar
Daniel Wortmann committed
134
          ka = 1.0/ws*wronskian(hybdat%bas1_MT(p,l,itype),hybdat%drbas1_MT(p,l,itype), hybdat%bas1_MT(2,l,itype),hybdat%drbas1_MT(2,l,itype))
135

Daniel Wortmann's avatar
Daniel Wortmann committed
136
          kb = 1.0/ws*wronskian(hybdat%bas1_MT(1,l,itype),hybdat%drbas1_MT(1,l,itype), hybdat%bas1_MT(p,l,itype),hybdat%drbas1_MT(p,l,itype))
137

Daniel Wortmann's avatar
Daniel Wortmann committed
138 139
          integrand    = hybdat%bas1(:,2,l,itype)*hybdat%bas1(:,2,l,itype) + hybdat%bas2(:,2,l,itype)*hybdat%bas2(:,2,l,itype) 
          olap_udot    = intgrf(integrand,atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype, hybdat%gridf)
140
          
Daniel Wortmann's avatar
Daniel Wortmann committed
141 142
          integrand    = hybdat%bas1(:,1,l,itype)*hybdat%bas1(:,p,l,itype) + hybdat%bas2(:,1,l,itype)*hybdat%bas2(:,p,l,itype) 
          olap_uulo    = intgrf(integrand,atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype, hybdat%gridf)
143
          
Daniel Wortmann's avatar
Daniel Wortmann committed
144 145
          integrand    = hybdat%bas1(:,2,l,itype)*hybdat%bas1(:,p,l,itype) + hybdat%bas2(:,2,l,itype)*hybdat%bas2(:,p,l,itype)
          olap_udotulo = intgrf(integrand,atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype, hybdat%gridf)
146 147 148 149 150 151

          rdum           =  ka**2 + (kb**2)*olap_udot + 1.0 +  2.0*ka*olap_uulo + 2.0*kb*olap_udotulo 
          clo(ilo,itype) = 1.0/sqrt( rdum  )
          alo(ilo,itype) = ka*clo(ilo,itype)
          blo(ilo,itype) = kb*clo(ilo,itype)
          
Daniel Wortmann's avatar
Daniel Wortmann committed
152
          u1_lo(:,ilo,itype) =  alo(ilo,itype)*hybdat%bas1(:,1,l,itype) +  blo(ilo,itype)*hybdat%bas1(:,2,l,itype) +  clo(ilo,itype)*hybdat%bas1(:,p,l,itype)
153

Daniel Wortmann's avatar
Daniel Wortmann committed
154
          u2_lo(:,ilo,itype) =  alo(ilo,itype)*hybdat%bas2(:,1,l,itype) +  blo(ilo,itype)*hybdat%bas2(:,2,l,itype) +  clo(ilo,itype)*hybdat%bas2(:,p,l,itype)
155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
         END DO
      END DO



      ! calculate lo wavefunction coefficients
      ALLOCATE( cmt_lo(dimension%neigd,-atoms%llod:atoms%llod,atoms%nlod,atoms%nat) )
      cmt_lo = 0
      iatom  = 0
      ic     = 0
      ibas   = lapw%nv(jsp)
      DO itype = 1,atoms%ntype
        ! the program is in hartree units, therefore 1/wronskian is
        ! (rmt**2)/2.
        const = fpi_const*(atoms%rmt(itype)**2)/2/sqrt(cell%omtil)
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
          IF( (atoms%invsat(iatom).eq.0) .or. (atoms%invsat(iatom) .eq. 1) ) THEN
            IF( atoms%invsat(iatom) .eq. 0 ) invsfct = 1
            IF( atoms%invsat(iatom) .eq. 1 ) THEN
              invsfct = 2
              iatom1 = sym%invsatnr(iatom)
            END IF
 
            DO ilo = 1,atoms%nlo(itype)
              l    = atoms%llo(ilo,itype)
              cdum = img**l*const
              DO ikvec = 1,invsfct*(2*l+1)
                ic    = ic   + 1
                ibas  = ibas + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
185
                kvec  = kpts%bk(:,nk) +(/ lapw%k1(hybdat%kveclo_eig(ic,nk),jsp), lapw%k2(hybdat%kveclo_eig(ic,nk),jsp), lapw%k3(hybdat%kveclo_eig(ic,nk),jsp)/)
186 187 188 189 190 191 192 193 194 195
                
               phase=exp(img*tpi_const*dot_product(atoms%taual(:,iatom),kvec))
                cdum1 = cdum*phase
                
                CALL ylm4(l,matmul(kvec,cell%bmat),ylm)
      
                lm = l**2
                DO M = -l,l
                  lm   = lm + 1
                  cdum2 = cdum1*conjg(ylm(lm))
Daniel Wortmann's avatar
Daniel Wortmann committed
196 197 198 199 200
                if (z%l_real) THEN
                   work_r  = z%data_r(ibas,:)
                   DO iband = 1,hybdat%nbands(nk)
                      cmt_lo(iband,M,ilo,iatom ) = cmt_lo(iband,M,ilo,iatom)  + cdum2*work_r(iband)
                      IF( invsfct .eq. 2 ) THEN
201 202
                      ! the factor (-1)**l is necessary as we do not calculate
                      ! the cmt_lo in the local coordinate system of the atom
Daniel Wortmann's avatar
Daniel Wortmann committed
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218
                         cmt_lo(iband,-M,ilo,iatom1) = cmt_lo(iband,-M,ilo,iatom1) + (-1)**(l+M)*conjg(cdum2)*work_r(iband)
                      END IF
                   END DO
                else
                   work_c  = z%data_c(ibas,:)
                   DO iband = 1,hybdat%nbands(nk)
                      cmt_lo(iband,M,ilo,iatom ) = cmt_lo(iband,M,ilo,iatom)  + cdum2*work_c(iband)
                      IF( invsfct .eq. 2 ) THEN
                      ! the factor (-1)**l is necessary as we do not calculate
                      ! the cmt_lo in the local coordinate system of the atom
                         cmt_lo(iband,-M,ilo,iatom1) = cmt_lo(iband,-M,ilo,iatom1) + (-1)**(l+M)*conjg(cdum2)*work_c(iband)
                      END IF
                   END DO
                end if
                   
             END DO
219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
                
              END DO  !ikvec
            END DO  ! ilo
        
          END IF
        
        END DO  !ieq
      END DO  !itype
      

      
      !
      ! calculate apw wavefunction coefficients up to lmax + 1
      ! note that the lo contribution is separated in cmt_lo
      ! 
    
      DO itype = 1,atoms%ntype
        lmp_start(itype) = sum( (/ (2*(2*l+1), l=0,atoms%lmax(itype)+1) /) )
      END DO
      idum = maxval(lmp_start)


      ALLOCATE( cmt_apw(dimension%neigd,idum,atoms%nat) )
      cmt_apw = 0
      DO i = 1,lapw%nv(jsp)
Daniel Wortmann's avatar
Daniel Wortmann committed
244
         kvec  = kpts%bk(:,nk) + (/ lapw%k1(i,jsp),lapw%k2(i,jsp),lapw%k3(i,jsp) /)
245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282
        kvecn = sqrt(dot_product(matmul(kvec,cell%bmat),matmul(kvec,cell%bmat)))
        
        iatom = 0
        DO itype = 1,atoms%ntype
          !calculate spherical sperical harmonics
          CALL ylm4(atoms%lmax(itype)+1,matmul(kvec,cell%bmat),ylm)

          !calculate spherical bessel function at |kvec|*R_MT(itype)
          CALL sphbes(atoms%lmax(itype)+1,kvecn*atoms%rmt(itype),sbes)
          
          !calculate radial derivative of spherical bessel function at |kvec|*R_MT(itype)
          CALL dsphbs(atoms%lmax(itype)+1,kvecn*atoms%rmt(itype),sbes,dsbes)
          dsbes = kvecn*dsbes
          
          DO ieq = 1,atoms%neq(itype)
            iatom = iatom + 1
                 
            phase = exp(img*tpi_const*dot_product(kvec,atoms%taual(:,iatom)))
            
            lm  = 0
            lmp = 0
            DO l = 0,atoms%lmax(itype)+1
              denom = wronskian(  bas1_MT_tmp(2,l,itype), drbas1_MT_tmp(2,l,itype),&
                                  bas1_MT_tmp(1,l,itype), drbas1_MT_tmp(1,l,itype) )
              cdum1 = fpi_const*img**l*  sbes(l) *phase /sqrt(cell%omtil)
              cdum2 = fpi_const*img**l* dsbes(l) *phase /sqrt(cell%omtil)
              DO M = -l,l
                lm  = lm + 1
                cj  = cdum1 *conjg( ylm(lm) )           
                cdj = cdum2 *conjg( ylm(lm) )  
                DO p = 1,2
                  lmp = lmp + 1
                  p1  = p + (-1)**(p-1)
                  
                  enum  = CMPLX(wronskian(  bas1_MT_tmp(p1,l,itype), drbas1_MT_tmp(p1,l,itype), REAL(cj),REAL(cdj)),&
                                wronskian(  bas1_MT_tmp(p1,l,itype), drbas1_MT_tmp(p1,l,itype), AIMAG(cj),AIMAG(cdj)))
          
                  cdum  = (-1)**(p+1)*enum/denom
Daniel Wortmann's avatar
Daniel Wortmann committed
283 284
                  if (z%l_real) THEN
                  work_r  = z%data_r(i,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
285
                  DO iband = 1,hybdat%nbands(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
286
                    cmt_apw(iband,lmp,iatom) = cmt_apw(iband,lmp,iatom) + cdum*work_r(iband)
287
                  END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
288 289 290 291 292 293
               else
                     work_c  = z%data_c(i,:)
                  DO iband = 1,hybdat%nbands(nk)
                    cmt_apw(iband,lmp,iatom) = cmt_apw(iband,lmp,iatom) + cdum*work_c(iband)
                  END DO
               end if
294 295 296 297 298 299 300 301 302 303 304 305 306 307
                END DO  !p
              END DO  !M
            END DO  !l

          END DO  !iatom
        END DO  ! itype
      END DO  ! i




      ! construct radial functions (complex) for the first order 
      ! incomplete basis set correction 
      
Daniel Wortmann's avatar
Daniel Wortmann committed
308
      ALLOCATE( u1(atoms%jmtd,3,mnobd,(atoms%lmaxd+1)**2,atoms%nat),stat=ok )!hybdat%nbands
309
      IF( ok .ne. 0 ) STOP 'kp_perturbation: failure allocation u1' 
Daniel Wortmann's avatar
Daniel Wortmann committed
310
      ALLOCATE( u2(atoms%jmtd,3,mnobd,(atoms%lmaxd+1)**2,atoms%nat),stat=ok )!hybdat%nbands
311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330
      IF( ok .ne. 0 ) STOP 'kp_perturbation: failure allocation u2' 
      u1 = 0; u2 = 0
      
      iatom = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1

          lm1   = 0
          DO l1 = 0,atoms%lmax(itype)! + 1
            DO m1 = -l1,l1
              lm1 = lm1 + 1

              DO p1 = 1,2

                carr2 = 0
                l2    = l1 + 1
                lmp2  = 2*l2**2 + p1 
                DO m2 = -l2,l2                  
                  carr = gauntvec(l1,m1,l2,m2,atoms)
Daniel Wortmann's avatar
Daniel Wortmann committed
331
                  DO iband = 1,mnobd! hybdat%nbands
332 333 334 335 336
                    carr2(1:3,iband) = carr2(1:3,iband) + carr*cmt_apw(iband,lmp2,iatom)
                  END DO
                  lmp2 = lmp2 + 2
                END DO
                  
Daniel Wortmann's avatar
Daniel Wortmann committed
337
                DO iband = 1,mnobd! hybdat%nbands
338 339 340 341 342 343 344 345 346 347 348 349 350 351 352
                  DO i = 1,3
                    DO ig = 1,atoms%jri(itype)
                      ! the r factor is already included in bas1
                      u1(ig,i,iband,lm1,iatom)= u1(ig,i,iband,lm1,iatom) - img*bas1_tmp(ig,p1,l2,itype)*carr2(i,iband)
                      u2(ig,i,iband,lm1,iatom)= u2(ig,i,iband,lm1,iatom) - img*bas2_tmp(ig,p1,l2,itype)*carr2(i,iband)
                    END DO
                  END DO
                END DO
                
                l2 = l1-1
                IF( l2 .ge. 0 ) THEN                
                  carr2 = 0
                  lmp2 = 2*l2**2 + p1
                  DO m2 = -l2,l2
                    carr = gauntvec(l1,m1,l2,m2,atoms)
Daniel Wortmann's avatar
Daniel Wortmann committed
353
                    DO iband = 1,mnobd! hybdat%nbands
354 355 356 357 358
                      carr2(1:3,iband) = carr2(1:3,iband) + carr*cmt_apw(iband,lmp2,iatom)
                    END DO
                    lmp2 = lmp2 + 2
                  END DO
                  
Daniel Wortmann's avatar
Daniel Wortmann committed
359
                  DO iband = 1,mnobd! hybdat%nbands
360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377
                    DO i = 1,3
                      DO ig = 1,atoms%jri(itype)
                        ! the r factor is already included in bas1
                        u1(ig,i,iband,lm1,iatom) = u1(ig,i,iband,lm1,iatom) - img*bas1_tmp(ig,p1,l2,itype)*carr2(i,iband)
                        u2(ig,i,iband,lm1,iatom) = u2(ig,i,iband,lm1,iatom) - img*bas2_tmp(ig,p1,l2,itype)*carr2(i,iband)
                      END DO
                    END DO
                  END DO
                
                END IF
                
                carr2 = 0
                l2    = l1 + 1
                lmp2  = 2*l2**2
                DO m2 = -l2,l2
                  carr = gauntvec(l1,m1,l2,m2,atoms)
                  DO p2 = 1,2
                    lmp2 = lmp2 + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
378 379
                    rdum = w(p1,l1,p2,l2,itype,bas1_MT_tmp, drbas1_MT_tmp,atoms%rmt)
                    DO iband = 1,mnobd! hybdat%nbands
380 381 382 383 384
                      carr2(1:3,iband) = carr2(1:3,iband) + img*carr*rdum*cmt_apw(iband,lmp2,iatom)
                    END DO
                  END DO
                END DO
                  
Daniel Wortmann's avatar
Daniel Wortmann committed
385
                DO iband = 1,mnobd! hybdat%nbands
386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401
                  DO i = 1,3
                    DO ig = 1,atoms%jri(itype)
                      u1(ig,i,iband,lm1,iatom)= u1(ig,i,iband,lm1,iatom) + bas1_tmp(ig,p1,l1,itype) * carr2(i,iband)/atoms%rmsh(ig,itype)
                      u2(ig,i,iband,lm1,iatom)= u2(ig,i,iband,lm1,iatom) + bas2_tmp(ig,p1,l1,itype) * carr2(i,iband)/atoms%rmsh(ig,itype)
                    END DO
                  END DO
                END DO
                
                l2 = l1 - 1
                IF( l2 .ge. 0 ) THEN 
                  carr2 = 0
                  lmp2  = 2*l2**2
                  DO m2 = -l2,l2
                    carr = gauntvec(l1,m1,l2,m2,atoms)
                    DO p2 = 1,2
                      lmp2 = lmp2 + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
402 403
                      rdum = w(p1,l1,p2,l2,itype,bas1_MT_tmp, drbas1_MT_tmp,atoms%rmt)
                      DO iband = 1,mnobd! hybdat%nbands
404 405 406 407 408
                        carr2(1:3,iband) = carr2(1:3,iband) + img*carr*rdum*cmt_apw(iband,lmp2,iatom)
                      END DO
                    END DO
                  END DO
                  
Daniel Wortmann's avatar
Daniel Wortmann committed
409
                  DO iband = 1,mnobd! hybdat%nbands
410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432
                    DO i = 1,3
                      DO ig = 1,atoms%jri(itype)
                        u1(ig,i,iband,lm1,iatom) = u1(ig,i,iband,lm1,iatom) &
                                         + bas1_tmp(ig,p1,l1,itype) * carr2(i,iband)/atoms%rmsh(ig,itype)
                        u2(ig,i,iband,lm1,iatom) = u2(ig,i,iband,lm1,iatom) &
                                         + bas2_tmp(ig,p1,l1,itype) * carr2(i,iband)/atoms%rmsh(ig,itype)
                      END DO
                    END DO
                  END DO
                END IF
                
              END DO  ! p1
            END DO  !m1
          END DO  !l1
        END DO  !ieq
      END DO  !iatom
      
      ! construct lo contribtution
      IF( any(atoms%llo .eq. atoms%lmaxd ) ) STOP 'ibs_correction: atoms%llo=atoms%lmaxd is not implemented'

      iatom = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
433
          hybrid%nindx = 2
434 435 436
          iatom = iatom + 1
          DO ilo = 1,atoms%nlo(itype)
            l1              = atoms%llo(ilo,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
437 438
            hybrid%nindx(l1,itype) = hybrid%nindx(l1,itype) + 1
            p1              = hybrid%nindx(l1,itype)
439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518
                                                  
            l2  = l1 + 1
            lm2 = l2**2
            DO m2 = -l2,l2
              lm2   = lm2 + 1
              carr2 = 0
  
              DO m1 = -l1,l1
                carr = gauntvec(l2,m2,l1,m1,atoms)
                DO iband = 1,mnobd
                  carr2(1:3,iband) = carr2(1:3,iband) + cmt_lo(iband,m1,ilo,iatom)*carr
                END DO
              END DO
              
              DO iband = 1,mnobd
                DO i = 1,3
                  DO ig = 1,atoms%jri(itype)
                    ! the r factor is already included in 
                    u1(ig,i,iband,lm2,iatom) = u1(ig,i,iband,lm2,iatom) - img*u1_lo(ig,ilo,itype)*carr2(i,iband)
                    u2(ig,i,iband,lm2,iatom) = u2(ig,i,iband,lm2,iatom) - img*u2_lo(ig,ilo,itype)*carr2(i,iband)
                  END DO
                END DO
              END DO  

            END DO
            
            l2  = l1 - 1
            IF( l2 .ge. 0 ) THEN
              lm2 = l2**2
              DO m2 = -l2,l2
                lm2   = lm2 + 1
                carr2 = 0
  
                DO m1 = -l1,l1
                  carr = gauntvec(l2,m2,l1,m1,atoms)
                  DO iband = 1,mnobd
                    carr2(1:3,iband) = carr2(1:3,iband) + cmt_lo(iband,m1,ilo,iatom)*carr
                  END DO
                END DO
              
                DO iband = 1,mnobd
                  DO i = 1,3
                    DO ig = 1,atoms%jri(itype)
                      ! the r factor is already included in 
                      u1(ig,i,iband,lm2,iatom)= u1(ig,i,iband,lm2,iatom) - img*u1_lo(ig,ilo,itype)*carr2(i,iband)
                      u2(ig,i,iband,lm2,iatom)= u2(ig,i,iband,lm2,iatom) - img*u2_lo(ig,ilo,itype)*carr2(i,iband)
                    END DO
                  END DO
                END DO  

              END DO
            END IF

          END DO
        END DO
      END DO
      
     

      !
      ! calculate projection < phi(n',k)|phi^1(n,k)>
      ! n' = iband1 , n= iband2
      !
      iatom     = 0
      proj_ibsc = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
          lm    = 0
          lmp   = 0
          DO l = 0,atoms%lmax(itype)
            DO M = -l,l
              lm = lm + 1
              ru1 = real( u1(:,:,:,lm,iatom) )
              iu1 = aimag( u1(:,:,:,lm,iatom) )
              ru2 = real( u2(:,:,:,lm,iatom) )
              iu2 = aimag( u2(:,:,:,lm,iatom) )
              DO p = 1,2
                lmp = lmp + 1
                
Daniel Wortmann's avatar
Daniel Wortmann committed
519
                DO iband = 1,mnobd! hybdat%nbands
520 521
                  DO i = 1,3
                    
Daniel Wortmann's avatar
Daniel Wortmann committed
522
                    rintegrand = atoms%rmsh(:,itype)*(hybdat%bas1(:,p,l,itype)*ru1(:,i,iband) +hybdat%bas2(:,p,l,itype)*ru2(:,i,iband))
523

Daniel Wortmann's avatar
Daniel Wortmann committed
524
                    iintegrand = atoms%rmsh(:,itype)*(hybdat%bas1(:,p,l,itype)*iu1(:,i,iband) +hybdat%bas2(:,p,l,itype)*iu2(:,i,iband))
525
                    
Daniel Wortmann's avatar
Daniel Wortmann committed
526 527
                    carr2(i,iband)=     intgrf(rintegrand,atoms%jri,atoms%jmtd,atoms%rmsh, atoms%dx,atoms%ntype,itype,hybdat%gridf)&
                                  + img*intgrf(iintegrand,atoms%jri,atoms%jmtd,atoms%rmsh, atoms%dx,atoms%ntype,itype,hybdat%gridf)
528 529 530 531
  
                  END DO
                END DO
                
Daniel Wortmann's avatar
Daniel Wortmann committed
532
                DO iband1 = 1,hybdat%nbands(nk)
533
                  cdum = conjg(cmt_apw(iband1,lmp,iatom))
Daniel Wortmann's avatar
Daniel Wortmann committed
534
                  DO iband2 = 1,mnobd! hybdat%nbands
535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559
                    proj_ibsc(1:3,iband2,iband1)= proj_ibsc(1:3,iband2,iband1)+cdum*carr2(1:3,iband2)
                  END DO
                END DO
                
              END DO!p
            END DO!M
          END DO!l
        
        END DO!ieq
      END DO!itype
      
      iatom = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
          DO ilo = 1,atoms%nlo(itype)
            l  = atoms%llo(ilo,itype)
            lm = l**2
            DO M = -l,l
              lm = lm + 1
              ru1 = real( u1(:,:,:,lm,iatom) )
              iu1 = aimag( u1(:,:,:,lm,iatom) )
              ru2 = real( u2(:,:,:,lm,iatom) )
              iu2 = aimag( u2(:,:,:,lm,iatom) )
              
Daniel Wortmann's avatar
Daniel Wortmann committed
560
              DO iband = 1,mnobd! hybdat%nbands
561 562 563 564 565 566
                DO i = 1,3
                    
                  rintegrand = atoms%rmsh(:,itype)*(u1_lo(:,ilo,itype)*ru1(:,i,iband) +u2_lo(:,ilo,itype)*ru2(:,i,iband) )

                  iintegrand = atoms%rmsh(:,itype)*(u1_lo(:,ilo,itype)*iu1(:,i,iband) +u2_lo(:,ilo,itype)*iu2(:,i,iband) )
                    
Daniel Wortmann's avatar
Daniel Wortmann committed
567 568
                  carr2(i,iband) =     intgrf(rintegrand,atoms%jri,atoms%jmtd,atoms%rmsh, atoms%dx,atoms%ntype,itype,hybdat%gridf)&
                                 + img*intgrf(iintegrand,atoms%jri,atoms%jmtd,atoms%rmsh, atoms%dx,atoms%ntype,itype,hybdat%gridf)
569 570 571 572
  
                END DO
              END DO
              
Daniel Wortmann's avatar
Daniel Wortmann committed
573
              DO iband1 = 1,hybdat%nbands(nk)
574
                cdum = conjg(cmt_lo(iband1,M,ilo,iatom))
Daniel Wortmann's avatar
Daniel Wortmann committed
575
                DO iband2 = 1,mnobd! hybdat%nbands
576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603
                  proj_ibsc(1:3,iband2,iband1) = proj_ibsc(1:3,iband2,iband1) + cdum*carr2(1:3,iband2)
                END DO
              END DO
              
            END DO
          
          END DO
        END DO
      END DO
      
      !
      ! calculate <phi^1(n1,k)|phi^1(n2,k)>
      ! n1 and n2 occupied
      !
      iatom     = 0
      olap_ibsc = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
          lm    = 0
          DO l = 0,atoms%lmax(itype)!+1
            DO M = -l,l
              lm = lm + 1
              ru1 = real( u1(:,:,:,lm,iatom) )
              iu1 = aimag( u1(:,:,:,lm,iatom) )
              ru2 = real( u2(:,:,:,lm,iatom) )
              iu2 = aimag( u2(:,:,:,lm,iatom) )
              
Daniel Wortmann's avatar
Daniel Wortmann committed
604
              DO iband1 = 1,mnobd ! hybdat%nbands
605 606 607 608 609 610 611 612 613 614 615
                DO iband2 = 1,mnobd!iband1 
                  DO i = 1,3
                    DO j = 1,3
                    
                      rintegrand = atoms%rmsh(:,itype)**2*(ru1(:,i,iband1)*ru1(:,j,iband2) + ru2(:,i,iband1)*ru2(:,j,iband2)&
                                    + iu1(:,i,iband1)*iu1(:,j,iband2) + iu2(:,i,iband1)*iu2(:,j,iband2) )

                      iintegrand = atoms%rmsh(:,itype)**2*(ru1(:,i,iband1)*iu1(:,j,iband2) + ru2(:,i,iband1)*iu2(:,j,iband2)&
                                    - iu1(:,i,iband1)*ru1(:,j,iband2) - iu2(:,i,iband1)*ru2(:,j,iband2) )
                  
                      olap_ibsc(i,j,iband2,iband1) = olap_ibsc(i,j,iband2,iband1) &
Daniel Wortmann's avatar
Daniel Wortmann committed
616 617
                    +     intgrf(rintegrand,atoms%jri,atoms%jmtd,atoms%rmsh, atoms%dx,atoms%ntype,itype,hybdat%gridf)&
                    + img*intgrf(iintegrand,atoms%jri,atoms%jmtd,atoms%rmsh, atoms%dx,atoms%ntype,itype,hybdat%gridf)
618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671
  
                    END DO
                  END DO
              
                END DO
              END DO
                          
            END DO
          END DO
        END DO
      END DO
      
      

      END SUBROUTINE  ibs_correction
      
      FUNCTION gauntvec(l1,m1,l2,m2,atoms)     

      USE m_constants  
      USE m_gaunt
      
      USE m_types
      IMPLICIT NONE
      TYPE(t_atoms),INTENT(IN)   :: atoms
      
      INTEGER , INTENT(IN)  ::  l1,m1,l2,m2
      
      COMPLEX               ::  gauntvec(-1:1)
      
      INTEGER               ::  j,mj
      INTEGER               ::  point(-1:1)
      REAL                  ::  rfac
      COMPLEX               ::  cfac
      COMPLEX , PARAMETER   ::  img = (0d0,1d0)
      
      rfac = sqrt(tpi_const/3)
      DO j = -1,1
        ! j = -1 corresponds to x-direction
        ! j = 0  corresponds to z-direction
        ! j = 1  corresponds to y-direction
        mj  = abs(j)
        cfac = img**( (j+mj)/2. ) * sqrt(2.)**(1-mj) * rfac
        gauntvec(j) = cfac* (   gaunt1(1,l1,l2,-mj,m1,m2,atoms%lmaxd+1) + j*gaunt1(1,l1,l2, mj,m1,m2,atoms%lmaxd+1) )
      END DO

      ! transform onto cartesian coordinates
      point(-1) = -1
      point( 0) =  1
      point( 1) =  0

      gauntvec = gauntvec(point)

      END FUNCTION gauntvec

Daniel Wortmann's avatar
Daniel Wortmann committed
672 673
      FUNCTION w(p1,l1,p2,l2,itype,bas1_mt,drbas1_mt,&
                 rmt)
674 675
        USE m_types
        IMPLICIT NONE
Daniel Wortmann's avatar
Daniel Wortmann committed
676 677 678 679

      INTEGER , INTENT(IN)       ::  p1,l1,p2,l2
      INTEGER , INTENT(IN)       ::  itype
      REAL,INTENT(IN)            :: rmt(:),bas1_mt(:,0:,:),drbas1_mt(:,0:,:)
680 681 682 683 684 685 686 687 688 689 690 691
      
      REAL                  ::  w
      
      INTEGER               ::  p
      REAL                  ::  denom,enum
      
      IF( p1 .gt. 2 .or. p2 .gt. 2 ) STOP 'w: the formalism is only valid for p<=2'
      
      denom = wronskian( bas1_MT(2,l1,itype),drbas1_MT(2,l1,itype), bas1_MT(1,l1,itype),drbas1_MT(1,l1,itype) )

      p = p1 + (-1)**(p1-1)
      
Daniel Wortmann's avatar
Daniel Wortmann committed
692
      enum  = bas1_MT(p,l1,itype)*bas1_MT(p2,l2,itype) + rmt(itype)*wronskian(  bas1_MT(p ,l1,itype),&
693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733
                                   drbas1_MT(p ,l1,itype), bas1_MT(p2,l2,itype), drbas1_MT(p2,l2,itype) )
      
      w = (-1)**(p1+1)*enum/denom
      
      END FUNCTION
      
      
      
      PURE FUNCTION wronskian(f,df,g,dg)
      IMPLICIT NONE
      REAL , INTENT(IN) ::  f,df,g,dg
      REAL              ::  wronskian
      
      wronskian = f*dg - df*g
           
      END FUNCTION

!     Calculates the derivative
!                                  ikr    s       s
!     dcprod(n',n,k,xyz) = d    < e    phi   | phi       > / sqrt(vol)
!                           xyz           qn      q+k,n'
!
!                                s             s
!                           < phi  | d    | phi    >
!                                nq   xyz      n'q
!                      = -i ------------------------ / sqrt(vol)  ,   s = ispin
!                                  s     s                        ,   n = occ.     ,   n' = unocc.
!                                 e    - e                        ,   bandi1 <= n <= bandf1 , bandi2 <= n' <= bandf2
!                                  n'q    nq
!
!     with kp perturbation theory and
!
!               d     d     d
!     d    =   ---,  ---,  --- .
!      xyz     dk    dk    dk
!                x     y     z
!
      SUBROUTINE dwavefproducts(&
                           dcprod,nk,bandi1,bandf1,bandi2,bandf2,lwrite,&
                           atoms,hybrid,&
                           cell,&
Daniel Wortmann's avatar
Daniel Wortmann committed
734 735 736
                           hybdat,kpts,nkpti,lapw,&
                           dimension,jsp,&
                           eig_irr)
737 738 739 740 741

      USE m_wrapper
      USE m_types
      IMPLICIT NONE

Daniel Wortmann's avatar
Daniel Wortmann committed
742 743
      TYPE(t_hybdat),INTENT(IN)   :: hybdat

744 745 746 747 748 749 750 751 752
      TYPE(t_dimension),INTENT(IN)   :: dimension
      TYPE(t_hybrid),INTENT(IN)   :: hybrid
      TYPE(t_cell),INTENT(IN)   :: cell
      TYPE(t_kpts),INTENT(IN)   :: kpts
      TYPE(t_atoms),INTENT(IN)   :: atoms
      TYPE(t_lapw),INTENT(IN)   :: lapw

!     - scalars -
      INTEGER,INTENT(IN)      ::  nk,bandi1,bandf1,bandi2,bandf2
Daniel Wortmann's avatar
Daniel Wortmann committed
753
      INTEGER,INTENT(IN)      :: nkpti
754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774
      INTEGER,INTENT(IN)      :: jsp   


!     - arrays -

      REAL,INTENT(IN)         ::  eig_irr(dimension%neigd,nkpti)
      COMPLEX,INTENT(OUT)     ::  dcprod(bandi2:bandf2,bandi1:bandf1,3)

!     - local scalars -
      INTEGER                 ::  ikpt,ikpt1,iband1,iband2
      REAL                    ::  rdum
      LOGICAL                 ::  lwrite


      !                                       __
      ! Get momentum-matrix elements -i < uj | \/ | ui >
      !
      CALL momentum_matrix(&
                     dcprod,nk,bandi1,bandf1,bandi2,bandf2,&
                     atoms,hybrid,&
                     cell,&
Daniel Wortmann's avatar
Daniel Wortmann committed
775 776
                     hybdat,kpts,lapw,&
                     dimension,jsp)
777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812


      !                                                __
      !  Calculate expansion coefficients -i < uj | \/ | ui > / ( ei - ej ) for periodic function ui 
      ! 
      DO iband1 = bandi1,bandf1
        DO iband2 = bandi2,bandf2
          rdum = eig_irr(iband2,nk) - eig_irr(iband1,nk)
          IF(abs(rdum).gt.1d-6) THEN !1d-6
            dcprod(iband2,iband1,:) = dcprod(iband2,iband1,:) / rdum
          ELSE
            dcprod(iband2,iband1,:) = 0d0
          END IF
        END DO
      END DO
  
      dcprod = dcprod / sqrt(cell%omtil)

     

      END SUBROUTINE dwavefproducts

!     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

!     Calculates the momentum matrix elements
!
!                                    s             s
!     momentum(n',n,q,xyz) = -i < phi  | d    | phi    >   ,   s  = ispin
!                                    nq   xyz      n'q     ,   n  = occ.    ( bandi1 <= n  <= bandf1 )  ,
!                                                              n' = unocc.  ( bandi2 <= n' <= bandf2 )
!
!
      SUBROUTINE momentum_matrix(&
                      momentum,nk,bandi1,bandf1,bandi2,bandf2,&
                      atoms,hybrid,&
                      cell,&
Daniel Wortmann's avatar
Daniel Wortmann committed
813 814
                      hybdat,kpts,lapw,&
                      dimension,jsp)
815 816 817 818 819 820 821 822 823



      USE m_olap
      USE m_wrapper
      USE m_util    ,only: derivative,intgrf_init,intgrf
      USE m_dr2fdr
      USE m_constants
      USE m_types
Daniel Wortmann's avatar
Daniel Wortmann committed
824
      USE m_io_hybrid
825
      IMPLICIT NONE
Daniel Wortmann's avatar
Daniel Wortmann committed
826 827

      TYPE(t_hybdat),INTENT(IN)   :: hybdat
828 829 830 831 832 833 834 835 836 837 838 839 840 841
      TYPE(t_dimension),INTENT(IN)   :: dimension
      TYPE(t_hybrid),INTENT(IN)   :: hybrid
      TYPE(t_cell),INTENT(IN)   :: cell
      TYPE(t_kpts),INTENT(IN)   :: kpts
      TYPE(t_atoms),INTENT(IN)   :: atoms
      TYPE(t_lapw),INTENT(IN)   :: lapw

!     - scalars -
      INTEGER,INTENT(IN)      ::  bandi1,bandf1,bandi2,bandf2
      INTEGER,INTENT(IN)      :: nk 
      INTEGER,INTENT(IN)      :: jsp


!     - arrays -
Daniel Wortmann's avatar
Daniel Wortmann committed
842
      TYPE(t_mat):: z
843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864
      COMPLEX,INTENT(OUT)     :: momentum(bandi2:bandf2,bandi1:bandf1,3)

!     - local scalars -
      INTEGER                 ::  itype,ieq,ic,i,j,l ,lm,n1,n2,ikpt, iband1,iband2,ll,mm
      INTEGER                 ::  lm_0,lm_1,lm0,lm1,lm2,lm3,n0,nn,n, l1,l2,m1,m2,ikpt1
      INTEGER                 ::  irecl_cmt,irecl_z,m
      COMPLEX                 ::  cdum
      COMPLEX                 ::  img=(0d0,1d0)

!     - local arrays -
      INTEGER                 ::  gpt(3,lapw%nv(jsp))

      REAL                    ::  fcoeff((atoms%lmaxd+1)**2,-1:1), gcoeff((atoms%lmaxd+1)**2,-1:1)
      REAL                    ::  qmat1(hybrid%maxindx,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype), dbas1(atoms%jmtd)
      REAL                    ::  qmat2(hybrid%maxindx,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype), dbas2(atoms%jmtd)
      REAL                    ::  qg(lapw%nv(jsp),3)

      COMPLEX                 ::  hlp(3,3)
      COMPLEX                 ::  cvec1(hybrid%maxlmindx),cvec2(hybrid%maxlmindx), cvec3(hybrid%maxlmindx)
      COMPLEX                 ::  cmt1(hybrid%maxlmindx,bandi1:bandf1), cmt2(hybrid%maxlmindx,bandi2:bandf2)
      COMPLEX                 ::  carr1(3),carr2(3)
      COMPLEX                 ::  cmt(dimension%neigd,hybrid%maxlmindx,atoms%nat)
Daniel Wortmann's avatar
Daniel Wortmann committed
865 866 867 868
      REAL                    ::  olap_r(lapw%nv(jsp)*(lapw%nv(jsp)+1)/2)
      COMPLEX                 ::  olap_c(lapw%nv(jsp)*(lapw%nv(jsp)+1)/2)
      REAL                    ::  vec1_r(lapw%nv(jsp)),vec2_r(lapw%nv(jsp)), vec3_r(lapw%nv(jsp))
      COMPLEX                 ::  vec1_c(lapw%nv(jsp)),vec2_c(lapw%nv(jsp)), vec3_c(lapw%nv(jsp))
869 870

      ! read in cmt coefficients from direct access file cmt at kpoint nk
Daniel Wortmann's avatar
Daniel Wortmann committed
871 872 873
     
      call read_cmt(cmt,nk)
      
874 875
      ! read in z coefficients from direct access file z at kpoint nk

Daniel Wortmann's avatar
Daniel Wortmann committed
876
      call read_z(z,nk)
877

Daniel Wortmann's avatar
Daniel Wortmann committed
878
      !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907

      gpt(1,1:lapw%nv(jsp)) = lapw%k1(1:lapw%nv(jsp),jsp)
      gpt(2,1:lapw%nv(jsp)) = lapw%k2(1:lapw%nv(jsp),jsp)
      gpt(3,1:lapw%nv(jsp)) = lapw%k3(1:lapw%nv(jsp),jsp)


!     Define coefficients F and G
      lm = 0
      DO l = 0,atoms%lmaxd
        DO M = -l,l
          lm = lm + 1
          fcoeff(lm,-1) = - sqrt( 1d0 * (l+M+1)*(l+M+2) / (2*(2*l+1)*(2*l+3)) )
          fcoeff(lm, 0) =   sqrt( 1d0 * (l-M+1)*(l+M+1) /   ((2*l+1)*(2*l+3)) )
          fcoeff(lm, 1) = - sqrt( 1d0 * (l-M+1)*(l-M+2) / (2*(2*l+1)*(2*l+3)) )
          gcoeff(lm,-1) =   sqrt( 1d0 * (l-M)*(l-M-1) / (2*(2*l-1)*(2*l+1)) )
          gcoeff(lm, 0) =   sqrt( 1d0 * (l-M)*(l+M) /   ((2*l-1)*(2*l+1)) )
          gcoeff(lm, 1) =   sqrt( 1d0 * (l+M)*(l+M-1) / (2*(2*l-1)*(2*l+1)) )
        END DO
      END DO


!     Calculate olap int r**2*u*u' + w * int r*u*u, w = -l,l+1 ( -> qmat1/2 )
      qmat1 = 0
      qmat2 = 0
      ic   = 0
      DO itype = 1,atoms%ntype
        DO l = 0,atoms%lmax(itype)
          DO n2 = 1,hybrid%nindx(l,itype)
            !ic = ic + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
908 909
            CALL derivative(dbas1,hybdat%bas1(:,n2,l,itype),atoms, itype)
            dbas1 = dbas1 - hybdat%bas1(:,n2,l,itype)/atoms%rmsh(:,itype)
910

Daniel Wortmann's avatar
Daniel Wortmann committed
911 912
            CALL derivative(dbas2,hybdat%bas2(:,n2,l,itype),atoms, itype)
            dbas2 = dbas2 - hybdat%bas2(:,n2,l,itype)/atoms%rmsh(:,itype)
913 914 915 916
            
            IF(l.ne.0) THEN
              DO n1 = 1,hybrid%nindx(l-1,itype)
                ic = ic + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
917 918 919 920
                qmat1(n1,n2,l,itype) = intgrf( dbas1(:) * hybdat%bas1(:,n1,l-1,itype) +&
                          dbas2(:) * hybdat%bas2(:,n1,l-1,itype) , atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,hybdat%gridf )&
                + intgrf( ( hybdat%bas1(:,n2,l,itype) * hybdat%bas1(:,n1,l-1,itype) + hybdat%bas2(:,n2,l,itype) * hybdat%bas2(:,n1,l-1,itype) ) &
                          / atoms%rmsh(:,itype) , atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,hybdat%gridf ) * (l+1)
921 922 923 924 925 926 927

              END DO
            END IF
            IF(l.ne.atoms%lmax(itype)) THEN
              DO n1 = 1,hybrid%nindx(l+1,itype)


Daniel Wortmann's avatar
Daniel Wortmann committed
928 929 930 931
                qmat2(n1,n2,l,itype) = intgrf( dbas1(:) * hybdat%bas1(:,n1,l+1,itype) + dbas2(:) * hybdat%bas2(:,n1,l+1,itype) , &
                          atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,hybdat%gridf )&
                - intgrf( ( hybdat%bas1(:,n2,l,itype) * hybdat%bas1(:,n1,l+1,itype) + hybdat%bas2(:,n2,l,itype) * hybdat%bas2(:,n1,l+1,itype) ) &
                          / atoms%rmsh(:,itype) , atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,hybdat%gridf ) * l
932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030


              END DO
            END IF


            
          END DO
        END DO
      END DO

      !                                                  __
      ! Calculate momentum matrix elements -i < uj | \/ | ui > wrt wave functions u (->momentum)
      !

      momentum = 0

      ! MT contribution

      ic = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          ic = ic + 1
          nn = sum( (/ ((2*l+1)*hybrid%nindx(l,itype),l=0,atoms%lmax(itype)) /) )
          DO iband1 = bandi1,bandf1 
            cmt1(:nn,iband1) = cmt(iband1,:nn,ic) 
          ENDDO
          DO iband2 = bandi2,bandf2 
            cmt2(:nn,iband2) = cmt(iband2,:nn,ic) 
          ENDDO
          DO iband1 = bandi1,bandf1

            cvec1 = 0 ; cvec2 = 0 ; cvec3 = 0
            ! build up left vector(s) ( -> cvec1/2/3 )
            lm_0 = 0              ! we start with s-functions (l=0)
            lm_1 = hybrid%nindx(0,itype) ! we start with p-functions (l=0+1)
            lm   = 0
            DO l = 0,atoms%lmax(itype)-1              
              n0 = hybrid%nindx(l,  itype)
              n1 = hybrid%nindx(l+1,itype)
              DO M = -l,l
                lm  = lm + 1
                lm0 = lm_0 + (M+l)    *n0
                lm1 = lm_1 + (M+1+l+1)*n1
                lm2 = lm_1 + (M  +l+1)*n1
                lm3 = lm_1 + (M-1+l+1)*n1
                cvec1(lm0+1:lm0+n0) = fcoeff(lm,-1) * matmul ( cmt1(lm1+1:lm1+n1,iband1) , qmat2(:n1,:n0,l,itype) )
                cvec2(lm0+1:lm0+n0) = fcoeff(lm, 0) * matmul ( cmt1(lm2+1:lm2+n1,iband1) , qmat2(:n1,:n0,l,itype) )
                cvec3(lm0+1:lm0+n0) = fcoeff(lm, 1) * matmul ( cmt1(lm3+1:lm3+n1,iband1) , qmat2(:n1,:n0,l,itype) )
              END DO
              lm_0 = lm_0 + (2*l+1)*n0
              lm_1 = lm_1 + (2*l+3)*n1
            END DO

            lm_0 = hybrid%nindx(0,itype) ! we start with p-functions (l=1)
            lm_1 = 0              ! we start with s-functions (l=1-1)
            lm   = 1
            DO l = 1,atoms%lmax(itype)
              n0 = hybrid%nindx(l,  itype)
              n1 = hybrid%nindx(l-1,itype)
              DO M = -l,l
                lm  = lm + 1
                lm0 = lm_0 + (M+l)    *n0
                lm1 = lm_1 + (M+1+l-1)*n1
                lm2 = lm_1 + (M  +l-1)*n1
                lm3 = lm_1 + (M-1+l-1)*n1
                IF(abs(M+1).le.l-1) cvec1(lm0+1:lm0+n0) = cvec1(lm0+1:lm0+n0) + gcoeff(lm,-1) * matmul ( cmt1(lm1+1:lm1+n1,iband1) , qmat1(:n1,:n0,l,itype) )
                IF(abs(M)  .le.l-1) cvec2(lm0+1:lm0+n0) = cvec2(lm0+1:lm0+n0) + gcoeff(lm, 0) * matmul ( cmt1(lm2+1:lm2+n1,iband1) , qmat1(:n1,:n0,l,itype) )
                IF(abs(M-1).le.l-1) cvec3(lm0+1:lm0+n0) = cvec3(lm0+1:lm0+n0) + gcoeff(lm, 1) * matmul ( cmt1(lm3+1:lm3+n1,iband1) , qmat1(:n1,:n0,l,itype) )
              END DO
              lm_0 = lm_0 + (2*l+1)*n0
              lm_1 = lm_1 + (2*l-1)*n1                
            END DO
            ! multiply with right vector
            DO iband2 = bandi2,bandf2
              momentum(iband2,iband1,1) = momentum(iband2,iband1,1) + dotprod(cvec1(:nn),cmt2(:nn,iband2))
              momentum(iband2,iband1,2) = momentum(iband2,iband1,2) + dotprod(cvec2(:nn),cmt2(:nn,iband2))
              momentum(iband2,iband1,3) = momentum(iband2,iband1,3) + dotprod(cvec3(:nn),cmt2(:nn,iband2))
            END DO ! iband2
          END DO ! iband1

        END DO ! ieq
      END DO ! itype

      ! Transform to cartesian coordinates
      hlp        =  0
      hlp( 1, 1) =  1d0/sqrt(2d0)
      hlp( 1, 3) = -1d0/sqrt(2d0)
      hlp( 2, 1) = -img/sqrt(2d0) 
      hlp( 2, 3) = -img/sqrt(2d0)
      hlp( 3, 2) =  1d0
      DO iband1 = bandi1,bandf1
        DO iband2 = bandi2,bandf2
          momentum(iband2,iband1,:) = -img*matmul(momentum(iband2,iband1,:),transpose(hlp))
        END DO
      END DO

      ! plane-wave contribution

Daniel Wortmann's avatar
Daniel Wortmann committed
1031
      CALL olap_pwp(z%l_real,olap_r,olap_c,gpt,lapw%nv(jsp),atoms,cell)
1032 1033 1034 1035 1036

      DO nn = 1,lapw%nv(jsp)
        qg(nn,:) = matmul(kpts%bk(:,nk)+gpt(:,nn),cell%bmat)
      END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
1037
      if (z%l_real) THEN
1038
      DO iband2 = bandi2,bandf2
Daniel Wortmann's avatar
Daniel Wortmann committed
1039 1040 1041
        vec1_r = matvec( olap_r , z%data_r(:lapw%nv(jsp),iband2) * qg(:,1) )
        vec2_r = matvec( olap_r , z%data_r(:lapw%nv(jsp),iband2) * qg(:,2) )
        vec3_r = matvec( olap_r , z%data_r(:lapw%nv(jsp),iband2) * qg(:,3) )
1042
        DO iband1 = bandi1,bandf1
Daniel Wortmann's avatar
Daniel Wortmann committed
1043 1044 1045
          momentum(iband2,iband1,1) = momentum(iband2,iband1,1) + dotprod(z%data_r(:lapw%nv(jsp),iband1),vec1_r)
          momentum(iband2,iband1,2) = momentum(iband2,iband1,2) + dotprod(z%data_r(:lapw%nv(jsp),iband1),vec2_r)
          momentum(iband2,iband1,3) = momentum(iband2,iband1,3) + dotprod(z%data_r(:lapw%nv(jsp),iband1),vec3_r)
1046 1047
        END DO
      END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060
   else
  DO iband2 = bandi2,bandf2
        vec1_c = matvec( olap_c , z%data_c(:lapw%nv(jsp),iband2) * qg(:,1) )
        vec2_c = matvec( olap_c , z%data_c(:lapw%nv(jsp),iband2) * qg(:,2) )
        vec3_c = matvec( olap_c , z%data_c(:lapw%nv(jsp),iband2) * qg(:,3) )
        DO iband1 = bandi1,bandf1
          momentum(iband2,iband1,1) = momentum(iband2,iband1,1) + dotprod(z%data_c(:lapw%nv(jsp),iband1),vec1_c)
          momentum(iband2,iband1,2) = momentum(iband2,iband1,2) + dotprod(z%data_c(:lapw%nv(jsp),iband1),vec2_c)
          momentum(iband2,iband1,3) = momentum(iband2,iband1,3) + dotprod(z%data_c(:lapw%nv(jsp),iband1),vec3_c)
        END DO
     END DO
  end if
      
1061 1062 1063
      END SUBROUTINE momentum_matrix
      
      END MODULE m_kp_perturbation