LCOV - code coverage report
Current view: top level - hybrid - kp_perturbation.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 468 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 6 0.0 %

          Line data    Source code
       1             : MODULE m_kp_perturbation
       2             : 
       3             :       CONTAINS
       4             : 
       5           0 :          SUBROUTINE ibs_correction( &
       6             :             nk, atoms, &
       7             :             dimension, input, jsp, &
       8             :             hybdat, hybrid, &
       9             :             lapw, kpts, nkpti, &
      10             :             cell, mnobd, &
      11             :             sym, &
      12           0 :             proj_ibsc, olap_ibsc)
      13             : 
      14             :             USE m_sphbes
      15             :             USE m_dsphbs
      16             :             USE m_constants
      17             :             USE m_ylm
      18             :             USE m_gaunt
      19             :             USE m_util
      20             :             USE m_types
      21             :             USE m_io_hybrid
      22             :             IMPLICIT NONE
      23             : 
      24             :             TYPE(t_hybdat), INTENT(IN)   :: hybdat
      25             :             TYPE(t_dimension), INTENT(IN)   :: dimension
      26             :             TYPE(t_hybrid), INTENT(INOUT)   :: hybrid
      27             :             TYPE(t_input), INTENT(IN)   :: input
      28             :             TYPE(t_sym), INTENT(IN)   :: sym
      29             :             TYPE(t_cell), INTENT(IN)   :: cell
      30             :             TYPE(t_kpts), INTENT(IN)   :: kpts
      31             :             TYPE(t_atoms), INTENT(IN)   :: atoms
      32             :             TYPE(t_lapw), INTENT(IN)   :: lapw
      33             : 
      34             :             ! - scalars -
      35             :             INTEGER, INTENT(IN)   :: jsp
      36             :             INTEGER, INTENT(IN)   ::  mnobd
      37             :             INTEGER, INTENT(IN)   ::  nk, nkpti
      38             : 
      39             : !       REAL   , INTENT(INOUT)::  ibs_corr(3,3,nbands)
      40             : 
      41             :             ! - arrays -
      42             : 
      43             :             COMPLEX, INTENT(INOUT)::  olap_ibsc(3, 3, mnobd, mnobd)
      44             :             COMPLEX, INTENT(INOUT)::  proj_ibsc(:, :, :)!(3,mnobd,hybrid%nbands(nk))
      45             :             ! - local scalars -
      46             :             INTEGER               ::  i, itype, ieq, iatom, iatom1, iband, iband1
      47             :             INTEGER               ::  iband2, ilo, ibas, ic, ikpt, ikvec, invsfct
      48             :             INTEGER               ::  irecl_cmt, irecl_z
      49             :             INTEGER               ::  j, m
      50             :             INTEGER               ::  l1, m1, p1, l2, m2, p2, l, p, lm, &
      51             :                                      lmp, lmp1, lmp2, lm1, lm2
      52             :             INTEGER               ::  ok, ig
      53             :             INTEGER               ::  idum
      54             :             REAL                  ::  const
      55             :             REAL                  ::  ka, kb
      56             :             REAL                  ::  kvecn
      57             :             REAL                  ::  olap_udot, olap_uulo, olap_udotulo
      58             :             REAL                  ::  rdum
      59             :             REAL                  ::  ws
      60             :             COMPLEX               ::  phase
      61             :             COMPLEX               ::  cj, cdj
      62             :             COMPLEX               ::  denom, enum
      63             :             COMPLEX               ::  cdum, cdum1, cdum2
      64             :             COMPLEX, PARAMETER    ::  img = (0.0, 1.0)
      65             :             ! - local arrays -
      66           0 :             INTEGER               ::  lmp_start(atoms%ntype)
      67           0 :             REAL                  ::  alo(atoms%nlod, atoms%ntype), blo(atoms%nlod, atoms%ntype), &
      68           0 :                                      clo(atoms%nlod, atoms%ntype)
      69           0 :             REAL                  ::  u1_lo(atoms%jmtd, atoms%nlod, atoms%ntype), &
      70           0 :                                      u2_lo(atoms%jmtd, atoms%nlod, atoms%ntype)
      71             :             REAL                  ::  kvec(3), qvec(3)
      72           0 :             REAL                  ::  sbes(0:atoms%lmaxd + 1), dsbes(0:atoms%lmaxd + 1)
      73           0 :             REAL                  ::  bas1_tmp(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd + 1, atoms%ntype), &
      74           0 :                                      bas2_tmp(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd + 1, atoms%ntype)
      75           0 :             REAL                  ::  bas1_MT_tmp(hybrid%maxindx, 0:atoms%lmaxd + 1, atoms%ntype), &
      76           0 :                                      drbas1_MT_tmp(hybrid%maxindx, 0:atoms%lmaxd + 1, atoms%ntype)
      77           0 :             REAL                  ::  ru1(atoms%jmtd, 3, mnobd), ru2(atoms%jmtd, 3, mnobd)
      78           0 :             REAL                  ::  iu1(atoms%jmtd, 3, mnobd), iu2(atoms%jmtd, 3, mnobd)
      79           0 :             REAL                  ::  rintegrand(atoms%jmtd), iintegrand(atoms%jmtd), &
      80           0 :                                      integrand(atoms%jmtd)
      81             : 
      82             :             COMPLEX               ::  f(atoms%jmtd, mnobd)
      83           0 :             COMPLEX               ::  carr(3), carr2(3, hybrid%nbands(nk))
      84           0 :             COMPLEX               ::  ylm((atoms%lmaxd + 2)**2)
      85           0 :             COMPLEX, ALLOCATABLE   ::  u1(:, :, :, :, :), u2(:, :, :, :, :)
      86           0 :             COMPLEX, ALLOCATABLE   ::  cmt_lo(:, :, :, :)
      87           0 :             COMPLEX, ALLOCATABLE   ::  cmt_apw(:, :, :)
      88           0 :             TYPE(t_mat)           ::  z
      89           0 :             REAL                  ::  work_r(dimension%neigd)
      90           0 :             COMPLEX               ::  work_c(dimension%neigd)
      91             : 
      92             :             !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
      93             : 
      94           0 :             bas1_tmp(:, :, 0:atoms%lmaxd, :) = hybdat%bas1(:, :, 0:atoms%lmaxd, :)
      95           0 :             bas2_tmp(:, :, 0:atoms%lmaxd, :) = hybdat%bas2(:, :, 0:atoms%lmaxd, :)
      96             : 
      97           0 :             bas1_MT_tmp(:, 0:atoms%lmaxd, :) = hybdat%bas1_MT(:, 0:atoms%lmaxd, :)
      98           0 :             drbas1_MT_tmp(:, 0:atoms%lmaxd, :) = hybdat%drbas1_MT(:, 0:atoms%lmaxd, :)
      99             : 
     100           0 :             bas1_tmp(:, :, atoms%lmaxd + 1, :) = hybdat%bas1(:, :, atoms%lmaxd, :)
     101           0 :             bas2_tmp(:, :, atoms%lmaxd + 1, :) = hybdat%bas2(:, :, atoms%lmaxd, :)
     102             : 
     103           0 :             bas1_MT_tmp(:, atoms%lmaxd + 1, :) = hybdat%bas1_MT(:, atoms%lmaxd, :)
     104           0 :             drbas1_MT_tmp(:, atoms%lmaxd + 1, :) = hybdat%drbas1_MT(:, atoms%lmaxd, :)
     105             : 
     106             :             ! read in z coefficient from direct access file z at k-point nk
     107             : 
     108           0 :             call read_z(z, nk)
     109             : 
     110             :             ! construct local orbital consisting of radial function times spherical harmonic
     111             :             ! where the radial function vanishes on the MT sphere boundary
     112             :             ! with this the local orbitals have a trivial k-dependence
     113             : 
     114             :             ! compute radial lo matching coefficients
     115           0 :             hybrid%nindx = 2
     116           0 :             DO itype = 1, atoms%ntype
     117           0 :                DO ilo = 1, atoms%nlo(itype)
     118           0 :                   l = atoms%llo(ilo, itype)
     119           0 :                   hybrid%nindx(l, itype) = hybrid%nindx(l, itype) + 1
     120           0 :                   p = hybrid%nindx(l, itype)
     121             : 
     122           0 :                   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))
     123             : 
     124           0 :                   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))
     125             : 
     126           0 :                   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))
     127             : 
     128           0 :                   integrand = hybdat%bas1(:, 2, l, itype)*hybdat%bas1(:, 2, l, itype) + hybdat%bas2(:, 2, l, itype)*hybdat%bas2(:, 2, l, itype)
     129           0 :                   olap_udot = intgrf(integrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
     130             : 
     131           0 :                   integrand = hybdat%bas1(:, 1, l, itype)*hybdat%bas1(:, p, l, itype) + hybdat%bas2(:, 1, l, itype)*hybdat%bas2(:, p, l, itype)
     132           0 :                   olap_uulo = intgrf(integrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
     133             : 
     134           0 :                   integrand = hybdat%bas1(:, 2, l, itype)*hybdat%bas1(:, p, l, itype) + hybdat%bas2(:, 2, l, itype)*hybdat%bas2(:, p, l, itype)
     135           0 :                   olap_udotulo = intgrf(integrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
     136             : 
     137           0 :                   rdum = ka**2 + (kb**2)*olap_udot + 1.0 + 2.0*ka*olap_uulo + 2.0*kb*olap_udotulo
     138           0 :                   clo(ilo, itype) = 1.0/sqrt(rdum)
     139           0 :                   alo(ilo, itype) = ka*clo(ilo, itype)
     140           0 :                   blo(ilo, itype) = kb*clo(ilo, itype)
     141             : 
     142           0 :                   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)
     143             : 
     144           0 :                   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)
     145             :                END DO
     146             :             END DO
     147             : 
     148             :             ! calculate lo wavefunction coefficients
     149           0 :             ALLOCATE (cmt_lo(dimension%neigd, -atoms%llod:atoms%llod, atoms%nlod, atoms%nat))
     150           0 :             cmt_lo = 0
     151           0 :             iatom = 0
     152           0 :             ic = 0
     153           0 :             ibas = lapw%nv(jsp)
     154           0 :             DO itype = 1, atoms%ntype
     155             :                ! the program is in hartree units, therefore 1/wronskian is
     156             :                ! (rmt**2)/2.
     157           0 :                const = fpi_const*(atoms%rmt(itype)**2)/2/sqrt(cell%omtil)
     158           0 :                DO ieq = 1, atoms%neq(itype)
     159           0 :                   iatom = iatom + 1
     160           0 :                   IF ((atoms%invsat(iatom) == 0) .or. (atoms%invsat(iatom) == 1)) THEN
     161           0 :                      IF (atoms%invsat(iatom) == 0) invsfct = 1
     162           0 :                      IF (atoms%invsat(iatom) == 1) THEN
     163           0 :                         invsfct = 2
     164           0 :                         iatom1 = sym%invsatnr(iatom)
     165             :                      END IF
     166             : 
     167           0 :                      DO ilo = 1, atoms%nlo(itype)
     168           0 :                         l = atoms%llo(ilo, itype)
     169           0 :                         cdum = img**l*const
     170           0 :                         DO ikvec = 1, invsfct*(2*l + 1)
     171           0 :                            ic = ic + 1
     172           0 :                            ibas = ibas + 1
     173           0 :                            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)/)
     174             : 
     175           0 :                            phase = exp(img*tpi_const*dot_product(atoms%taual(:, iatom), kvec))
     176           0 :                            cdum1 = cdum*phase
     177             : 
     178           0 :                            CALL ylm4(l, matmul(kvec, cell%bmat), ylm)
     179             : 
     180           0 :                            lm = l**2
     181           0 :                            DO M = -l, l
     182           0 :                               lm = lm + 1
     183           0 :                               cdum2 = cdum1*conjg(ylm(lm))
     184           0 :                               if (z%l_real) THEN
     185           0 :                                  work_r = z%data_r(ibas, :)
     186           0 :                                  DO iband = 1, hybrid%nbands(nk)
     187           0 :                                     cmt_lo(iband, M, ilo, iatom) = cmt_lo(iband, M, ilo, iatom) + cdum2*work_r(iband)
     188           0 :                                     IF (invsfct == 2) THEN
     189             :                                        ! the factor (-1)**l is necessary as we do not calculate
     190             :                                        ! the cmt_lo in the local coordinate system of the atom
     191           0 :                                        cmt_lo(iband, -M, ilo, iatom1) = cmt_lo(iband, -M, ilo, iatom1) + (-1)**(l + M)*conjg(cdum2)*work_r(iband)
     192             :                                     END IF
     193             :                                  END DO
     194             :                               else
     195           0 :                                  work_c = z%data_c(ibas, :)
     196           0 :                                  DO iband = 1, hybrid%nbands(nk)
     197           0 :                                     cmt_lo(iband, M, ilo, iatom) = cmt_lo(iband, M, ilo, iatom) + cdum2*work_c(iband)
     198           0 :                                     IF (invsfct == 2) THEN
     199             :                                        ! the factor (-1)**l is necessary as we do not calculate
     200             :                                        ! the cmt_lo in the local coordinate system of the atom
     201           0 :                                        cmt_lo(iband, -M, ilo, iatom1) = cmt_lo(iband, -M, ilo, iatom1) + (-1)**(l + M)*conjg(cdum2)*work_c(iband)
     202             :                                     END IF
     203             :                                  END DO
     204             :                               end if
     205             : 
     206             :                            END DO
     207             : 
     208             :                         END DO  !ikvec
     209             :                      END DO  ! ilo
     210             : 
     211             :                   END IF
     212             : 
     213             :                END DO  !ieq
     214             :             END DO  !itype
     215             : 
     216             :             !
     217             :             ! calculate apw wavefunction coefficients up to lmax + 1
     218             :             ! note that the lo contribution is separated in cmt_lo
     219             :             !
     220             : 
     221           0 :             DO itype = 1, atoms%ntype
     222           0 :                lmp_start(itype) = sum((/(2*(2*l + 1), l=0, atoms%lmax(itype) + 1)/))
     223             :             END DO
     224           0 :             idum = maxval(lmp_start)
     225             : 
     226           0 :             ALLOCATE (cmt_apw(dimension%neigd, idum, atoms%nat))
     227           0 :             cmt_apw = 0
     228           0 :             DO i = 1, lapw%nv(jsp)
     229           0 :                kvec = kpts%bk(:, nk) + (/lapw%k1(i, jsp), lapw%k2(i, jsp), lapw%k3(i, jsp)/)
     230           0 :                kvecn = sqrt(dot_product(matmul(kvec, cell%bmat), matmul(kvec, cell%bmat)))
     231             : 
     232           0 :                iatom = 0
     233           0 :                DO itype = 1, atoms%ntype
     234             :                   !calculate spherical sperical harmonics
     235           0 :                   CALL ylm4(atoms%lmax(itype) + 1, matmul(kvec, cell%bmat), ylm)
     236             : 
     237             :                   !calculate spherical bessel function at |kvec|*R_MT(itype)
     238           0 :                   CALL sphbes(atoms%lmax(itype) + 1, kvecn*atoms%rmt(itype), sbes)
     239             : 
     240             :                   !calculate radial derivative of spherical bessel function at |kvec|*R_MT(itype)
     241           0 :                   CALL dsphbs(atoms%lmax(itype) + 1, kvecn*atoms%rmt(itype), sbes, dsbes)
     242           0 :                   dsbes = kvecn*dsbes
     243             : 
     244           0 :                   DO ieq = 1, atoms%neq(itype)
     245           0 :                      iatom = iatom + 1
     246             : 
     247           0 :                      phase = exp(img*tpi_const*dot_product(kvec, atoms%taual(:, iatom)))
     248             : 
     249           0 :                      lm = 0
     250           0 :                      lmp = 0
     251           0 :                      DO l = 0, atoms%lmax(itype) + 1
     252             :                         denom = wronskian(bas1_MT_tmp(2, l, itype), drbas1_MT_tmp(2, l, itype), &
     253           0 :                                           bas1_MT_tmp(1, l, itype), drbas1_MT_tmp(1, l, itype))
     254           0 :                         cdum1 = fpi_const*img**l*sbes(l)*phase/sqrt(cell%omtil)
     255           0 :                         cdum2 = fpi_const*img**l*dsbes(l)*phase/sqrt(cell%omtil)
     256           0 :                         DO M = -l, l
     257           0 :                            lm = lm + 1
     258           0 :                            cj = cdum1*conjg(ylm(lm))
     259           0 :                            cdj = cdum2*conjg(ylm(lm))
     260           0 :                            DO p = 1, 2
     261           0 :                               lmp = lmp + 1
     262           0 :                               p1 = p + (-1)**(p - 1)
     263             : 
     264             :                               enum = CMPLX(wronskian(bas1_MT_tmp(p1, l, itype), drbas1_MT_tmp(p1, l, itype), REAL(cj), REAL(cdj)), &
     265           0 :                                            wronskian(bas1_MT_tmp(p1, l, itype), drbas1_MT_tmp(p1, l, itype), AIMAG(cj), AIMAG(cdj)))
     266             : 
     267           0 :                               cdum = (-1)**(p + 1)*enum/denom
     268           0 :                               if (z%l_real) THEN
     269           0 :                                  work_r = z%data_r(i, :)
     270           0 :                                  DO iband = 1, hybrid%nbands(nk)
     271           0 :                                     cmt_apw(iband, lmp, iatom) = cmt_apw(iband, lmp, iatom) + cdum*work_r(iband)
     272             :                                  END DO
     273             :                               else
     274           0 :                                  work_c = z%data_c(i, :)
     275           0 :                                  DO iband = 1, hybrid%nbands(nk)
     276           0 :                                     cmt_apw(iband, lmp, iatom) = cmt_apw(iband, lmp, iatom) + cdum*work_c(iband)
     277             :                                  END DO
     278             :                               end if
     279             :                            END DO  !p
     280             :                         END DO  !M
     281             :                      END DO  !l
     282             : 
     283             :                   END DO  !iatom
     284             :                END DO  ! itype
     285             :             END DO  ! i
     286             : 
     287             :             ! construct radial functions (complex) for the first order
     288             :             ! incomplete basis set correction
     289             : 
     290           0 :             ALLOCATE (u1(atoms%jmtd, 3, mnobd, (atoms%lmaxd + 1)**2, atoms%nat), stat=ok)!hybrid%nbands
     291           0 :             IF (ok /= 0) STOP 'kp_perturbation: failure allocation u1'
     292           0 :             ALLOCATE (u2(atoms%jmtd, 3, mnobd, (atoms%lmaxd + 1)**2, atoms%nat), stat=ok)!hybrid%nbands
     293           0 :             IF (ok /= 0) STOP 'kp_perturbation: failure allocation u2'
     294           0 :             u1 = 0; u2 = 0
     295             : 
     296           0 :             iatom = 0
     297           0 :             DO itype = 1, atoms%ntype
     298           0 :                DO ieq = 1, atoms%neq(itype)
     299           0 :                   iatom = iatom + 1
     300             : 
     301           0 :                   lm1 = 0
     302           0 :                   DO l1 = 0, atoms%lmax(itype)! + 1
     303           0 :                      DO m1 = -l1, l1
     304           0 :                         lm1 = lm1 + 1
     305             : 
     306           0 :                         DO p1 = 1, 2
     307             : 
     308           0 :                            carr2 = 0
     309           0 :                            l2 = l1 + 1
     310           0 :                            lmp2 = 2*l2**2 + p1
     311           0 :                            DO m2 = -l2, l2
     312           0 :                               carr = gauntvec(l1, m1, l2, m2, atoms)
     313           0 :                               DO iband = 1, mnobd! hybrid%nbands
     314           0 :                                  carr2(1:3, iband) = carr2(1:3, iband) + carr*cmt_apw(iband, lmp2, iatom)
     315             :                               END DO
     316           0 :                               lmp2 = lmp2 + 2
     317             :                            END DO
     318             : 
     319           0 :                            DO iband = 1, mnobd! hybrid%nbands
     320           0 :                               DO i = 1, 3
     321           0 :                                  DO ig = 1, atoms%jri(itype)
     322             :                                     ! the r factor is already included in bas1
     323           0 :                                     u1(ig, i, iband, lm1, iatom) = u1(ig, i, iband, lm1, iatom) - img*bas1_tmp(ig, p1, l2, itype)*carr2(i, iband)
     324           0 :                                     u2(ig, i, iband, lm1, iatom) = u2(ig, i, iband, lm1, iatom) - img*bas2_tmp(ig, p1, l2, itype)*carr2(i, iband)
     325             :                                  END DO
     326             :                               END DO
     327             :                            END DO
     328             : 
     329           0 :                            l2 = l1 - 1
     330           0 :                            IF (l2 >= 0) THEN
     331           0 :                               carr2 = 0
     332           0 :                               lmp2 = 2*l2**2 + p1
     333           0 :                               DO m2 = -l2, l2
     334           0 :                                  carr = gauntvec(l1, m1, l2, m2, atoms)
     335           0 :                                  DO iband = 1, mnobd! hybrid%nbands
     336           0 :                                     carr2(1:3, iband) = carr2(1:3, iband) + carr*cmt_apw(iband, lmp2, iatom)
     337             :                                  END DO
     338           0 :                                  lmp2 = lmp2 + 2
     339             :                               END DO
     340             : 
     341           0 :                               DO iband = 1, mnobd! hybrid%nbands
     342           0 :                                  DO i = 1, 3
     343           0 :                                     DO ig = 1, atoms%jri(itype)
     344             :                                        ! the r factor is already included in bas1
     345           0 :                                        u1(ig, i, iband, lm1, iatom) = u1(ig, i, iband, lm1, iatom) - img*bas1_tmp(ig, p1, l2, itype)*carr2(i, iband)
     346           0 :                                        u2(ig, i, iband, lm1, iatom) = u2(ig, i, iband, lm1, iatom) - img*bas2_tmp(ig, p1, l2, itype)*carr2(i, iband)
     347             :                                     END DO
     348             :                                  END DO
     349             :                               END DO
     350             : 
     351             :                            END IF
     352             : 
     353           0 :                            carr2 = 0
     354           0 :                            l2 = l1 + 1
     355           0 :                            lmp2 = 2*l2**2
     356           0 :                            DO m2 = -l2, l2
     357           0 :                               carr = gauntvec(l1, m1, l2, m2, atoms)
     358           0 :                               DO p2 = 1, 2
     359           0 :                                  lmp2 = lmp2 + 1
     360           0 :                                  rdum = w(p1, l1, p2, l2, itype, bas1_MT_tmp, drbas1_MT_tmp, atoms%rmt)
     361           0 :                                  DO iband = 1, mnobd! hybrid%nbands
     362           0 :                                     carr2(1:3, iband) = carr2(1:3, iband) + img*carr*rdum*cmt_apw(iband, lmp2, iatom)
     363             :                                  END DO
     364             :                               END DO
     365             :                            END DO
     366             : 
     367           0 :                            DO iband = 1, mnobd! hybrid%nbands
     368           0 :                               DO i = 1, 3
     369           0 :                                  DO ig = 1, atoms%jri(itype)
     370           0 :                                     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)
     371           0 :                                     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)
     372             :                                  END DO
     373             :                               END DO
     374             :                            END DO
     375             : 
     376           0 :                            l2 = l1 - 1
     377           0 :                            IF (l2 >= 0) THEN
     378           0 :                               carr2 = 0
     379           0 :                               lmp2 = 2*l2**2
     380           0 :                               DO m2 = -l2, l2
     381           0 :                                  carr = gauntvec(l1, m1, l2, m2, atoms)
     382           0 :                                  DO p2 = 1, 2
     383           0 :                                     lmp2 = lmp2 + 1
     384           0 :                                     rdum = w(p1, l1, p2, l2, itype, bas1_MT_tmp, drbas1_MT_tmp, atoms%rmt)
     385           0 :                                     DO iband = 1, mnobd! hybrid%nbands
     386           0 :                                        carr2(1:3, iband) = carr2(1:3, iband) + img*carr*rdum*cmt_apw(iband, lmp2, iatom)
     387             :                                     END DO
     388             :                                  END DO
     389             :                               END DO
     390             : 
     391           0 :                               DO iband = 1, mnobd! hybrid%nbands
     392           0 :                                  DO i = 1, 3
     393           0 :                                     DO ig = 1, atoms%jri(itype)
     394             :                                        u1(ig, i, iband, lm1, iatom) = u1(ig, i, iband, lm1, iatom) &
     395           0 :                                                                       + bas1_tmp(ig, p1, l1, itype)*carr2(i, iband)/atoms%rmsh(ig, itype)
     396             :                                        u2(ig, i, iband, lm1, iatom) = u2(ig, i, iband, lm1, iatom) &
     397           0 :                                                                       + bas2_tmp(ig, p1, l1, itype)*carr2(i, iband)/atoms%rmsh(ig, itype)
     398             :                                     END DO
     399             :                                  END DO
     400             :                               END DO
     401             :                            END IF
     402             : 
     403             :                         END DO  ! p1
     404             :                      END DO  !m1
     405             :                   END DO  !l1
     406             :                END DO  !ieq
     407             :             END DO  !iatom
     408             : 
     409             :             ! construct lo contribtution
     410           0 :             IF (any(atoms%llo == atoms%lmaxd)) STOP 'ibs_correction: atoms%llo=atoms%lmaxd is not implemented'
     411             : 
     412           0 :             iatom = 0
     413           0 :             DO itype = 1, atoms%ntype
     414           0 :                DO ieq = 1, atoms%neq(itype)
     415           0 :                   hybrid%nindx = 2
     416           0 :                   iatom = iatom + 1
     417           0 :                   DO ilo = 1, atoms%nlo(itype)
     418           0 :                      l1 = atoms%llo(ilo, itype)
     419           0 :                      hybrid%nindx(l1, itype) = hybrid%nindx(l1, itype) + 1
     420           0 :                      p1 = hybrid%nindx(l1, itype)
     421             : 
     422           0 :                      l2 = l1 + 1
     423           0 :                      lm2 = l2**2
     424           0 :                      DO m2 = -l2, l2
     425           0 :                         lm2 = lm2 + 1
     426           0 :                         carr2 = 0
     427             : 
     428           0 :                         DO m1 = -l1, l1
     429           0 :                            carr = gauntvec(l2, m2, l1, m1, atoms)
     430           0 :                            DO iband = 1, mnobd
     431           0 :                               carr2(1:3, iband) = carr2(1:3, iband) + cmt_lo(iband, m1, ilo, iatom)*carr
     432             :                            END DO
     433             :                         END DO
     434             : 
     435           0 :                         DO iband = 1, mnobd
     436           0 :                            DO i = 1, 3
     437           0 :                               DO ig = 1, atoms%jri(itype)
     438             :                                  ! the r factor is already included in
     439           0 :                                  u1(ig, i, iband, lm2, iatom) = u1(ig, i, iband, lm2, iatom) - img*u1_lo(ig, ilo, itype)*carr2(i, iband)
     440           0 :                                  u2(ig, i, iband, lm2, iatom) = u2(ig, i, iband, lm2, iatom) - img*u2_lo(ig, ilo, itype)*carr2(i, iband)
     441             :                               END DO
     442             :                            END DO
     443             :                         END DO
     444             : 
     445             :                      END DO
     446             : 
     447           0 :                      l2 = l1 - 1
     448           0 :                      IF (l2 >= 0) THEN
     449           0 :                         lm2 = l2**2
     450           0 :                         DO m2 = -l2, l2
     451           0 :                            lm2 = lm2 + 1
     452           0 :                            carr2 = 0
     453             : 
     454           0 :                            DO m1 = -l1, l1
     455           0 :                               carr = gauntvec(l2, m2, l1, m1, atoms)
     456           0 :                               DO iband = 1, mnobd
     457           0 :                                  carr2(1:3, iband) = carr2(1:3, iband) + cmt_lo(iband, m1, ilo, iatom)*carr
     458             :                               END DO
     459             :                            END DO
     460             : 
     461           0 :                            DO iband = 1, mnobd
     462           0 :                               DO i = 1, 3
     463           0 :                                  DO ig = 1, atoms%jri(itype)
     464             :                                     ! the r factor is already included in
     465           0 :                                     u1(ig, i, iband, lm2, iatom) = u1(ig, i, iband, lm2, iatom) - img*u1_lo(ig, ilo, itype)*carr2(i, iband)
     466           0 :                                     u2(ig, i, iband, lm2, iatom) = u2(ig, i, iband, lm2, iatom) - img*u2_lo(ig, ilo, itype)*carr2(i, iband)
     467             :                                  END DO
     468             :                               END DO
     469             :                            END DO
     470             : 
     471             :                         END DO
     472             :                      END IF
     473             : 
     474             :                   END DO
     475             :                END DO
     476             :             END DO
     477             : 
     478             :             !
     479             :             ! calculate projection < phi(n',k)|phi^1(n,k)>
     480             :             ! n' = iband1 , n= iband2
     481             :             !
     482             :             iatom = 0
     483           0 :             proj_ibsc = 0
     484           0 :             DO itype = 1, atoms%ntype
     485           0 :                DO ieq = 1, atoms%neq(itype)
     486           0 :                   iatom = iatom + 1
     487           0 :                   lm = 0
     488           0 :                   lmp = 0
     489           0 :                   DO l = 0, atoms%lmax(itype)
     490           0 :                      DO M = -l, l
     491           0 :                         lm = lm + 1
     492           0 :                         ru1 = real(u1(:, :, :, lm, iatom))
     493           0 :                         iu1 = aimag(u1(:, :, :, lm, iatom))
     494           0 :                         ru2 = real(u2(:, :, :, lm, iatom))
     495           0 :                         iu2 = aimag(u2(:, :, :, lm, iatom))
     496           0 :                         DO p = 1, 2
     497           0 :                            lmp = lmp + 1
     498             : 
     499           0 :                            DO iband = 1, mnobd! hybrid%nbands
     500           0 :                               DO i = 1, 3
     501             : 
     502           0 :                                  rintegrand = atoms%rmsh(:, itype)*(hybdat%bas1(:, p, l, itype)*ru1(:, i, iband) + hybdat%bas2(:, p, l, itype)*ru2(:, i, iband))
     503             : 
     504           0 :                                  iintegrand = atoms%rmsh(:, itype)*(hybdat%bas1(:, p, l, itype)*iu1(:, i, iband) + hybdat%bas2(:, p, l, itype)*iu2(:, i, iband))
     505             : 
     506             :                                  carr2(i, iband) = intgrf(rintegrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf) &
     507           0 :                                                    + img*intgrf(iintegrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
     508             : 
     509             :                               END DO
     510             :                            END DO
     511             : 
     512           0 :                            DO iband1 = 1, hybrid%nbands(nk)
     513           0 :                               cdum = conjg(cmt_apw(iband1, lmp, iatom))
     514           0 :                               DO iband2 = 1, mnobd! hybrid%nbands
     515           0 :                                  proj_ibsc(1:3, iband2, iband1) = proj_ibsc(1:3, iband2, iband1) + cdum*carr2(1:3, iband2)
     516             :                               END DO
     517             :                            END DO
     518             : 
     519             :                         END DO!p
     520             :                      END DO!M
     521             :                   END DO!l
     522             : 
     523             :                END DO!ieq
     524             :             END DO!itype
     525             : 
     526           0 :             iatom = 0
     527           0 :             DO itype = 1, atoms%ntype
     528           0 :                DO ieq = 1, atoms%neq(itype)
     529           0 :                   iatom = iatom + 1
     530           0 :                   DO ilo = 1, atoms%nlo(itype)
     531           0 :                      l = atoms%llo(ilo, itype)
     532           0 :                      lm = l**2
     533           0 :                      DO M = -l, l
     534           0 :                         lm = lm + 1
     535           0 :                         ru1 = real(u1(:, :, :, lm, iatom))
     536           0 :                         iu1 = aimag(u1(:, :, :, lm, iatom))
     537           0 :                         ru2 = real(u2(:, :, :, lm, iatom))
     538           0 :                         iu2 = aimag(u2(:, :, :, lm, iatom))
     539             : 
     540           0 :                         DO iband = 1, mnobd! hybrid%nbands
     541           0 :                            DO i = 1, 3
     542             : 
     543           0 :                               rintegrand = atoms%rmsh(:, itype)*(u1_lo(:, ilo, itype)*ru1(:, i, iband) + u2_lo(:, ilo, itype)*ru2(:, i, iband))
     544             : 
     545           0 :                               iintegrand = atoms%rmsh(:, itype)*(u1_lo(:, ilo, itype)*iu1(:, i, iband) + u2_lo(:, ilo, itype)*iu2(:, i, iband))
     546             : 
     547             :                               carr2(i, iband) = intgrf(rintegrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf) &
     548           0 :                                                 + img*intgrf(iintegrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
     549             : 
     550             :                            END DO
     551             :                         END DO
     552             : 
     553           0 :                         DO iband1 = 1, hybrid%nbands(nk)
     554           0 :                            cdum = conjg(cmt_lo(iband1, M, ilo, iatom))
     555           0 :                            DO iband2 = 1, mnobd! hybrid%nbands
     556           0 :                               proj_ibsc(1:3, iband2, iband1) = proj_ibsc(1:3, iband2, iband1) + cdum*carr2(1:3, iband2)
     557             :                            END DO
     558             :                         END DO
     559             : 
     560             :                      END DO
     561             : 
     562             :                   END DO
     563             :                END DO
     564             :             END DO
     565             : 
     566             :             !
     567             :             ! calculate <phi^1(n1,k)|phi^1(n2,k)>
     568             :             ! n1 and n2 occupied
     569             :             !
     570             :             iatom = 0
     571           0 :             olap_ibsc = 0
     572           0 :             DO itype = 1, atoms%ntype
     573           0 :                DO ieq = 1, atoms%neq(itype)
     574           0 :                   iatom = iatom + 1
     575           0 :                   lm = 0
     576           0 :                   DO l = 0, atoms%lmax(itype)!+1
     577           0 :                      DO M = -l, l
     578           0 :                         lm = lm + 1
     579           0 :                         ru1 = real(u1(:, :, :, lm, iatom))
     580           0 :                         iu1 = aimag(u1(:, :, :, lm, iatom))
     581           0 :                         ru2 = real(u2(:, :, :, lm, iatom))
     582           0 :                         iu2 = aimag(u2(:, :, :, lm, iatom))
     583             : 
     584           0 :                         DO iband1 = 1, mnobd ! hybrid%nbands
     585           0 :                            DO iband2 = 1, mnobd!iband1
     586           0 :                               DO i = 1, 3
     587           0 :                                  DO j = 1, 3
     588             : 
     589             :                                     rintegrand = atoms%rmsh(:, itype)**2*(ru1(:, i, iband1)*ru1(:, j, iband2) + ru2(:, i, iband1)*ru2(:, j, iband2) &
     590           0 :                                                                           + iu1(:, i, iband1)*iu1(:, j, iband2) + iu2(:, i, iband1)*iu2(:, j, iband2))
     591             : 
     592             :                                     iintegrand = atoms%rmsh(:, itype)**2*(ru1(:, i, iband1)*iu1(:, j, iband2) + ru2(:, i, iband1)*iu2(:, j, iband2) &
     593           0 :                                                                           - iu1(:, i, iband1)*ru1(:, j, iband2) - iu2(:, i, iband1)*ru2(:, j, iband2))
     594             : 
     595             :                                     olap_ibsc(i, j, iband2, iband1) = olap_ibsc(i, j, iband2, iband1) &
     596             :                                                                       + intgrf(rintegrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf) &
     597           0 :                                                                       + img*intgrf(iintegrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
     598             : 
     599             :                                  END DO
     600             :                               END DO
     601             : 
     602             :                            END DO
     603             :                         END DO
     604             : 
     605             :                      END DO
     606             :                   END DO
     607             :                END DO
     608             :             END DO
     609             : 
     610           0 :          END SUBROUTINE ibs_correction
     611             : 
     612           0 :          FUNCTION gauntvec(l1, m1, l2, m2, atoms)
     613             : 
     614             :             USE m_constants
     615             :             USE m_gaunt
     616             : 
     617             :             USE m_types
     618             :             IMPLICIT NONE
     619             :             TYPE(t_atoms), INTENT(IN)   :: atoms
     620             : 
     621             :             INTEGER, INTENT(IN)  ::  l1, m1, l2, m2
     622             : 
     623             :             COMPLEX               ::  gauntvec(-1:1)
     624             : 
     625             :             INTEGER               ::  j, mj
     626             :             INTEGER               ::  point(-1:1)
     627             :             REAL                  ::  rfac
     628             :             COMPLEX               ::  cfac
     629             :             COMPLEX, PARAMETER   ::  img = (0.0, 1.0)
     630             : 
     631           0 :             rfac = sqrt(tpi_const/3)
     632           0 :             DO j = -1, 1
     633             :                ! j = -1 corresponds to x-direction
     634             :                ! j = 0  corresponds to z-direction
     635             :                ! j = 1  corresponds to y-direction
     636           0 :                mj = abs(j)
     637           0 :                cfac = img**((j + mj)/2.)*sqrt(2.)**(1 - mj)*rfac
     638           0 :                gauntvec(j) = cfac*(gaunt1(1, l1, l2, -mj, m1, m2, atoms%lmaxd + 1) + j*gaunt1(1, l1, l2, mj, m1, m2, atoms%lmaxd + 1))
     639             :             END DO
     640             : 
     641             :             ! transform onto cartesian coordinates
     642           0 :             point(-1) = -1
     643           0 :             point(0) = 1
     644           0 :             point(1) = 0
     645             : 
     646           0 :             gauntvec = gauntvec(point)
     647             : 
     648             :          END FUNCTION gauntvec
     649             : 
     650           0 :          FUNCTION w(p1, l1, p2, l2, itype, bas1_mt, drbas1_mt, &
     651           0 :                     rmt)
     652             :             USE m_types
     653             :             IMPLICIT NONE
     654             : 
     655             :             INTEGER, INTENT(IN)       ::  p1, l1, p2, l2
     656             :             INTEGER, INTENT(IN)       ::  itype
     657             :             REAL, INTENT(IN)            :: rmt(:), bas1_mt(:, 0:, :), drbas1_mt(:, 0:, :)
     658             : 
     659             :             REAL                  ::  w
     660             : 
     661             :             INTEGER               ::  p
     662             :             REAL                  ::  denom, enum
     663             : 
     664           0 :             IF (p1 > 2 .or. p2 > 2) STOP 'w: the formalism is only valid for p<=2'
     665             : 
     666           0 :             denom = wronskian(bas1_MT(2, l1, itype), drbas1_MT(2, l1, itype), bas1_MT(1, l1, itype), drbas1_MT(1, l1, itype))
     667             : 
     668           0 :             p = p1 + (-1)**(p1 - 1)
     669             : 
     670             :             enum = bas1_MT(p, l1, itype)*bas1_MT(p2, l2, itype) + rmt(itype)*wronskian(bas1_MT(p, l1, itype), &
     671           0 :                                                                                        drbas1_MT(p, l1, itype), bas1_MT(p2, l2, itype), drbas1_MT(p2, l2, itype))
     672             : 
     673           0 :             w = (-1)**(p1 + 1)*enum/denom
     674             : 
     675           0 :          END FUNCTION
     676             : 
     677           0 :          PURE FUNCTION wronskian(f, df, g, dg)
     678             :             IMPLICIT NONE
     679             :             REAL, INTENT(IN) ::  f, df, g, dg
     680             :             REAL              ::  wronskian
     681             : 
     682           0 :             wronskian = f*dg - df*g
     683             : 
     684           0 :          END FUNCTION
     685             : 
     686             : !     Calculates the derivative
     687             : !                                  ikr    s       s
     688             : !     dcprod(n',n,k,xyz) = d    < e    phi   | phi       > / sqrt(vol)
     689             : !                           xyz           qn      q+k,n'
     690             : !
     691             : !                                s             s
     692             : !                           < phi  | d    | phi    >
     693             : !                                nq   xyz      n'q
     694             : !                      = -i ------------------------ / sqrt(vol)  ,   s = ispin
     695             : !                                  s     s                        ,   n = occ.     ,   n' = unocc.
     696             : !                                 e    - e                        ,   bandi1 <= n <= bandf1 , bandi2 <= n' <= bandf2
     697             : !                                  n'q    nq
     698             : !
     699             : !     with kp perturbation theory and
     700             : !
     701             : !               d     d     d
     702             : !     d    =   ---,  ---,  --- .
     703             : !      xyz     dk    dk    dk
     704             : !                x     y     z
     705             : !
     706           0 :          SUBROUTINE dwavefproducts( &
     707           0 :             dcprod, nk, bandi1, bandf1, bandi2, bandf2, lwrite, &
     708             :             atoms, hybrid, &
     709             :             cell, &
     710             :             hybdat, kpts, nkpti, lapw, &
     711             :             dimension, jsp, &
     712           0 :             eig_irr)
     713             : 
     714             :             USE m_wrapper
     715             :             USE m_types
     716             :             IMPLICIT NONE
     717             : 
     718             :             TYPE(t_hybdat), INTENT(IN)   :: hybdat
     719             : 
     720             :             TYPE(t_dimension), INTENT(IN)   :: dimension
     721             :             TYPE(t_hybrid), INTENT(IN)   :: hybrid
     722             :             TYPE(t_cell), INTENT(IN)   :: cell
     723             :             TYPE(t_kpts), INTENT(IN)   :: kpts
     724             :             TYPE(t_atoms), INTENT(IN)   :: atoms
     725             :             TYPE(t_lapw), INTENT(IN)   :: lapw
     726             : 
     727             : !     - scalars -
     728             :             INTEGER, INTENT(IN)      ::  nk, bandi1, bandf1, bandi2, bandf2
     729             :             INTEGER, INTENT(IN)      :: nkpti
     730             :             INTEGER, INTENT(IN)      :: jsp
     731             : 
     732             : !     - arrays -
     733             : 
     734             :             REAL, INTENT(IN)         ::  eig_irr(dimension%neigd, nkpti)
     735             :             COMPLEX, INTENT(OUT)     ::  dcprod(bandi2:bandf2, bandi1:bandf1, 3)
     736             : 
     737             : !     - local scalars -
     738             :             INTEGER                 ::  ikpt, ikpt1, iband1, iband2
     739             :             REAL                    ::  rdum
     740             :             LOGICAL                 ::  lwrite
     741             : 
     742             :             !                                       __
     743             :             ! Get momentum-matrix elements -i < uj | \/ | ui >
     744             :             !
     745             :             CALL momentum_matrix( &
     746             :                dcprod, nk, bandi1, bandf1, bandi2, bandf2, &
     747             :                atoms, hybrid, &
     748             :                cell, &
     749             :                hybdat, kpts, lapw, &
     750           0 :                dimension, jsp)
     751             : 
     752             :             !                                                __
     753             :             !  Calculate expansion coefficients -i < uj | \/ | ui > / ( ei - ej ) for periodic function ui
     754             :             !
     755           0 :             DO iband1 = bandi1, bandf1
     756           0 :                DO iband2 = bandi2, bandf2
     757           0 :                   rdum = eig_irr(iband2, nk) - eig_irr(iband1, nk)
     758           0 :                   IF (abs(rdum) > 1e-6) THEN !10.0**-6
     759           0 :                      dcprod(iband2, iband1, :) = dcprod(iband2, iband1, :)/rdum
     760             :                   ELSE
     761           0 :                      dcprod(iband2, iband1, :) = 0.0
     762             :                   END IF
     763             :                END DO
     764             :             END DO
     765             : 
     766           0 :             dcprod = dcprod/sqrt(cell%omtil)
     767             : 
     768           0 :          END SUBROUTINE dwavefproducts
     769             : 
     770             : !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     771             : 
     772             : !     Calculates the momentum matrix elements
     773             : !
     774             : !                                    s             s
     775             : !     momentum(n',n,q,xyz) = -i < phi  | d    | phi    >   ,   s  = ispin
     776             : !                                    nq   xyz      n'q     ,   n  = occ.    ( bandi1 <= n  <= bandf1 )  ,
     777             : !                                                              n' = unocc.  ( bandi2 <= n' <= bandf2 )
     778             : !
     779             : !
     780           0 :          SUBROUTINE momentum_matrix( &
     781           0 :             momentum, nk, bandi1, bandf1, bandi2, bandf2, &
     782             :             atoms, hybrid, &
     783             :             cell, &
     784             :             hybdat, kpts, lapw, &
     785             :             dimension, jsp)
     786             : 
     787             :             USE m_olap
     788             :             USE m_wrapper
     789             :             USE m_util, only: derivative, intgrf_init, intgrf
     790             :             USE m_dr2fdr
     791             :             USE m_constants
     792             :             USE m_types
     793             :             USE m_io_hybrid
     794             :             IMPLICIT NONE
     795             : 
     796             :             TYPE(t_hybdat), INTENT(IN)   :: hybdat
     797             :             TYPE(t_dimension), INTENT(IN)   :: dimension
     798             :             TYPE(t_hybrid), INTENT(IN)   :: hybrid
     799             :             TYPE(t_cell), INTENT(IN)   :: cell
     800             :             TYPE(t_kpts), INTENT(IN)   :: kpts
     801             :             TYPE(t_atoms), INTENT(IN)   :: atoms
     802             :             TYPE(t_lapw), INTENT(IN)   :: lapw
     803             : 
     804             : !     - scalars -
     805             :             INTEGER, INTENT(IN)      ::  bandi1, bandf1, bandi2, bandf2
     806             :             INTEGER, INTENT(IN)      :: nk
     807             :             INTEGER, INTENT(IN)      :: jsp
     808             : 
     809             : !     - arrays -
     810           0 :             TYPE(t_mat):: z
     811             :             COMPLEX, INTENT(OUT)     :: momentum(bandi2:bandf2, bandi1:bandf1, 3)
     812             : 
     813             : !     - local scalars -
     814             :             INTEGER                 ::  itype, ieq, ic, i, j, l, lm, n1, n2, ikpt, iband1, iband2, ll, mm
     815             :             INTEGER                 ::  lm_0, lm_1, lm0, lm1, lm2, lm3, n0, nn, n, l1, l2, m1, m2, ikpt1
     816             :             INTEGER                 ::  irecl_cmt, irecl_z, m
     817             :             COMPLEX                 ::  cdum
     818             :             COMPLEX                 ::  img = (0.0, 1.0)
     819             : 
     820             : !     - local arrays -
     821           0 :             INTEGER                 ::  gpt(3, lapw%nv(jsp))
     822             : 
     823           0 :             REAL                    ::  fcoeff((atoms%lmaxd + 1)**2, -1:1), gcoeff((atoms%lmaxd + 1)**2, -1:1)
     824           0 :             REAL                    ::  qmat1(hybrid%maxindx, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype), dbas1(atoms%jmtd)
     825           0 :             REAL                    ::  qmat2(hybrid%maxindx, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype), dbas2(atoms%jmtd)
     826           0 :             REAL                    ::  qg(lapw%nv(jsp), 3)
     827             : 
     828             :             COMPLEX                 ::  hlp(3, 3)
     829           0 :             COMPLEX                 ::  cvec1(hybrid%maxlmindx), cvec2(hybrid%maxlmindx), cvec3(hybrid%maxlmindx)
     830           0 :             COMPLEX                 ::  cmt1(hybrid%maxlmindx, bandi1:bandf1), cmt2(hybrid%maxlmindx, bandi2:bandf2)
     831             :             COMPLEX                 ::  carr1(3), carr2(3)
     832           0 :             COMPLEX                 ::  cmt(dimension%neigd, hybrid%maxlmindx, atoms%nat)
     833           0 :             REAL                    ::  olap_r(lapw%nv(jsp)*(lapw%nv(jsp) + 1)/2)
     834           0 :             COMPLEX                 ::  olap_c(lapw%nv(jsp)*(lapw%nv(jsp) + 1)/2)
     835           0 :             REAL                    ::  vec1_r(lapw%nv(jsp)), vec2_r(lapw%nv(jsp)), vec3_r(lapw%nv(jsp))
     836           0 :             COMPLEX                 ::  vec1_c(lapw%nv(jsp)), vec2_c(lapw%nv(jsp)), vec3_c(lapw%nv(jsp))
     837             : 
     838             :             ! read in cmt coefficients from direct access file cmt at kpoint nk
     839             : 
     840           0 :             call read_cmt(cmt, nk)
     841             : 
     842             :             ! read in z coefficients from direct access file z at kpoint nk
     843             : 
     844           0 :             call read_z(z, nk)
     845             : 
     846             :             !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
     847             : 
     848           0 :             gpt(1, 1:lapw%nv(jsp)) = lapw%k1(1:lapw%nv(jsp), jsp)
     849           0 :             gpt(2, 1:lapw%nv(jsp)) = lapw%k2(1:lapw%nv(jsp), jsp)
     850           0 :             gpt(3, 1:lapw%nv(jsp)) = lapw%k3(1:lapw%nv(jsp), jsp)
     851             : 
     852             : !     Define coefficients F and G
     853             :             lm = 0
     854           0 :             DO l = 0, atoms%lmaxd
     855           0 :                DO M = -l, l
     856           0 :                   lm = lm + 1
     857           0 :                   fcoeff(lm, -1) = -sqrt(1.0*(l + M + 1)*(l + M + 2)/(2*(2*l + 1)*(2*l + 3)))
     858           0 :                   fcoeff(lm, 0) = sqrt(1.0*(l - M + 1)*(l + M + 1)/((2*l + 1)*(2*l + 3)))
     859           0 :                   fcoeff(lm, 1) = -sqrt(1.0*(l - M + 1)*(l - M + 2)/(2*(2*l + 1)*(2*l + 3)))
     860           0 :                   gcoeff(lm, -1) = sqrt(1.0*(l - M)*(l - M - 1)/(2*(2*l - 1)*(2*l + 1)))
     861           0 :                   gcoeff(lm, 0) = sqrt(1.0*(l - M)*(l + M)/((2*l - 1)*(2*l + 1)))
     862           0 :                   gcoeff(lm, 1) = sqrt(1.0*(l + M)*(l + M - 1)/(2*(2*l - 1)*(2*l + 1)))
     863             :                END DO
     864             :             END DO
     865             : 
     866             : !     Calculate olap int r**2*u*u' + w * int r*u*u, w = -l,l+1 ( -> qmat1/2 )
     867           0 :             qmat1 = 0
     868           0 :             qmat2 = 0
     869           0 :             ic = 0
     870           0 :             DO itype = 1, atoms%ntype
     871           0 :                DO l = 0, atoms%lmax(itype)
     872           0 :                   DO n2 = 1, hybrid%nindx(l, itype)
     873             :                      !ic = ic + 1
     874           0 :                      CALL derivative(dbas1, hybdat%bas1(:, n2, l, itype), atoms, itype)
     875           0 :                      dbas1 = dbas1 - hybdat%bas1(:, n2, l, itype)/atoms%rmsh(:, itype)
     876             : 
     877           0 :                      CALL derivative(dbas2, hybdat%bas2(:, n2, l, itype), atoms, itype)
     878           0 :                      dbas2 = dbas2 - hybdat%bas2(:, n2, l, itype)/atoms%rmsh(:, itype)
     879             : 
     880           0 :                      IF (l /= 0) THEN
     881           0 :                         DO n1 = 1, hybrid%nindx(l - 1, itype)
     882           0 :                            ic = ic + 1
     883             :                            qmat1(n1, n2, l, itype) = intgrf(dbas1(:)*hybdat%bas1(:, n1, l - 1, itype) + &
     884             :                                                             dbas2(:)*hybdat%bas2(:, n1, l - 1, itype), atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf) &
     885             :                                                      + intgrf((hybdat%bas1(:, n2, l, itype)*hybdat%bas1(:, n1, l - 1, itype) + hybdat%bas2(:, n2, l, itype)*hybdat%bas2(:, n1, l - 1, itype)) &
     886           0 :                                                               /atoms%rmsh(:, itype), atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)*(l + 1)
     887             : 
     888             :                         END DO
     889             :                      END IF
     890           0 :                      IF (l /= atoms%lmax(itype)) THEN
     891           0 :                         DO n1 = 1, hybrid%nindx(l + 1, itype)
     892             : 
     893             :                            qmat2(n1, n2, l, itype) = intgrf(dbas1(:)*hybdat%bas1(:, n1, l + 1, itype) + dbas2(:)*hybdat%bas2(:, n1, l + 1, itype), &
     894             :                                                             atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf) &
     895             :                                                      - intgrf((hybdat%bas1(:, n2, l, itype)*hybdat%bas1(:, n1, l + 1, itype) + hybdat%bas2(:, n2, l, itype)*hybdat%bas2(:, n1, l + 1, itype)) &
     896           0 :                                                               /atoms%rmsh(:, itype), atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)*l
     897             : 
     898             :                         END DO
     899             :                      END IF
     900             : 
     901             :                   END DO
     902             :                END DO
     903             :             END DO
     904             : 
     905             :             !                                                  __
     906             :             ! Calculate momentum matrix elements -i < uj | \/ | ui > wrt wave functions u (->momentum)
     907             :             !
     908             : 
     909           0 :             momentum = 0
     910             : 
     911             :             ! MT contribution
     912             : 
     913           0 :             ic = 0
     914           0 :             DO itype = 1, atoms%ntype
     915           0 :                DO ieq = 1, atoms%neq(itype)
     916           0 :                   ic = ic + 1
     917           0 :                   nn = sum((/((2*l + 1)*hybrid%nindx(l, itype), l=0, atoms%lmax(itype))/))
     918           0 :                   DO iband1 = bandi1, bandf1
     919           0 :                      cmt1(:nn, iband1) = cmt(iband1, :nn, ic)
     920             :                   ENDDO
     921           0 :                   DO iband2 = bandi2, bandf2
     922           0 :                      cmt2(:nn, iband2) = cmt(iband2, :nn, ic)
     923             :                   ENDDO
     924           0 :                   DO iband1 = bandi1, bandf1
     925             : 
     926           0 :                      cvec1 = 0; cvec2 = 0; cvec3 = 0
     927             :                      ! build up left vector(s) ( -> cvec1/2/3 )
     928           0 :                      lm_0 = 0              ! we start with s-functions (l=0)
     929           0 :                      lm_1 = hybrid%nindx(0, itype) ! we start with p-functions (l=0+1)
     930           0 :                      lm = 0
     931           0 :                      DO l = 0, atoms%lmax(itype) - 1
     932           0 :                         n0 = hybrid%nindx(l, itype)
     933           0 :                         n1 = hybrid%nindx(l + 1, itype)
     934           0 :                         DO M = -l, l
     935           0 :                            lm = lm + 1
     936           0 :                            lm0 = lm_0 + (M + l)*n0
     937           0 :                            lm1 = lm_1 + (M + 1 + l + 1)*n1
     938           0 :                            lm2 = lm_1 + (M + l + 1)*n1
     939           0 :                            lm3 = lm_1 + (M - 1 + l + 1)*n1
     940           0 :                            cvec1(lm0 + 1:lm0 + n0) = fcoeff(lm, -1)*matmul(cmt1(lm1 + 1:lm1 + n1, iband1), qmat2(:n1, :n0, l, itype))
     941           0 :                            cvec2(lm0 + 1:lm0 + n0) = fcoeff(lm, 0)*matmul(cmt1(lm2 + 1:lm2 + n1, iband1), qmat2(:n1, :n0, l, itype))
     942           0 :                            cvec3(lm0 + 1:lm0 + n0) = fcoeff(lm, 1)*matmul(cmt1(lm3 + 1:lm3 + n1, iband1), qmat2(:n1, :n0, l, itype))
     943             :                         END DO
     944           0 :                         lm_0 = lm_0 + (2*l + 1)*n0
     945           0 :                         lm_1 = lm_1 + (2*l + 3)*n1
     946             :                      END DO
     947             : 
     948           0 :                      lm_0 = hybrid%nindx(0, itype) ! we start with p-functions (l=1)
     949           0 :                      lm_1 = 0              ! we start with s-functions (l=1-1)
     950           0 :                      lm = 1
     951           0 :                      DO l = 1, atoms%lmax(itype)
     952           0 :                         n0 = hybrid%nindx(l, itype)
     953           0 :                         n1 = hybrid%nindx(l - 1, itype)
     954           0 :                         DO M = -l, l
     955           0 :                            lm = lm + 1
     956           0 :                            lm0 = lm_0 + (M + l)*n0
     957           0 :                            lm1 = lm_1 + (M + 1 + l - 1)*n1
     958           0 :                            lm2 = lm_1 + (M + l - 1)*n1
     959           0 :                            lm3 = lm_1 + (M - 1 + l - 1)*n1
     960           0 :                            IF (abs(M + 1) <= 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))
     961           0 :                            IF (abs(M) <= 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))
     962           0 :                            IF (abs(M - 1) <= 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))
     963             :                         END DO
     964           0 :                         lm_0 = lm_0 + (2*l + 1)*n0
     965           0 :                         lm_1 = lm_1 + (2*l - 1)*n1
     966             :                      END DO
     967             :                      ! multiply with right vector
     968           0 :                      DO iband2 = bandi2, bandf2
     969           0 :                         momentum(iband2, iband1, 1) = momentum(iband2, iband1, 1) + dotprod(cvec1(:nn), cmt2(:nn, iband2))
     970           0 :                         momentum(iband2, iband1, 2) = momentum(iband2, iband1, 2) + dotprod(cvec2(:nn), cmt2(:nn, iband2))
     971           0 :                         momentum(iband2, iband1, 3) = momentum(iband2, iband1, 3) + dotprod(cvec3(:nn), cmt2(:nn, iband2))
     972             :                      END DO ! iband2
     973             :                   END DO ! iband1
     974             : 
     975             :                END DO ! ieq
     976             :             END DO ! itype
     977             : 
     978             :             ! Transform to cartesian coordinates
     979           0 :             hlp = 0
     980           0 :             hlp(1, 1) = 1.0/sqrt(2.0)
     981           0 :             hlp(1, 3) = -1.0/sqrt(2.0)
     982           0 :             hlp(2, 1) = -img/sqrt(2.0)
     983           0 :             hlp(2, 3) = -img/sqrt(2.0)
     984           0 :             hlp(3, 2) = 1.0
     985           0 :             DO iband1 = bandi1, bandf1
     986           0 :                DO iband2 = bandi2, bandf2
     987           0 :                   momentum(iband2, iband1, :) = -img*matmul(momentum(iband2, iband1, :), transpose(hlp))
     988             :                END DO
     989             :             END DO
     990             : 
     991             :             ! plane-wave contribution
     992             : 
     993           0 :             CALL olap_pwp(z%l_real, olap_r, olap_c, gpt, lapw%nv(jsp), atoms, cell)
     994             : 
     995           0 :             DO nn = 1, lapw%nv(jsp)
     996           0 :                qg(nn, :) = matmul(kpts%bk(:, nk) + gpt(:, nn), cell%bmat)
     997             :             END DO
     998             : 
     999           0 :             if (z%l_real) THEN
    1000           0 :             DO iband2 = bandi2, bandf2
    1001           0 :                vec1_r = matvec(olap_r, z%data_r(:lapw%nv(jsp), iband2)*qg(:, 1))
    1002           0 :                vec2_r = matvec(olap_r, z%data_r(:lapw%nv(jsp), iband2)*qg(:, 2))
    1003           0 :                vec3_r = matvec(olap_r, z%data_r(:lapw%nv(jsp), iband2)*qg(:, 3))
    1004           0 :                DO iband1 = bandi1, bandf1
    1005           0 :                   momentum(iband2, iband1, 1) = momentum(iband2, iband1, 1) + dotprod(z%data_r(:lapw%nv(jsp), iband1), vec1_r)
    1006           0 :                   momentum(iband2, iband1, 2) = momentum(iband2, iband1, 2) + dotprod(z%data_r(:lapw%nv(jsp), iband1), vec2_r)
    1007           0 :                   momentum(iband2, iband1, 3) = momentum(iband2, iband1, 3) + dotprod(z%data_r(:lapw%nv(jsp), iband1), vec3_r)
    1008             :                END DO
    1009             :             END DO
    1010             :             else
    1011           0 :             DO iband2 = bandi2, bandf2
    1012           0 :                vec1_c = matvec(olap_c, z%data_c(:lapw%nv(jsp), iband2)*qg(:, 1))
    1013           0 :                vec2_c = matvec(olap_c, z%data_c(:lapw%nv(jsp), iband2)*qg(:, 2))
    1014           0 :                vec3_c = matvec(olap_c, z%data_c(:lapw%nv(jsp), iband2)*qg(:, 3))
    1015           0 :                DO iband1 = bandi1, bandf1
    1016           0 :                   momentum(iband2, iband1, 1) = momentum(iband2, iband1, 1) + dotprod(z%data_c(:lapw%nv(jsp), iband1), vec1_c)
    1017           0 :                   momentum(iband2, iband1, 2) = momentum(iband2, iband1, 2) + dotprod(z%data_c(:lapw%nv(jsp), iband1), vec2_c)
    1018           0 :                   momentum(iband2, iband1, 3) = momentum(iband2, iband1, 3) + dotprod(z%data_c(:lapw%nv(jsp), iband1), vec3_c)
    1019             :                END DO
    1020             :             END DO
    1021             :             end if
    1022             : 
    1023           0 :          END SUBROUTINE momentum_matrix
    1024             : 
    1025             :       END MODULE m_kp_perturbation

Generated by: LCOV version 1.13