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

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : MODULE m_wavefproducts
       7             :    USE m_judft
       8             :    PRIVATE
       9             :    PUBLIC wavefproducts_noinv, wavefproducts_noinv5
      10             :    PUBLIC wavefproducts_inv, wavefproducts_inv5
      11             : CONTAINS
      12             : 
      13           0 :    SUBROUTINE wavefproducts_noinv(bandi, bandf, nk, iq, dimension, input, jsp,&                  !cprod,&
      14             :   &                 cell, atoms, hybrid, hybdat,&
      15             :   &                 kpts, mnobd,&
      16           0 :   &                 lapw, sym, noco, nbasm_mt, nkqpt, cprod)
      17             : 
      18             :       USE m_constants
      19             :       USE m_util, ONLY: modulo1
      20             :       USE m_wrapper
      21             :       USE m_types
      22             :       USE m_io_hybrid
      23             :       IMPLICIT NONE
      24             : 
      25             :       TYPE(t_input), INTENT(IN)       :: input
      26             :       TYPE(t_dimension), INTENT(IN)   :: dimension
      27             :       TYPE(t_hybrid), INTENT(IN)      :: hybrid
      28             :       TYPE(t_sym), INTENT(IN)         :: sym
      29             :       TYPE(t_noco), INTENT(IN)        :: noco
      30             :       TYPE(t_cell), INTENT(IN)        :: cell
      31             :       TYPE(t_kpts), INTENT(IN)        :: kpts
      32             :       TYPE(t_atoms), INTENT(IN)       :: atoms
      33             :       TYPE(t_lapw), INTENT(IN)        :: lapw
      34             :       TYPE(t_hybdat), INTENT(INOUT)   :: hybdat
      35             : 
      36             : !     - scalars -
      37             :       INTEGER, INTENT(IN)      ::  nk, iq, jsp
      38             :       INTEGER, INTENT(IN)      :: mnobd
      39             :       INTEGER, INTENT(IN)      :: nbasm_mt
      40             :       INTEGER, INTENT(IN)      ::  bandi, bandf
      41             :       INTEGER, INTENT(OUT)     ::  nkqpt
      42             : 
      43             : !     - arrays -
      44             : 
      45             :       COMPLEX, INTENT(OUT)    ::  cprod(hybrid%maxbasm1, mnobd, bandf - bandi + 1)
      46             : 
      47             : !     - local scalars -
      48             :       INTEGER                 ::  ic, l, n, l1, l2, n1, n2, lm_0, lm1_0, lm2_0, lm, lm1, lm2, m1, m2, i, j, ll
      49             :       INTEGER                 ::  itype, ieq, ikpt, ikpt1, ikpt2, igpt, igptp, igpt1, igpt2, iband, iband1, iband2
      50             :       INTEGER                 ::  k, ic1, ioffset, ibando
      51             :       INTEGER                 ::  q, idum, m
      52             :       INTEGER                 ::  nbasm_ir
      53             :       INTEGER                 ::  nbasmmt, nbasfcn
      54             :       INTEGER                 ::  ok
      55             :       REAL                    ::  rdum, svol, s2, pi
      56             :       REAL                    ::  mtthr = 0
      57             :       COMPLEX                 ::  cdum, cdum0
      58             :       COMPLEX, PARAMETER       ::  img = (0.0, 1.0)
      59             : 
      60             :       LOGICAL                 ::  offdiag
      61             : !      - local arrays -
      62           0 :       INTEGER                 ::  iarr(lapw%nv(jsp))
      63           0 :       INTEGER                 ::  g(3), ghelp(3), lmstart(0:atoms%lmaxd, atoms%ntype)
      64             :       INTEGER                 ::  gpthlp(3, dimension%nvd), nvhlp(input%jspins)
      65           0 :       INTEGER                 ::  gpt_nk(3, lapw%nv(jsp))
      66             :       INTEGER                 ::  g_t(3)
      67             :       INTEGER                 ::  gsum(3)
      68             :       REAL                    ::  bkpt(3)
      69             :       REAL                    ::  bkhlp(3)
      70             :       REAL                    ::  kqpt(3), kqpthlp(3)
      71           0 :       COMPLEX                 ::  carr(1:mnobd, bandf - bandi + 1)
      72             : !      COMPLEX                 :: chelp(maxbasm,mnobd,bandf-bandi+1,nkpt_EIBZ)
      73             :       COMPLEX                 ::  cexp
      74           0 :       COMPLEX                 ::  z_help(lapw%nv(jsp))
      75           0 :       COMPLEX                 ::  cmt(dimension%neigd, hybrid%maxlmindx, atoms%nat)
      76           0 :       COMPLEX                 ::  cmt_nk(dimension%neigd, hybrid%maxlmindx, atoms%nat)
      77           0 :       COMPLEX, ALLOCATABLE     ::  cprod_ir(:, :, :)
      78           0 :       TYPE(t_mat)             :: z_nk, z_kqpt
      79           0 :       TYPE(t_lapw)            :: lapw_nkqpt
      80           0 :       CALL timestart("wavefproducts_noinv")
      81             : 
      82             :       ! preparations
      83             : 
      84           0 :       svol = sqrt(cell%omtil)
      85           0 :       s2 = sqrt(2.0)
      86             : 
      87           0 :       gpt_nk(1, :) = lapw%k1(:lapw%nv(jsp), jsp)
      88           0 :       gpt_nk(2, :) = lapw%k2(:lapw%nv(jsp), jsp)
      89           0 :       gpt_nk(3, :) = lapw%k3(:lapw%nv(jsp), jsp)
      90             : 
      91             :       !
      92             :       ! compute k+q point for given q point in EIBZ(k)
      93             :       !
      94           0 :       kqpthlp = kpts%bkf(:, nk) + kpts%bkf(:, iq)
      95             :       ! k+q can lie outside the first BZ, transfer
      96             :       ! it back into the 1. BZ
      97           0 :       kqpt = modulo1(kqpthlp, kpts%nkpt3)
      98           0 :       g_t(:) = nint(kqpt - kqpthlp)
      99             :       ! determine number of kqpt
     100           0 :       nkqpt = 0
     101           0 :       DO ikpt = 1, kpts%nkptf
     102           0 :          IF (maxval(abs(kqpt - kpts%bkf(:, ikpt))) <= 1E-06) THEN
     103           0 :             nkqpt = ikpt
     104           0 :             EXIT
     105             :          END IF
     106             :       END DO
     107           0 :       IF (nkqpt == 0) STOP 'wavefproducts: k-point not found'
     108             : 
     109             :       ! lmstart = lm start index for each l-quantum number and atom type (for cmt-coefficients)
     110           0 :       DO itype = 1, atoms%ntype
     111           0 :          DO l = 0, atoms%lmax(itype)
     112           0 :             lmstart(l, itype) = sum((/(hybrid%nindx(ll, itype)*(2*ll + 1), ll=0, l - 1)/))
     113             :          END DO
     114             :       END DO
     115             : 
     116           0 :       nbasm_ir = maxval(hybrid%ngptm)
     117           0 :       ALLOCATE (cprod_ir(bandf - bandi + 1, mnobd, nbasm_ir), stat=ok)
     118           0 :       IF (ok /= 0) STOP 'wavefproducts: failure allocation cprod_ir'
     119           0 :       cprod_ir = 0
     120             : 
     121           0 :       cprod = 0
     122             : 
     123           0 :       CALL lapw_nkqpt%init(input, noco, kpts, atoms, sym, nkqpt, cell, sym%zrfs)
     124           0 :       nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
     125           0 :       call z_nk%alloc(.false., nbasfcn, dimension%neigd)
     126           0 :       nbasfcn = MERGE(lapw_nkqpt%nv(1) + lapw_nkqpt%nv(2) + 2*atoms%nlotot, lapw_nkqpt%nv(1) + atoms%nlotot, noco%l_noco)
     127           0 :       call z_kqpt%alloc(.false., nbasfcn, dimension%neigd)
     128             : 
     129           0 :       call read_z(z_nk, nk)
     130           0 :       call read_z(z_kqpt, nkqpt)
     131             : 
     132             :       ! read in cmt coefficients from direct access file cmt
     133           0 :       call read_cmt(cmt_nk, nk)
     134             : 
     135             :       ! IR contribution
     136             : 
     137           0 :       CALL timestart("wavefproducts_noinv IR")
     138             : 
     139           0 :       DO igpt = 1, hybrid%ngptm(iq)
     140           0 :          igptp = hybrid%pgptm(igpt, iq)
     141           0 :          ghelp = hybrid%gptm(:, igptp) - g_t(:)
     142           0 :          DO i = 1, lapw%nv(jsp)
     143           0 :             gsum(:) = ghelp + gpt_nk(:, i)
     144           0 :             IF (all(abs(gsum) <= hybdat%pntgptd)) THEN
     145           0 :                iarr(i) = hybdat%pntgpt(gsum(1), gsum(2), gsum(3), nkqpt)
     146             :             ELSE
     147           0 :                iarr(i) = 0
     148             :             END IF
     149             : 
     150             :          END DO
     151             : 
     152           0 :          DO iband1 = 1, hybrid%nobd(nkqpt)
     153           0 :             where (iarr > 0)
     154             :             z_help(:) = z_kqpt%data_c(iarr(:), iband1)
     155             :             elsewhere
     156             :             z_help = 0.0
     157             :             end where
     158             : 
     159           0 :             DO iband = bandi, bandf
     160           0 :                cprod_ir(iband, iband1, igpt) = 1/svol*dotprod(z_nk%data_c(:lapw%nv(jsp), iband), z_help)
     161             :             END DO !iband
     162             : 
     163             :          END DO  !iband1
     164             : 
     165             :       END DO !igpt
     166           0 :       CALL timestop("wavefproducts_noinv IR")
     167             : 
     168             :       !
     169             :       ! MT contribution
     170             :       !
     171           0 :       call read_cmt(cmt, nkqpt)
     172           0 :       lm_0 = 0
     173           0 :       ic = 0
     174             : 
     175           0 :       DO itype = 1, atoms%ntype
     176           0 :          DO ieq = 1, atoms%neq(itype)
     177           0 :             ic = ic + 1
     178           0 :             ic1 = 0
     179             : 
     180           0 :             cexp = exp(-2*img*pi_const*dot_product(kpts%bkf(:, iq), atoms%taual(:, ic)))
     181             : 
     182           0 :             DO l = 0, hybrid%lcutm1(itype)
     183           0 :                DO n = 1, hybdat%nindxp1(l, itype) ! loop over basis-function products
     184             : 
     185           0 :                   l1 = hybdat%prod(n, l, itype)%l1 !
     186           0 :                   l2 = hybdat%prod(n, l, itype)%l2 ! current basis-function product
     187           0 :                   n1 = hybdat%prod(n, l, itype)%n1 ! = bas(:,n1,l1,itype)*bas(:,n2,l2,itype) = b1*b2
     188           0 :                   n2 = hybdat%prod(n, l, itype)%n2 !
     189             : 
     190           0 :                   IF (mod(l1 + l2 + l, 2) /= 0) cycle
     191             : 
     192           0 :                   offdiag = l1 /= l2 .or. n1 /= n2 ! offdiag=true means that b1*b2 and b2*b1 are different combinations
     193             :                   !(leading to the same basis-function product)
     194             : 
     195           0 :                   lm1_0 = lmstart(l1, itype) ! start at correct lm index of cmt-coefficients
     196           0 :                   lm2_0 = lmstart(l2, itype) ! (corresponding to l1 and l2)
     197             : 
     198           0 :                   lm = lm_0
     199           0 :                   DO m = -l, l
     200             : 
     201           0 :                      carr = 0.0
     202             : 
     203           0 :                      lm1 = lm1_0 + n1 ! go to lm index for m1=-l1
     204           0 :                      DO m1 = -l1, l1
     205           0 :                         m2 = m1 + m ! Gaunt condition -m1+m2-m=0
     206           0 :                         IF (abs(m2) <= l2) THEN
     207           0 :                            lm2 = lm2_0 + n2 + (m2 + l2)*hybrid%nindx(l2, itype)
     208           0 :                            rdum = hybdat%gauntarr(1, l1, l2, l, m1, m) ! precalculated Gaunt coefficient
     209           0 :                            IF (rdum /= 0) THEN
     210           0 :                               DO iband = bandi, bandf
     211           0 :                                  cdum = rdum*conjg(cmt_nk(iband, lm1, ic)) !nk
     212           0 :                                  DO iband1 = 1, mnobd
     213           0 :                                     carr(iband1, iband) = carr(iband1, iband) + cdum*cmt(iband1, lm2, ic) !ikpt
     214             : 
     215             :                                  END DO
     216             :                               END DO
     217             :                            END IF
     218             :                         END IF
     219             : 
     220           0 :                         m2 = m1 - m ! switch role of b1 and b2
     221           0 :                         IF (abs(m2) <= l2 .and. offdiag) THEN
     222           0 :                            lm2 = lm2_0 + n2 + (m2 + l2)*hybrid%nindx(l2, itype)
     223           0 :                            rdum = hybdat%gauntarr(2, l1, l2, l, m1, m) ! precalculated Gaunt coefficient
     224           0 :                            IF (rdum /= 0) THEN
     225           0 :                               DO iband = bandi, bandf
     226           0 :                                  cdum = rdum*conjg(cmt_nk(iband, lm2, ic)) !nk
     227           0 :                                  DO iband1 = 1, mnobd
     228           0 :                                     carr(iband1, iband) = carr(iband1, iband) + cdum*cmt(iband1, lm1, ic)
     229             :                                  END DO
     230             :                               END DO
     231             :                            END IF
     232             :                         END IF
     233             : 
     234           0 :                         lm1 = lm1 + hybrid%nindx(l1, itype) ! go to lm start index for next m1-quantum number
     235             : 
     236             :                      END DO  !m1
     237             : 
     238           0 :                      DO iband = bandi, bandf
     239           0 :                         DO iband1 = 1, mnobd
     240           0 :                            cdum = carr(iband1, iband)*cexp
     241           0 :                            DO i = 1, hybrid%nindxm1(l, itype)
     242           0 :                               j = lm + i
     243           0 :                               cprod(j, iband1, iband) = cprod(j, iband1, iband) + hybdat%prodm(i, n, l, itype)*cdum
     244             :                            END DO
     245             : 
     246             :                         END DO
     247             :                      END DO
     248             : 
     249           0 :                      lm = lm + hybrid%nindxm1(l, itype) ! go to lm start index for next m-quantum number
     250             : 
     251             :                   END DO
     252             : 
     253             :                END DO
     254           0 :                lm_0 = lm_0 + hybrid%nindxm1(l, itype)*(2*l + 1) ! go to the lm start index of the next l-quantum number
     255           0 :                IF (lm /= lm_0) STOP 'wavefproducts: counting of lm-index incorrect (bug?)'
     256             :             END DO
     257             :          END DO
     258             :       END DO
     259             : 
     260           0 :       CALL timestop("wavefproducts_noinv")
     261             : 
     262           0 :       ic = nbasm_mt
     263           0 :       DO igpt = 1, hybrid%ngptm(iq)
     264           0 :          ic = ic + 1
     265           0 :          DO ibando = 1, mnobd
     266           0 :             DO iband = bandi, bandf
     267           0 :                cprod(ic, ibando, iband) = cprod_ir(iband, ibando, igpt)
     268             :             END DO
     269             :          END DO
     270             :       END DO
     271             : 
     272           0 :    END SUBROUTINE wavefproducts_noinv
     273             : 
     274           0 :    SUBROUTINE wavefproducts_inv(&
     275             :   &                  bandi, bandf, dimension, input, jsp, atoms,&
     276             :   &                  lapw, kpts,&
     277             :   &                  nk, iq, hybdat, mnobd, hybrid,&
     278             :   &                  parent, cell,&
     279             :   &                  nbasm_mt, sym, noco,&
     280           0 :   &                  nkqpt, cprod)
     281             : 
     282             :       USE m_util, ONLY: modulo1
     283             :       USE m_wrapper
     284             :       USE m_constants
     285             :       USE m_types
     286             :       USE m_io_hybrid
     287             :       IMPLICIT NONE
     288             :       TYPE(t_hybdat), INTENT(IN)   :: hybdat
     289             :       TYPE(t_dimension), INTENT(IN)   :: dimension
     290             :       TYPE(t_input), INTENT(IN)   :: input
     291             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
     292             :       TYPE(t_sym), INTENT(IN)   :: sym
     293             :       TYPE(t_noco), INTENT(IN)   :: noco
     294             :       TYPE(t_cell), INTENT(IN)   :: cell
     295             :       TYPE(t_kpts), INTENT(IN)   :: kpts
     296             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     297             :       TYPE(t_lapw), INTENT(IN)   :: lapw
     298             : 
     299             :       ! - scalars -
     300             :       INTEGER, INTENT(IN)      ::    bandi, bandf
     301             :       INTEGER, INTENT(IN)      ::    jsp, nk, iq
     302             :       INTEGER, INTENT(IN)      ::    mnobd
     303             :       INTEGER, INTENT(IN)      ::    nbasm_mt
     304             :       INTEGER, INTENT(OUT)     ::    nkqpt
     305             : 
     306             :       ! - arrays -
     307             :       INTEGER, INTENT(IN)      ::    parent(kpts%nkptf)
     308             : 
     309             :       REAL, INTENT(OUT)       ::    cprod(hybrid%maxbasm1, mnobd, bandf - bandi + 1)
     310             : 
     311             :       ! - local scalars -
     312             :       INTEGER                 ::    i, ikpt, ic, iband, iband1, igpt, igptp, ibando, iatom, iiatom, itype, ieq, ishift, ioffset, iatom1, iatom2
     313             :       INTEGER                 ::    l, p, l1, m1, l2, m2, p1, p2, n, ok
     314             :       INTEGER                 ::    lm, lm1, lm2, lm_0, lm_00, lm1_0, lm2_0, lmp1, lmp2, lmp3, lmp4, lp1, lp2
     315             :       INTEGER                 ::    j, ll, m, nbasfcn
     316             :       INTEGER                 :: nbasm_ir
     317             :       REAL                    ::    svol, sr2
     318             :       REAL                    ::    rdum, rfac, rfac1, rfac2, rdum1, rdum2
     319             :       REAL                    ::    sin1, sin2, cos1, cos2, add1, add2
     320             :       REAL                    ::    fac, fac1, fac2
     321             :       REAL                    ::    monepl1, monepl2, monepl, monepm1, monepm, moneplm, monepl1m1
     322             :       COMPLEX, PARAMETER       ::    img = (0.0, 1.0)
     323             :       COMPLEX                 ::    fexp
     324             :       COMPLEX                 ::    cdum, cconst, cfac
     325             :       LOGICAL                 ::    offdiag
     326             : 
     327             :       ! - local arrays -
     328           0 :       INTEGER                 ::    iarr(lapw%nv(jsp))
     329           0 :       INTEGER                 ::    gpt_nk(3, lapw%nv(jsp)), ghelp(3)
     330             :       INTEGER                 ::    gsum(3)
     331             :       INTEGER                 ::    g_t(3)
     332           0 :       INTEGER                 ::    lmstart(0:atoms%lmaxd, atoms%ntype)
     333             :       INTEGER                 ::    lmstart2(0:hybrid%maxlcutm1, atoms%nat)
     334             :       REAL                    ::    kqpt(3), kqpthlp(3)
     335             : 
     336           0 :       REAL, ALLOCATABLE        ::    cprod_ir(:, :, :)
     337             : 
     338           0 :       REAL                    ::    z_help(lapw%nv(jsp))
     339             : 
     340           0 :       REAL                    ::    cmt_nk(dimension%neigd, hybrid%maxlmindx, atoms%nat)
     341           0 :       REAL                    ::    cmt(dimension%neigd, hybrid%maxlmindx, atoms%nat)
     342             : 
     343           0 :       COMPLEX                 ::    ccmt_nk(dimension%neigd, hybrid%maxlmindx, atoms%nat)
     344           0 :       COMPLEX                 ::    ccmt(dimension%neigd, hybrid%maxlmindx, atoms%nat)
     345             : 
     346           0 :       REAL                    ::    rarr1(1:mnobd, bandf - bandi + 1)
     347           0 :       REAL                    ::    rarr(2, 1:mnobd, bandf - bandi + 1)
     348             :       COMPLEX                 ::    cmthlp(dimension%neigd), cmthlp1(dimension%neigd)
     349           0 :       COMPLEX                 ::    cexp(atoms%nat), cexp_nk(atoms%nat)
     350           0 :       TYPE(t_mat)             :: z_nk, z_kqpt
     351           0 :       TYPE(t_lapw)            :: lapw_nkqpt
     352             : 
     353           0 :       CALL timestart("wavefproducts_inv")
     354           0 :       CALL timestart("wavefproducts_inv IR")
     355           0 :       svol = sqrt(cell%omtil)
     356           0 :       sr2 = sqrt(2.0)
     357             : 
     358           0 :       nbasm_ir = maxval(hybrid%ngptm)
     359           0 :       ALLOCATE (cprod_ir(bandf - bandi + 1, mnobd, nbasm_ir))
     360           0 :       cprod_ir = 0
     361           0 :       gpt_nk(1, :) = lapw%k1(:lapw%nv(jsp), jsp)
     362           0 :       gpt_nk(2, :) = lapw%k2(:lapw%nv(jsp), jsp)
     363           0 :       gpt_nk(3, :) = lapw%k3(:lapw%nv(jsp), jsp)
     364             : 
     365             :       !
     366             :       ! compute k+q point for q (iq) in EIBZ(k)
     367             :       !
     368             : 
     369           0 :       kqpthlp = kpts%bkf(:, nk) + kpts%bkf(:, iq)
     370             :       ! kqpt can lie outside the first BZ, transfer it back
     371           0 :       kqpt = modulo1(kqpthlp, kpts%nkpt3)
     372           0 :       g_t(:) = nint(kqpt - kqpthlp)
     373             :       ! determine number of kqpt
     374           0 :       nkqpt = 0
     375           0 :       DO ikpt = 1, kpts%nkptf
     376           0 :          IF (maxval(abs(kqpt - kpts%bkf(:, ikpt))) <= 1E-06) THEN
     377           0 :             nkqpt = ikpt
     378           0 :             EXIT
     379             :          END IF
     380             :       END DO
     381           0 :       IF (nkqpt == 0) STOP 'wavefproducts: k-point not found'
     382             : 
     383             :       ! read in z at current k-point nk
     384             : 
     385           0 :       CALL lapw_nkqpt%init(input, noco, kpts, atoms, sym, nkqpt, cell, sym%zrfs)
     386           0 :       nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
     387           0 :       call z_nk%alloc(.true., nbasfcn, dimension%neigd)
     388           0 :       nbasfcn = MERGE(lapw_nkqpt%nv(1) + lapw_nkqpt%nv(2) + 2*atoms%nlotot, lapw_nkqpt%nv(1) + atoms%nlotot, noco%l_noco)
     389           0 :       call z_kqpt%alloc(.true., nbasfcn, dimension%neigd)
     390             : 
     391           0 :       call read_z(z_nk, nk)
     392           0 :       call read_z(z_kqpt, nkqpt)
     393             : 
     394           0 :       DO igpt = 1, hybrid%ngptm(iq)
     395           0 :          igptp = hybrid%pgptm(igpt, iq)
     396           0 :          ghelp = hybrid%gptm(:, igptp) - g_t(:)
     397           0 :          DO i = 1, lapw%nv(jsp)
     398           0 :             gsum(:) = ghelp + gpt_nk(:, i)
     399           0 :             IF (all(abs(gsum) <= hybdat%pntgptd)) THEN
     400           0 :                iarr(i) = hybdat%pntgpt(gsum(1), gsum(2), gsum(3), nkqpt)
     401             :             ELSE
     402           0 :                iarr(i) = 0
     403             :             END IF
     404             : 
     405             :          END DO
     406             : 
     407           0 :          DO iband1 = 1, hybrid%nobd(nkqpt)
     408           0 :             where (iarr > 0)
     409             :             z_help(:) = z_kqpt%data_r(iarr(:), iband1)
     410             :             elsewhere
     411             :             z_help = 0.0
     412             :             end where
     413           0 :             DO iband = bandi, bandf
     414           0 :                cprod_ir(iband, iband1, igpt) = 1/svol*dotprod(z_nk%data_r(:lapw%nv(jsp), iband), z_help)
     415             : 
     416             :             END DO !iband
     417             : 
     418             :          END DO  !iband1
     419             : 
     420             :       END DO !igpt
     421             : 
     422           0 :       CALL timestop("wavefproducts_inv IR")
     423             : 
     424             :       ! lmstart = lm start index for each l-quantum number and atom type (for cmt-coefficients)
     425           0 :       DO itype = 1, atoms%ntype
     426           0 :          DO l = 0, atoms%lmax(itype)
     427           0 :             lmstart(l, itype) = sum((/(hybrid%nindx(ll, itype)*(2*ll + 1), ll=0, l - 1)/))
     428             :          END DO
     429             :       END DO
     430             : 
     431             :       ! read in cmt coefficient at k-point nk
     432             : 
     433           0 :       CALL read_cmt(ccmt_nk, nk)
     434             : 
     435             :       !read in cmt coefficients at k+q point
     436           0 :       call read_cmt(ccmt, nkqpt)
     437             : 
     438           0 :       iatom = 0
     439           0 :       DO itype = 1, atoms%ntype
     440           0 :          DO ieq = 1, atoms%neq(itype)
     441           0 :             iatom = iatom + 1
     442             : 
     443           0 :             cexp(iatom) = exp((-img)*tpi_const*dotprod(kpts%bkf(:, iq) + kpts%bkf(:, nk), atoms%taual(:, iatom)))
     444             : 
     445           0 :             cexp_nk(iatom) = exp((-img)*tpi_const*dotprod(kpts%bkf(:, nk), atoms%taual(:, iatom)))
     446             :          END DO
     447             :       END DO
     448             : 
     449             :       rfac = 1./sr2
     450             :       cfac = -img/sr2
     451             :       iatom = 0
     452           0 :       DO itype = 1, atoms%ntype
     453           0 :          DO ieq = 1, atoms%neq(itype)
     454           0 :             iatom = iatom + 1
     455             :             ! determine the number of the inverse symmetric atom belonging to iatom
     456           0 :             IF (sym%invsatnr(iatom) == 0) THEN
     457             :                iiatom = iatom
     458             :             ELSE
     459           0 :                iiatom = sym%invsatnr(iatom)
     460             :             END IF
     461             :             ! the cmt coefficients at iatom and iiatom are made real in one step
     462           0 :             IF (iiatom < iatom) CYCLE
     463           0 :             lm1 = 0
     464           0 :             DO l = 0, atoms%lmax(itype)
     465           0 :                DO m = -l, l
     466           0 :                   DO p = 1, hybrid%nindx(l, itype)
     467           0 :                      lm1 = lm1 + 1
     468             :                      ! lm index at l,-m
     469           0 :                      lm2 = lm1 - 2*m*hybrid%nindx(l, itype)
     470             : 
     471           0 :                      IF (iatom == iiatom) THEN
     472           0 :                         IF (m < 0) THEN
     473           0 :                            cmt(:, lm1, iatom) = (ccmt(:, lm1, iatom) + (-1)**(l + m)*ccmt(:, lm2, iiatom))*cexp(iatom)*rfac
     474             : 
     475           0 :                            cmt_nk(:, lm1, iatom) = (ccmt_nk(:, lm1, iatom) + (-1)**(l + m)*ccmt_nk(:, lm2, iiatom))*cexp_nk(iatom)*rfac
     476           0 :                         ELSE IF (m > 0) THEN
     477             : 
     478           0 :                            cmt(:, lm1, iatom) = (ccmt(:, lm1, iatom) - (-1)**(l + m)*ccmt(:, lm2, iiatom))*cexp(iatom)*cfac
     479             : 
     480           0 :                            cmt_nk(:, lm1, iatom) = (ccmt_nk(:, lm1, iatom) - (-1)**(l + m)*ccmt_nk(:, lm2, iiatom))*cexp_nk(iatom)*cfac
     481             :                         ELSE
     482           0 :                            IF (mod(l, 2) == 0) THEN
     483           0 :                               cmt(:, lm1, iatom) = ccmt(:, lm1, iatom)*cexp(iatom)
     484           0 :                               cmt_nk(:, lm1, iatom) = ccmt_nk(:, lm1, iatom)*cexp_nk(iatom)
     485             :                            ELSE
     486           0 :                               cmt(:, lm1, iatom) = ccmt(:, lm1, iatom)*(-img)*cexp(iatom)
     487           0 :                               cmt_nk(:, lm1, iatom) = ccmt_nk(:, lm1, iatom)*(-img)*cexp_nk(iatom)
     488             :                            END IF
     489             :                         END IF
     490             :                      ELSE
     491           0 :                         cmt(:, lm1, iatom) = (ccmt(:, lm1, iatom) + (-1)**(l + m)*ccmt(:, lm2, iiatom))*rfac
     492             : 
     493           0 :                         cmt(:, lm1, iiatom) = (ccmt(:, lm1, iatom) - (-1)**(l + m)*ccmt(:, lm2, iiatom))*cfac
     494             : 
     495           0 :                         cmt_nk(:, lm1, iatom) = (ccmt_nk(:, lm1, iatom) + (-1)**(l + m)*ccmt_nk(:, lm2, iiatom))*rfac
     496             : 
     497           0 :                         cmt_nk(:, lm1, iiatom) = (ccmt_nk(:, lm1, iatom) - (-1)**(l + m)*ccmt_nk(:, lm2, iiatom))*cfac
     498             :                      END IF
     499             : 
     500             :                   END DO
     501             :                END DO
     502             :             END DO
     503             : 
     504             :          END DO
     505             :       END DO
     506             : 
     507           0 :       cprod = 0.0
     508             : 
     509           0 :       lm_0 = 0
     510           0 :       lm_00 = 0
     511           0 :       iatom1 = 0
     512           0 :       iiatom = 0
     513             : 
     514           0 :       DO itype = 1, atoms%ntype
     515           0 :          ioffset = sum((/((2*ll + 1)*hybrid%nindxm1(ll, itype), ll=0, hybrid%lcutm1(itype))/))
     516           0 :          lm_0 = lm_00
     517           0 :          DO ieq = 1, atoms%neq(itype)
     518           0 :             iatom1 = iatom1 + 1
     519           0 :             IF (sym%invsatnr(iatom1) == 0) THEN
     520             :                iatom2 = iatom1
     521             :             ELSE
     522           0 :                iatom2 = sym%invsatnr(iatom1)
     523             :             END IF
     524           0 :             IF (iatom1 > iatom2) CYCLE
     525             : 
     526           0 :             IF (iatom1 /= iatom2) THEN
     527             :                ! loop over l of mixed basis
     528           0 :                DO l = 0, hybrid%lcutm1(itype)
     529             :                   ! loop over basis functions products, which belong to l
     530           0 :                   DO n = 1, hybdat%nindxp1(l, itype)
     531             : 
     532             :                      ! determine l1,p1 and l2,p2 for the basis functions, which can generate l
     533           0 :                      l1 = hybdat%prod(n, l, itype)%l1
     534           0 :                      l2 = hybdat%prod(n, l, itype)%l2
     535           0 :                      p1 = hybdat%prod(n, l, itype)%n1
     536           0 :                      p2 = hybdat%prod(n, l, itype)%n2
     537             : 
     538             :                      ! condition for Gaunt coefficients
     539           0 :                      IF (mod(l + l1 + l2, 2) /= 0) CYCLE
     540             : 
     541           0 :                      offdiag = l1 /= l2 .or. p1 /= p2 ! offdiag=true means that b1*b2 and b2*b1 are different combinations
     542             :                      !(leading to the same basis-function product)
     543             : 
     544           0 :                      lm1_0 = lmstart(l1, itype) ! start at correct lm index of cmt-coefficients
     545           0 :                      lm2_0 = lmstart(l2, itype) ! (corresponding to l1 and l2)
     546             : 
     547           0 :                      lm = lm_0
     548           0 :                      lp1 = lm1_0 + p1
     549           0 :                      lp2 = lm2_0 + p2
     550             : 
     551             :                      ! sum over different m of mixed basis functions with qn l
     552           0 :                      DO m = -l, l
     553           0 :                         rarr = 0.0
     554             : 
     555             :                         ! go to lm index for m1=-l1
     556           0 :                         lmp1 = lm1_0 + p1
     557             : 
     558           0 :                         DO m1 = -l1, l1
     559             :                            ! Gaunt condition -m1+m2-m=0
     560           0 :                            m2 = m1 + m
     561           0 :                            IF (abs(m2) <= l2) THEN
     562           0 :                               lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
     563             :                               ! precalculated Gaunt coefficient
     564           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, m)
     565           0 :                               IF (rdum /= 0) THEN
     566           0 :                                  DO iband = bandi, bandf
     567           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)
     568           0 :                                     rdum2 = rdum*cmt_nk(iband, lmp1, iatom2)
     569             :                                     ! loop over occupied bands
     570           0 :                                     DO ibando = 1, mnobd!hybrid%nobd(peibz(ikpt))
     571             : 
     572           0 :                                        rarr(1, ibando, iband) = rarr(1, ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1) + rdum2*cmt(ibando, lmp2, iatom2)
     573             : 
     574           0 :                                        rarr(2, ibando, iband) = rarr(2, ibando, iband) + rdum1*cmt(ibando, lmp2, iatom2) - rdum2*cmt(ibando, lmp2, iatom1)
     575             : 
     576             :                                     END DO  !ibando
     577             :                                  END DO  !iband
     578             :                               END IF  ! rdum
     579             :                            END IF  ! abs(m2) .le. l2
     580             : 
     581           0 :                            m2 = m1 - m ! switch role of b1 and b2
     582           0 :                            IF (abs(m2) <= l2 .and. offdiag) THEN
     583           0 :                               lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
     584           0 :                               rdum = hybdat%gauntarr(2, l1, l2, l, m1, m) ! precalculated Gaunt coefficient
     585           0 :                               IF (rdum /= 0) THEN
     586           0 :                                  DO iband = bandi, bandf
     587           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)
     588           0 :                                     rdum2 = rdum*cmt_nk(iband, lmp2, iatom2)
     589             :                                     ! loop over occupied bands
     590           0 :                                     DO ibando = 1, mnobd!hybrid%nobd(peibz(ikpt)
     591           0 :                                        rarr(1, ibando, iband) = rarr(1, ibando, iband) + rdum1*cmt(ibando, lmp1, iatom1) + rdum2*cmt(ibando, lmp1, iatom2)
     592             : 
     593           0 :                                        rarr(2, ibando, iband) = rarr(2, ibando, iband) + rdum1*cmt(ibando, lmp1, iatom2) - rdum2*cmt(ibando, lmp1, iatom1)
     594             :                                     END DO  !ibando
     595             :                                  END DO  !iband
     596             :                               END IF  ! rdum .ne. 0
     597             :                            END IF  ! abs(m2) .le. l2 .and. offdiag
     598             : 
     599             :                            ! go to lmp start index for next m1-quantum number
     600           0 :                            lmp1 = lmp1 + hybrid%nindx(l1, itype)
     601             : 
     602             :                         END DO  !m1
     603             : 
     604           0 :                         ishift = -2*m*hybrid%nindxm1(l, itype)
     605             : 
     606             :                         ! go to lm mixed basis startindx for l and m
     607           0 :                         lm1 = lm + (iatom1 - 1 - iiatom)*ioffset
     608           0 :                         lm2 = lm + (iatom2 - 1 - iiatom)*ioffset + ishift
     609             : 
     610           0 :                         rdum = tpi_const*dotprod(kpts%bkf(:, iq), atoms%taual(:, iatom1))
     611           0 :                         rfac1 = sin(rdum)/sr2
     612           0 :                         rfac2 = cos(rdum)/sr2
     613           0 :                         DO iband = bandi, bandf
     614           0 :                            DO ibando = 1, mnobd
     615           0 :                               rdum1 = rarr(1, ibando, iband)
     616           0 :                               rdum2 = rarr(2, ibando, iband)
     617             : !                       sin1  = rdum1*rfac1
     618             : !                       cos1  = rdum1*rfac2
     619             : !                       sin2  = rdum2*rfac1
     620             : !                       cos2  = rdum2*rfac2
     621           0 :                               add1 = rdum1*rfac2 + rdum2*rfac1
     622           0 :                               add2 = rdum2*rfac2 - rdum1*rfac1
     623           0 :                               DO i = 1, hybrid%nindxm1(l, itype)
     624           0 :                                  j = lm1 + i
     625           0 :                                  cprod(j, ibando, iband) = cprod(j, ibando, iband) + hybdat%prodm(i, n, l, itype)*add1!( cos1 + sin2 )
     626           0 :                                  j = lm2 + i
     627           0 :                                  cprod(j, ibando, iband) = cprod(j, ibando, iband) + hybdat%prodm(i, n, l, itype)*add2!( cos2 - sin1 )
     628             : 
     629             :                               END DO  !i -> loop over mixed basis functions
     630             :                            END DO  !ibando
     631             :                         END DO  !iband
     632             : 
     633             :                         ! go to lm start index for next m-quantum number
     634           0 :                         lm = lm + hybrid%nindxm1(l, itype)
     635             : 
     636             :                      END DO  !m
     637             : 
     638             :                   END DO !n
     639           0 :                   lm_0 = lm_0 + hybrid%nindxm1(l, itype)*(2*l + 1) ! go to the lm start index of the next l-quantum number
     640           0 :                   IF (lm /= lm_0) STOP 'wavefproducts: counting of lm-index incorrect (bug?)'
     641             :                END DO !l
     642             : 
     643             :             ELSE !case: iatom1==iatom2
     644             : 
     645             :                ! loop over l of mixed basis
     646           0 :                monepl = -1
     647           0 :                DO l = 0, hybrid%lcutm1(itype)
     648           0 :                   monepl = -monepl
     649             :                   ! loop over basis functions products, which belong to l
     650           0 :                   DO n = 1, hybdat%nindxp1(l, itype)
     651             : 
     652             :                      ! determine l1,p1 and l2,p2 for the basis functions, which can generate l
     653           0 :                      l1 = hybdat%prod(n, l, itype)%l1
     654           0 :                      l2 = hybdat%prod(n, l, itype)%l2
     655           0 :                      p1 = hybdat%prod(n, l, itype)%n1
     656           0 :                      p2 = hybdat%prod(n, l, itype)%n2
     657             : 
     658             :                      ! condition for Gaunt coefficients
     659           0 :                      IF (mod(l + l1 + l2, 2) /= 0) CYCLE
     660             : 
     661           0 :                      offdiag = l1 /= l2 .or. p1 /= p2 ! offdiag=true means that b1*b2 and b2*b1 are different combinations
     662             :                      !(leading to the same basis-function product)
     663             : 
     664           0 :                      lm1_0 = lmstart(l1, itype) ! start at correct lm index of cmt-coefficients
     665           0 :                      lm2_0 = lmstart(l2, itype) ! (corresponding to l1 and l2)
     666             : 
     667           0 :                      lm = lm_0
     668           0 :                      lp1 = lm1_0 + p1
     669           0 :                      lp2 = lm2_0 + p2
     670             : 
     671             :                      ! calculate (-1)**l1 and (-1)**l2 (monepl = minus one power l)
     672           0 :                      monepl1 = (-1)**l1
     673           0 :                      monepl2 = (-1)**l2
     674             : 
     675             :                      ! sum over different m of mixed basis functions with qn l
     676             : 
     677             :                      !
     678             :                      !case m<0
     679             :                      !
     680             : 
     681           0 :                      monepm = -monepl
     682           0 :                      DO m = -l, -1
     683           0 :                         monepm = -monepm
     684           0 :                         moneplm = monepl*monepm
     685             : 
     686             :                         ! calculate the contributions which are identical for m>0 and m <0
     687           0 :                         rarr1 = 0.0
     688           0 :                         IF (abs(m) <= l2) THEN
     689           0 :                            lmp1 = lp1 + l1*hybrid%nindx(l1, itype)
     690           0 :                            IF (mod(l1, 2) == 0) THEN
     691           0 :                               lmp2 = lp2 + (m + l2)*hybrid%nindx(l2, itype)
     692             :                            ELSE
     693           0 :                               lmp2 = lp2 + (-m + l2)*hybrid%nindx(l2, itype)
     694             :                            END IF
     695             : 
     696           0 :                            rdum = hybdat%gauntarr(1, l1, l2, l, 0, m)
     697           0 :                            IF (rdum /= 0) THEN
     698           0 :                               DO iband = bandi, bandf
     699           0 :                                  rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)
     700           0 :                                  IF (mod(l1, 2) /= 0) rdum1 = moneplm*rdum1
     701           0 :                                  DO ibando = 1, mnobd
     702           0 :                                     rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
     703             :                                  END DO  ! ibando
     704             :                               END DO  ! iband
     705             :                            END IF  ! rdum .ne. 0
     706             : 
     707           0 :                            IF (offdiag) THEN
     708           0 :                               rdum = hybdat%gauntarr(1, l2, l1, l, -m, m)
     709           0 :                               IF (rdum /= 0) THEN
     710           0 :                                  DO iband = bandi, bandf
     711           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)
     712           0 :                                     IF (mod(l1, 2) == 0) rdum1 = moneplm*rdum1
     713           0 :                                     DO ibando = 1, mnobd
     714           0 :                                        rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp1, iatom1)
     715             :                                     END DO  ! ibando
     716             :                                  END DO  ! iband
     717             :                               END IF  ! rdum .ne. 0
     718             :                            END IF  ! offdiag
     719             : 
     720             :                         END IF  ! abs(m) .le. l2
     721             : 
     722           0 :                         IF (abs(m) <= l1) THEN
     723           0 :                            IF (mod(l2, 2) == 0) THEN
     724           0 :                               lmp3 = lp1 + (m + l1)*hybrid%nindx(l1, itype)
     725             :                            ELSE
     726           0 :                               lmp3 = lp1 + (-m + l1)*hybrid%nindx(l1, itype)
     727             :                            END IF
     728           0 :                            lmp2 = lp2 + l2*hybrid%nindx(l2, itype)
     729             : 
     730           0 :                            rdum = hybdat%gauntarr(1, l1, l2, l, -m, m)
     731           0 :                            IF (rdum /= 0) THEN
     732           0 :                               DO iband = bandi, bandf
     733           0 :                                  rdum1 = rdum*cmt_nk(iband, lmp3, iatom1)
     734           0 :                                  IF (mod(l2, 2) == 0) rdum1 = moneplm*rdum1
     735           0 :                                  DO ibando = 1, mnobd
     736           0 :                                     rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
     737             :                                  END DO  ! ibando
     738             :                               END DO  ! iband
     739             :                            END IF  ! rdum .ne. 0
     740             : 
     741           0 :                            IF (offdiag) THEN
     742           0 :                               rdum = hybdat%gauntarr(1, l2, l1, l, 0, m)
     743           0 :                               IF (rdum /= 0) THEN
     744           0 :                                  DO iband = bandi, bandf
     745           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)
     746           0 :                                     IF (mod(l2, 2) /= 0) rdum1 = moneplm*rdum1
     747           0 :                                     DO ibando = 1, mnobd
     748           0 :                                        rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
     749             :                                     END DO  ! ibando
     750             :                                  END DO  ! iband
     751             :                               END IF  ! rdum .ne. 0
     752             :                            END IF  ! offdiag
     753             : 
     754             :                         END IF  ! abs(m) .le. l2
     755             : 
     756             :                         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     757             : 
     758             :                         !go to lm index for m1=-l1
     759           0 :                         lmp1 = lp1
     760           0 :                         monepm1 = -monepl1
     761           0 :                         DO m1 = -l1, l1
     762           0 :                            monepm1 = -monepm1
     763           0 :                            IF (m1 == 0) THEN
     764           0 :                               lmp1 = lmp1 + hybrid%nindx(l1, itype)
     765           0 :                               CYCLE
     766             :                            END IF
     767             :                            ! (-1)**(l1+m1)
     768           0 :                            monepl1m1 = monepl1*monepm1
     769           0 :                            m2 = m1 + m
     770           0 :                            IF (abs(m2) <= l2 .and. m2 /= 0) THEN
     771           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, m)
     772           0 :                               IF (rdum /= 0) THEN
     773           0 :                                  IF (sign(1, m2) + sign(1, m1) /= 0) THEN
     774           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
     775             :                                  ELSE
     776           0 :                                     lmp2 = lp2 + (-m2 + l2)*hybrid%nindx(l2, itype)
     777           0 :                                     fac = 1/2.*moneplm*monepl1m1*(sign(1, m2) - sign(1, m1))
     778             :                                  END IF
     779           0 :                                  rdum = rdum/sr2
     780           0 :                                  DO iband = bandi, bandf
     781           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)!rdum*cmt_nk(iband,lmp1,iatom1)/sr2
     782           0 :                                     IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
     783           0 :                                     DO ibando = 1, mnobd
     784           0 :                                        rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
     785             :                                     END DO  ! ibando
     786             :                                  END DO  ! iband
     787             :                               END IF  ! rdum .ne. 0
     788             : 
     789           0 :                               IF (offdiag) THEN
     790           0 :                                  rdum = hybdat%gauntarr(1, l2, l1, l, m2, -m)
     791           0 :                                  IF (rdum /= 0) THEN
     792           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
     793           0 :                                     IF (sign(1, m2) + sign(1, m1) /= 0) THEN
     794             :                                        lmp3 = lmp1
     795             :                                     ELSE
     796           0 :                                        lmp3 = lmp1 - 2*m1*hybrid%nindx(l1, itype)
     797           0 :                                        fac = 1/2.*monepl1m1*(sign(1, m1) - sign(1, m2))
     798             :                                     END IF
     799           0 :                                     rdum = moneplm*rdum/sr2
     800           0 :                                     DO iband = bandi, bandf
     801           0 :                                        rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)!moneplm*rdum*cmt_nk(iband,lmp2,iatom1)/sr2
     802           0 :                                        IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
     803           0 :                                        DO ibando = 1, mnobd
     804           0 :                                           rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
     805             :                                        END DO  ! ibando
     806             :                                     END DO  ! iband
     807             :                                  END IF  ! rdum .ne. 0
     808             :                               END IF  ! offdiag
     809             : 
     810             :                            END IF  ! abs(m2) .le. l2 .and. m2 .ne. 0
     811             : 
     812           0 :                            m2 = m1 - m
     813           0 :                            IF (abs(m2) <= l2 .and. m2 /= 0) THEN
     814             : 
     815           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, -m)
     816           0 :                               IF (rdum /= 0) THEN
     817             : 
     818           0 :                                  IF (sign(1, m2) + sign(1, m1) /= 0) THEN
     819           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
     820             :                                  ELSE
     821           0 :                                     lmp2 = lp2 + (-m2 + l2)*hybrid%nindx(l2, itype)
     822           0 :                                     fac = 1/2.*moneplm*monepl1m1*(sign(1, m2) - sign(1, m1))
     823             :                                  END IF
     824           0 :                                  rdum = moneplm*rdum/sr2
     825           0 :                                  DO iband = bandi, bandf
     826           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)!moneplm*rdum*cmt_nk(iband,lmp1,iatom1)/sr2
     827           0 :                                     IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
     828           0 :                                     DO ibando = 1, mnobd
     829           0 :                                        rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
     830             :                                     END DO  ! ibando
     831             :                                  END DO  ! iband
     832             : 
     833             :                               END IF  ! rdum .ne. 0
     834             : 
     835           0 :                               IF (offdiag) THEN
     836           0 :                                  rdum = hybdat%gauntarr(1, l2, l1, l, m2, m)
     837           0 :                                  IF (rdum /= 0) THEN
     838           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
     839           0 :                                     IF (sign(1, m1) + sign(1, m2) /= 0) THEN
     840             :                                        lmp3 = lmp1
     841             :                                     ELSE
     842           0 :                                        lmp3 = lmp1 - 2*m1*hybrid%nindx(l1, itype)
     843           0 :                                        fac = 1/2.*monepl1m1*(sign(1, m1) - sign(1, m2))
     844             :                                     END IF
     845           0 :                                     rdum = rdum/sr2
     846           0 :                                     DO iband = bandi, bandf
     847           0 :                                        rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)!rdum*cmt_nk(iband,lmp2,iatom1)/sr2
     848           0 :                                        IF (sign(1, m1) + sign(1, m2) == 0) rdum1 = fac*rdum1
     849           0 :                                        DO ibando = 1, mnobd
     850           0 :                                           rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
     851             :                                        END DO  ! ibando
     852             :                                     END DO  ! iband
     853             :                                  END IF  ! rdum .ne. 0
     854             :                               END IF  ! offdiag
     855             : 
     856             :                            END IF  ! abs(m2) .le. l2 .and. m1 .ne. 0
     857             : 
     858             :                            !go to lmp start index for next m1-quantum number
     859           0 :                            lmp1 = lmp1 + hybrid%nindx(l1, itype)
     860             :                         END DO  ! m1
     861             : 
     862             :                         ! go to lm mixed basis startindx for l and m
     863           0 :                         lm1 = lm + (iatom1 - 1 - iiatom)*ioffset
     864           0 :                         DO iband = bandi, bandf
     865           0 :                            DO ibando = 1, mnobd
     866           0 :                               rdum = rarr1(ibando, iband)
     867           0 :                               DO i = 1, hybrid%nindxm1(l, itype)
     868           0 :                                  j = lm1 + i
     869           0 :                                  cprod(j, ibando, iband) = cprod(j, ibando, iband) + hybdat%prodm(i, n, l, itype)*rdum
     870             :                               END DO  !i -> loop over mixed basis functions
     871             :                            END DO  !ibando
     872             :                         END DO  !iband
     873             : 
     874             :                         ! go to lm start index for next m-quantum number
     875           0 :                         lm = lm + hybrid%nindxm1(l, itype)
     876             : 
     877             :                      END DO  ! m=-l,-1
     878             : 
     879             :                      !
     880             :                      !case m=0
     881             :                      !
     882             : 
     883             :                      m = 0
     884           0 :                      rarr1 = 0.0
     885           0 :                      lmp1 = lp1
     886             : 
     887           0 :                      monepm1 = -monepl1
     888             : 
     889           0 :                      DO m1 = -l1, l1
     890           0 :                         m2 = m1
     891           0 :                         monepm1 = -monepm1
     892           0 :                         IF (abs(m2) <= l2) THEN
     893             : 
     894           0 :                            IF (mod(l, 2) == 0) THEN
     895           0 :                               lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
     896             :                               !lmp3 and lmp4 are variables, which avoid an if clause in the loop
     897           0 :                               lmp3 = lmp2
     898           0 :                               lmp4 = lmp1
     899             :                            ELSE
     900           0 :                               lmp2 = lp2 + (-m2 + l2)*hybrid%nindx(l2, itype)
     901             :                               !lmp3 and lmp3 are variables, which avoid an if clause in the loop
     902           0 :                               lmp3 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
     903           0 :                               lmp4 = lmp1 - 2*m1*hybrid%nindx(l1, itype)
     904             : 
     905           0 :                               fac1 = monepl1*monepm1 ! (-1)**(l1+m1)
     906           0 :                               fac2 = monepl2*monepm1 ! (-1)**(l2+m1)
     907             :                            END IF
     908             : 
     909             :                            !precalculated Gaunt coefficient
     910           0 :                            IF (mod(l, 2) == 0) THEN
     911           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, m)
     912             :                            ELSE
     913           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, m)*fac1
     914             :                            END IF
     915           0 :                            IF (rdum /= 0) THEN
     916           0 :                               DO iband = bandi, bandf
     917           0 :                                  rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)
     918           0 :                                  DO ibando = 1, mnobd
     919           0 :                                     rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
     920             :                                  END DO  ! ibando
     921             :                               END DO  ! iband
     922             :                            END IF  ! rdum.ne.0
     923             : 
     924             :                            !change role of b1 and b2
     925           0 :                            IF (offdiag) THEN
     926           0 :                               IF (mod(l, 2) == 0) THEN
     927           0 :                                  rdum = hybdat%gauntarr(2, l1, l2, l, m1, m)
     928             :                               ELSE
     929           0 :                                  rdum = hybdat%gauntarr(2, l1, l2, l, m1, m)*fac2
     930             :                               END IF
     931           0 :                               IF (rdum /= 0) THEN
     932           0 :                                  DO iband = bandi, bandf
     933           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp3, iatom1)
     934           0 :                                     DO ibando = 1, mnobd
     935           0 :                                        rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp4, iatom1)
     936             :                                     END DO  ! ibando
     937             :                                  END DO  ! iband
     938             :                               END IF  ! rdum.ne.0
     939             :                            END IF  ! offdiag
     940             : 
     941             :                         END IF  ! abs(m2).le.l2
     942             : 
     943             :                         ! go to lmp start index for next m1-quantum number
     944           0 :                         lmp1 = lmp1 + hybrid%nindx(l1, itype)
     945             :                      END DO  !m1
     946             : 
     947             :                      ! go to lm mixed basis startindx for l and m
     948           0 :                      lm1 = lm + (iatom1 - 1 - iiatom)*ioffset
     949           0 :                      DO iband = bandi, bandf
     950           0 :                         DO ibando = 1, mnobd
     951           0 :                            rdum = rarr1(ibando, iband)
     952           0 :                            DO i = 1, hybrid%nindxm1(l, itype)
     953           0 :                               j = lm1 + i
     954           0 :                               cprod(j, ibando, iband) = cprod(j, ibando, iband) + hybdat%prodm(i, n, l, itype)*rdum
     955             :                            END DO  !i -> loop over mixed basis functions
     956             :                         END DO  !ibando
     957             :                      END DO  !iband
     958             : 
     959             :                      ! go to lm start index for next m-quantum number
     960           0 :                      lm = lm + hybrid%nindxm1(l, itype)
     961             : 
     962             :                      !
     963             :                      ! case: m>0
     964             :                      !
     965             : 
     966           0 :                      rarr1 = 0.0
     967             :                      monepm = 1
     968           0 :                      DO m = 1, l
     969           0 :                         monepm = -monepm
     970           0 :                         moneplm = monepl*monepm
     971             : 
     972             :                         ! calculate the contributions which are identical for m>0 and m <0
     973           0 :                         rarr1 = 0.0
     974           0 :                         IF (abs(m) <= l2) THEN
     975           0 :                            lmp1 = lp1 + l1*hybrid%nindx(l1, itype)
     976           0 :                            IF (mod(l1, 2) == 0) THEN
     977           0 :                               lmp2 = lp2 + (m + l2)*hybrid%nindx(l2, itype)
     978             :                            ELSE
     979           0 :                               lmp2 = lp2 + (-m + l2)*hybrid%nindx(l2, itype)
     980             :                            END IF
     981             : 
     982           0 :                            rdum = hybdat%gauntarr(1, l1, l2, l, 0, m)
     983           0 :                            IF (rdum /= 0) THEN
     984           0 :                               DO iband = bandi, bandf
     985           0 :                                  rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)
     986           0 :                                  IF (mod(l1, 2) /= 0) rdum1 = moneplm*rdum1
     987           0 :                                  DO ibando = 1, mnobd
     988           0 :                                     rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
     989             :                                  END DO  ! ibando
     990             :                               END DO  ! iband
     991             :                            END IF  ! rdum .ne. 0
     992             : 
     993           0 :                            IF (offdiag) THEN
     994           0 :                               rdum = hybdat%gauntarr(1, l2, l1, l, -m, m)
     995           0 :                               IF (rdum /= 0) THEN
     996           0 :                                  DO iband = bandi, bandf
     997           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)
     998           0 :                                     IF (mod(l1, 2) == 0) rdum1 = moneplm*rdum1
     999           0 :                                     DO ibando = 1, mnobd
    1000           0 :                                        rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp1, iatom1)
    1001             :                                     END DO  ! ibando
    1002             :                                  END DO  ! iband
    1003             :                               END IF  ! rdum .ne. 0
    1004             :                            END IF  ! offdiag
    1005             : 
    1006             :                         END IF  ! abs(m) .le. l2
    1007             : 
    1008           0 :                         IF (abs(m) <= l1) THEN
    1009           0 :                            IF (mod(l2, 2) == 0) THEN
    1010           0 :                               lmp3 = lp1 + (m + l1)*hybrid%nindx(l1, itype)
    1011             :                            ELSE
    1012           0 :                               lmp3 = lp1 + (-m + l1)*hybrid%nindx(l1, itype)
    1013             :                            END IF
    1014           0 :                            lmp2 = lp2 + l2*hybrid%nindx(l2, itype)
    1015             : 
    1016           0 :                            rdum = hybdat%gauntarr(1, l1, l2, l, -m, m)
    1017           0 :                            IF (rdum /= 0) THEN
    1018           0 :                               DO iband = bandi, bandf
    1019           0 :                                  rdum1 = rdum*cmt_nk(iband, lmp3, iatom1)
    1020           0 :                                  IF (mod(l2, 2) == 0) rdum1 = moneplm*rdum1
    1021           0 :                                  DO ibando = 1, mnobd
    1022           0 :                                     rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    1023             :                                  END DO  ! ibando
    1024             :                               END DO  ! iband
    1025             :                            END IF  ! rdum .ne. 0
    1026             : 
    1027           0 :                            IF (offdiag) THEN
    1028           0 :                               rdum = hybdat%gauntarr(1, l2, l1, l, 0, m)
    1029           0 :                               IF (rdum /= 0) THEN
    1030           0 :                                  DO iband = bandi, bandf
    1031           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)
    1032           0 :                                     IF (mod(l2, 2) /= 0) rdum1 = moneplm*rdum1
    1033           0 :                                     DO ibando = 1, mnobd
    1034           0 :                                        rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
    1035             :                                     END DO  ! ibando
    1036             :                                  END DO  ! iband
    1037             :                               END IF  ! rdum .ne. 0
    1038             :                            END IF  ! offdiag
    1039             : 
    1040             :                         END IF  ! abs(m) .le. l2
    1041             : 
    1042             :                         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1043             : 
    1044             :                         !go to lm index for m1=-l1
    1045             :                         lmp1 = lp1
    1046             :                         monepm1 = -monepl1
    1047           0 :                         DO m1 = -l1, l1
    1048           0 :                            monepm1 = -monepm1
    1049           0 :                            IF (m1 == 0) THEN
    1050           0 :                               lmp1 = lmp1 + hybrid%nindx(l1, itype)
    1051           0 :                               CYCLE
    1052             :                            END IF
    1053           0 :                            m2 = m1 + m
    1054             :                            ! (-1)**(l1+m1)
    1055           0 :                            monepl1m1 = monepl1*monepm1
    1056           0 :                            IF (abs(m2) <= l2 .and. m2 /= 0) THEN
    1057           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, m)
    1058           0 :                               IF (rdum /= 0) THEN
    1059             : 
    1060           0 :                                  IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    1061           0 :                                     lmp2 = lp2 + (-m2 + l2)*hybrid%nindx(l2, itype)
    1062             :                                  ELSE
    1063           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1064           0 :                                     fac = -moneplm*monepl1m1*(sign(1, m2) - sign(1, m1))/2
    1065             :                                  END IF
    1066             : 
    1067           0 :                                  rdum = -moneplm*monepl1m1*rdum/sr2
    1068           0 :                                  DO iband = bandi, bandf
    1069           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)!-moneplm*monepl1m1*rdum*cmt_nk(iband,lmp1,iatom1)/sr2
    1070           0 :                                     IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
    1071           0 :                                     DO ibando = 1, mnobd
    1072           0 :                                        rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    1073             :                                     END DO  ! ibando
    1074             :                                  END DO  ! iband
    1075             : 
    1076             :                               END IF  ! rdum .ne. 0
    1077             : 
    1078           0 :                               IF (offdiag) THEN
    1079           0 :                                  rdum = hybdat%gauntarr(2, l1, l2, l, m1, -m)
    1080           0 :                                  IF (rdum /= 0) THEN
    1081           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1082           0 :                                     IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    1083           0 :                                        lmp3 = lmp1 - 2*m1*hybrid%nindx(l1, itype)
    1084             :                                     ELSE
    1085           0 :                                        lmp3 = lmp1
    1086           0 :                                        fac = 1/2.*monepl1m1*(sign(1, m2) - sign(1, m1))
    1087             :                                     END IF
    1088           0 :                                     rdum = monepl1m1*moneplm*rdum/sr2
    1089           0 :                                     DO iband = bandi, bandf
    1090           0 :                                        rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)!monepl1m1*moneplm*rdum*cmt_nk(iband,lmp2,iatom1)/sr2
    1091           0 :                                        IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
    1092           0 :                                        DO ibando = 1, mnobd
    1093           0 :                                           rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
    1094             :                                        END DO  ! ibando
    1095             :                                     END DO  ! iband
    1096             : 
    1097             :                                  END IF  ! rdum
    1098             : 
    1099             :                               END IF  ! offdiag
    1100             :                            END IF  ! abs(m2) .le. l2 .and. m2 .ne. 0
    1101             : 
    1102           0 :                            m2 = m1 - m
    1103           0 :                            IF (abs(m2) <= l2 .and. m2 /= 0) THEN
    1104             : 
    1105           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, -m)
    1106           0 :                               IF (rdum /= 0) THEN
    1107             : 
    1108           0 :                                  IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    1109           0 :                                     lmp2 = lp2 + (-m2 + l2)*hybrid%nindx(l2, itype)
    1110             :                                  ELSE
    1111           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1112           0 :                                     fac = 1/2.*moneplm*monepl1m1*(sign(1, m1) - sign(1, m2))
    1113             :                                  END IF
    1114           0 :                                  rdum = monepl1m1*rdum/sr2
    1115           0 :                                  DO iband = bandi, bandf
    1116           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)!monepl1m1*rdum*cmt_nk(iband,lmp1,iatom1)/sr2
    1117           0 :                                     IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = rdum1*fac
    1118           0 :                                     DO ibando = 1, mnobd
    1119           0 :                                        rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    1120             :                                     END DO  ! ibando
    1121             :                                  END DO  ! iband
    1122             :                               END IF  ! rdum .ne. 0
    1123             : 
    1124           0 :                               IF (offdiag) THEN
    1125           0 :                                  rdum = hybdat%gauntarr(2, l1, l2, l, m1, m)
    1126           0 :                                  IF (rdum /= 0) THEN
    1127           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1128           0 :                                     IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    1129           0 :                                        lmp3 = lmp1 - 2*m1*hybrid%nindx(l1, itype)
    1130             :                                     ELSE
    1131           0 :                                        lmp3 = lmp1
    1132           0 :                                        fac = -monepl1m1*(sign(1, m1) - sign(1, m2))/2
    1133             :                                     END IF
    1134           0 :                                     rdum = -monepl1m1*rdum/sr2
    1135           0 :                                     DO iband = bandi, bandf
    1136           0 :                                        rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)!-monepl1m1*rdum*cmt_nk(iband,lmp2,iatom1)/sr2
    1137           0 :                                        IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
    1138           0 :                                        DO ibando = 1, mnobd
    1139           0 :                                           rarr1(ibando, iband) = rarr1(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
    1140             :                                        END DO  ! ibando
    1141             :                                     END DO  ! iband
    1142             : 
    1143             :                                  END IF  ! rdum .ne. 0
    1144             :                               END IF  ! offdiag
    1145             : 
    1146             :                            END IF  !  abs(m2) .le. l2 .and. m2 .ne. 0
    1147             : 
    1148             :                            !go to lmp start index for next m1-quantum number
    1149           0 :                            lmp1 = lmp1 + hybrid%nindx(l1, itype)
    1150             :                         END DO  ! m1
    1151             : 
    1152             :                         ! multiply rarr1 by (-1)**(l+m+1)
    1153           0 :                         rarr1(:, :) = (-1)*moneplm*rarr1(:, :)
    1154             : 
    1155             :                         ! go to lm mixed basis startindx for l and m
    1156           0 :                         lm1 = lm + (iatom1 - 1 - iiatom)*ioffset
    1157             : 
    1158           0 :                         DO iband = bandi, bandf
    1159           0 :                            DO ibando = 1, mnobd
    1160           0 :                               rdum = rarr1(ibando, iband)
    1161           0 :                               DO i = 1, hybrid%nindxm1(l, itype)
    1162           0 :                                  j = lm1 + i
    1163           0 :                                  cprod(j, ibando, iband) = cprod(j, ibando, iband) + hybdat%prodm(i, n, l, itype)*rdum
    1164             :                               END DO  !i -> loop over mixed basis functions
    1165             :                            END DO  !ibando
    1166             :                         END DO  !iband
    1167             : 
    1168             :                         ! go to lm start index for next m-quantum number
    1169           0 :                         lm = lm + hybrid%nindxm1(l, itype)
    1170             : 
    1171             :                      END DO  ! m=1,l
    1172             : 
    1173             :                   END DO !n
    1174           0 :                   lm_0 = lm_0 + hybrid%nindxm1(l, itype)*(2*l + 1) ! go to the m start index of the next l-quantum number
    1175           0 :                   IF (lm /= lm_0) STOP 'wavefproducts: counting of lm-index incorrect (bug?)'
    1176             :                END DO !l
    1177             : 
    1178             :             END IF  ! iatom1 .ne. iatom2
    1179             : 
    1180           0 :             lm_0 = lm_00
    1181             :          END DO !ieq
    1182           0 :          iiatom = iiatom + atoms%neq(itype)
    1183           0 :          lm_00 = lm_00 + atoms%neq(itype)*ioffset
    1184             :       END DO  !itype
    1185             : 
    1186           0 :       ic = nbasm_mt
    1187           0 :       DO igpt = 1, hybrid%ngptm(iq)
    1188           0 :          ic = ic + 1
    1189           0 :          DO ibando = 1, mnobd
    1190           0 :             DO iband = bandi, bandf
    1191           0 :                cprod(ic, ibando, iband) = cprod_ir(iband, ibando, igpt)
    1192             :             END DO
    1193             :          END DO
    1194             :       END DO
    1195             : 
    1196           0 :       CALL timestop("wavefproducts_inv")
    1197             : 
    1198           0 :    END SUBROUTINE wavefproducts_inv
    1199             : 
    1200           0 :    SUBROUTINE wavefproducts_inv5(&
    1201             :   &                    bandi, bandf, bandoi, bandof,&
    1202             :   &                    dimension, input, jsp, atoms,&
    1203             :   &                    lapw, kpts,&
    1204             :   &                    nk, iq, hybdat, mnobd, hybrid,&
    1205             :   &                    parent, cell,&
    1206             :   &                    nbasm_mt, sym,&
    1207             :   &                    noco,&
    1208           0 :   &                    nkqpt, cprod)
    1209             : 
    1210             :       USE m_util, ONLY: modulo1
    1211             :       USE m_olap, ONLY: gptnorm
    1212             :       USE m_wrapper
    1213             :       USE m_constants
    1214             :       USE m_types
    1215             :       USE m_io_hybrid
    1216             :       IMPLICIT NONE
    1217             :       TYPE(t_dimension), INTENT(IN) :: dimension
    1218             :       TYPE(t_hybrid), INTENT(IN)    :: hybrid
    1219             :       TYPE(t_input), INTENT(IN)     :: input
    1220             :       TYPE(t_noco), INTENT(IN)      :: noco
    1221             :       TYPE(t_sym), INTENT(IN)       :: sym
    1222             :       TYPE(t_cell), INTENT(IN)      :: cell
    1223             :       TYPE(t_kpts), INTENT(IN)      :: kpts
    1224             :       TYPE(t_atoms), INTENT(IN)     :: atoms
    1225             :       TYPE(t_lapw), INTENT(IN)      :: lapw
    1226             :       TYPE(t_hybdat), INTENT(INOUT) :: hybdat
    1227             : 
    1228             :       ! - scalars -
    1229             :       INTEGER, INTENT(IN)      :: bandi, bandf, bandoi, bandof
    1230             :       INTEGER, INTENT(IN)      :: jsp, nk, iq
    1231             :       INTEGER, INTENT(IN)      :: mnobd
    1232             :       INTEGER, INTENT(IN)      :: nbasm_mt
    1233             :       INTEGER, INTENT(OUT)     :: nkqpt
    1234             : 
    1235             :       ! - arrays -
    1236             :       INTEGER, INTENT(IN)      ::    parent(kpts%nkptf)
    1237             : 
    1238             :       REAL, INTENT(OUT)        ::    cprod(hybrid%maxbasm1, bandoi:bandof, bandf - bandi + 1)
    1239             : 
    1240             :       ! - local scalars -
    1241             :       INTEGER                 ::    i, ikpt, ic, iband, iband1, igpt, igptp, ig, ig2, ig1
    1242             :       INTEGER                 ::    iatom, iiatom, itype, ieq, ishift
    1243             :       INTEGER                 ::    ibando, iatom1, iatom2, ioffset
    1244             :       INTEGER                 ::    k, l, p, l1, m1, l2, m2, p1, p2, n, ok
    1245             :       INTEGER                 ::    igptm, iigptm
    1246             :       INTEGER                 ::    lm, lm1, lm2, lm_0, lm_00, lm1_0, lm2_0, lmp1, lmp2, lmp3, lmp4, lp1, lp2
    1247             :       INTEGER                 ::    j, ll, m
    1248             :       INTEGER                 ::    nbasm_ir, ngpt0
    1249             :       INTEGER                 ::    n1, n2, nbasfcn
    1250             :       REAL                    ::    svol, sr2
    1251             :       REAL                    ::    rdum, rfac, rfac1, rfac2, rdum1, rdum2
    1252             :       REAL                    ::    sin1, sin2, cos1, cos2, add1, add2
    1253             :       REAL                    ::    fac, fac1, fac2
    1254             :       REAL                    ::    monepl1, monepl2, monepl, monepm1, monepm, moneplm, monepl1m1
    1255             :       COMPLEX                 ::    fexp
    1256             :       COMPLEX                 ::    cdum, cconst, cfac
    1257             :       COMPLEX, PARAMETER       ::    img = (0.0, 1.0)
    1258             :       LOGICAL                 ::    offdiag
    1259           0 :       TYPE(t_lapw)            ::    lapw_nkqpt
    1260             : 
    1261             :       ! - local arrays -
    1262             :       INTEGER                 ::    g(3), g_t(3)
    1263           0 :       INTEGER                 ::    lmstart(0:atoms%lmaxd, atoms%ntype)
    1264           0 :       INTEGER, ALLOCATABLE    ::    gpt0(:, :)
    1265           0 :       INTEGER, ALLOCATABLE    ::    pointer(:, :, :)
    1266             : 
    1267             :       REAL                    ::    kqpt(3), kqpthlp(3)
    1268             :       REAL                    ::    bkpt(3)
    1269           0 :       REAL                    ::    cmt_nk(dimension%neigd, hybrid%maxlmindx, atoms%nat)
    1270           0 :       REAL                    ::    cmt(dimension%neigd, hybrid%maxlmindx, atoms%nat)
    1271           0 :       REAL                    ::    rarr1(bandoi:bandof)
    1272           0 :       REAL                    ::    rarr2(bandoi:bandof, bandf - bandi + 1)
    1273           0 :       REAL                    ::    rarr3(2, bandoi:bandof, bandf - bandi + 1)
    1274           0 :       REAL, ALLOCATABLE       ::    z0(:, :)
    1275             : 
    1276           0 :       COMPLEX                 ::    cexp(atoms%nat), cexp_nk(atoms%nat)
    1277           0 :       COMPLEX, ALLOCATABLE     ::    ccmt_nk(:, :, :)
    1278           0 :       COMPLEX, ALLOCATABLE     ::    ccmt(:, :, :)
    1279           0 :       TYPE(t_mat)             :: z_nk, z_kqpt
    1280             : 
    1281           0 :       CALL timestart("wavefproducts_inv5")
    1282           0 :       CALL timestart("wavefproducts_inv5 IR")
    1283             : 
    1284           0 :       cprod = 0
    1285           0 :       svol = sqrt(cell%omtil)
    1286           0 :       sr2 = sqrt(2.0)
    1287             : 
    1288           0 :       nbasm_ir = maxval(hybrid%ngptm)
    1289             : 
    1290             :       !
    1291             :       ! compute k+q point for q (iq) in EIBZ(k)
    1292             :       !
    1293             : 
    1294           0 :       kqpthlp = kpts%bkf(:, nk) + kpts%bkf(:, iq)
    1295             :       ! kqpt can lie outside the first BZ, transfer it back
    1296           0 :       kqpt = modulo1(kqpthlp, kpts%nkpt3)
    1297           0 :       g_t(:) = nint(kqpt - kqpthlp)
    1298             :       ! determine number of kqpt
    1299           0 :       nkqpt = 0
    1300           0 :       DO ikpt = 1, kpts%nkptf
    1301           0 :          IF (maxval(abs(kqpt - kpts%bkf(:, ikpt))) <= 1E-06) THEN
    1302           0 :             nkqpt = ikpt
    1303           0 :             EXIT
    1304             :          END IF
    1305             :       END DO
    1306           0 :       IF (nkqpt == 0) STOP 'wavefproducts_inv5: k-point not found'
    1307             : 
    1308             :       !
    1309             :       ! compute G's fulfilling |bk(:,nkqpt) + G| <= rkmax
    1310             :       !
    1311           0 :       CALL lapw_nkqpt%init(input, noco, kpts, atoms, sym, nkqpt, cell, sym%zrfs)
    1312             : 
    1313           0 :       nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
    1314           0 :       call z_nk%alloc(.true., nbasfcn, dimension%neigd)
    1315           0 :       nbasfcn = MERGE(lapw_nkqpt%nv(1) + lapw_nkqpt%nv(2) + 2*atoms%nlotot, lapw_nkqpt%nv(1) + atoms%nlotot, noco%l_noco)
    1316           0 :       call z_kqpt%alloc(.true., nbasfcn, dimension%neigd)
    1317             : 
    1318             :       ! read in z at k-point nk and nkqpt
    1319           0 :       call timestart("read_z")
    1320           0 :       CALL read_z(z_nk, nk)
    1321           0 :       call read_z(z_kqpt, nkqpt)
    1322           0 :       call timestop("read_z")
    1323             : 
    1324             :       g(1) = maxval(abs(lapw%k1(:lapw%nv(jsp), jsp))) &
    1325             :      &     + maxval(abs(lapw_nkqpt%k1(:lapw_nkqpt%nv(jsp), jsp)))&
    1326           0 :      &     + maxval(abs(hybrid%gptm(1, hybrid%pgptm(:hybrid%ngptm(iq), iq)))) + 1
    1327             :       g(2) = maxval(abs(lapw%k2(:lapw%nv(jsp), jsp)))&
    1328             :      &     + maxval(abs(lapw_nkqpt%k2(:lapw_nkqpt%nv(jsp), jsp)))&
    1329           0 :      &     + maxval(abs(hybrid%gptm(2, hybrid%pgptm(:hybrid%ngptm(iq), iq)))) + 1
    1330             :       g(3) = maxval(abs(lapw%k3(:lapw%nv(jsp), jsp)))&
    1331             :      &     + maxval(abs(lapw_nkqpt%k3(:lapw_nkqpt%nv(jsp), jsp)))&
    1332           0 :      &     + maxval(abs(hybrid%gptm(3, hybrid%pgptm(:hybrid%ngptm(iq), iq)))) + 1
    1333             : 
    1334           0 :       ALLOCATE (pointer(-g(1):g(1), -g(2):g(2), -g(3):g(3)), stat=ok)
    1335           0 :       IF (ok /= 0) STOP 'wavefproducts_inv5: error allocation pointer'
    1336           0 :       ALLOCATE (gpt0(3, size(pointer)), stat=ok)
    1337           0 :       IF (ok /= 0) STOP 'wavefproducts_inv5: error allocation gpt0'
    1338             : 
    1339           0 :       if (.not. allocated(hybdat%stepfunc_r)) then
    1340           0 :          call timestart("setup stepfunction")
    1341           0 :          ALLOCATE (hybdat%stepfunc_r(-g(1):g(1), -g(2):g(2), -g(3):g(3)), stat=ok)
    1342           0 :          IF (ok /= 0) then
    1343           0 :             call juDFT_error('wavefproducts_inv5: error allocation stepfunc_r')
    1344             :          endif
    1345             : 
    1346           0 :          DO i = -g(1), g(1)
    1347           0 :             DO j = -g(2), g(2)
    1348           0 :                DO k = -g(3), g(3)
    1349           0 :                   hybdat%stepfunc_r(i, j, k) = stepfunction(cell, atoms, (/i, j, k/))
    1350             :                END DO
    1351             :             END DO
    1352             :          END DO
    1353           0 :          call timestop("setup stepfunction")
    1354             :       endif
    1355             : 
    1356             :       !
    1357             :       ! convolute phi(n,k) with the step function and store in cpw0
    1358             :       !
    1359             : 
    1360             :       !(1) prepare list of G vectors
    1361           0 :       call timestart("prep list of Gvec")
    1362           0 :       pointer = 0
    1363           0 :       ic = 0
    1364           0 :       DO ig1 = 1, lapw%nv(jsp)
    1365           0 :          DO igptm = 1, hybrid%ngptm(iq)
    1366           0 :             iigptm = hybrid%pgptm(igptm, iq)
    1367           0 :             g(1) = lapw%k1(ig1, jsp) + hybrid%gptm(1, iigptm) - g_t(1)
    1368           0 :             g(2) = lapw%k2(ig1, jsp) + hybrid%gptm(2, iigptm) - g_t(2)
    1369           0 :             g(3) = lapw%k3(ig1, jsp) + hybrid%gptm(3, iigptm) - g_t(3)
    1370           0 :             IF (pointer(g(1), g(2), g(3)) == 0) THEN
    1371           0 :                ic = ic + 1
    1372           0 :                gpt0(:, ic) = g
    1373           0 :                pointer(g(1), g(2), g(3)) = ic
    1374             :             END IF
    1375             :          END DO
    1376             :       END DO
    1377           0 :       ngpt0 = ic
    1378           0 :       call timestop("prep list of Gvec")
    1379             : 
    1380             :       !(2) calculate convolution
    1381           0 :       call timestart("calc convolution")
    1382           0 :       ALLOCATE (z0(bandoi:bandof, ngpt0), stat=ok)
    1383           0 :       IF (ok /= 0) STOP 'wavefproducts_inv5: error allocation z0'
    1384           0 :       z0 = 0
    1385           0 :       call timestart("step function")
    1386           0 :       DO ig2 = 1, lapw_nkqpt%nv(jsp)
    1387           0 :          rarr1 = z_kqpt%data_r(ig2, bandoi:bandof)
    1388           0 :          DO ig = 1, ngpt0
    1389           0 :             g(1) = gpt0(1, ig) - lapw_nkqpt%k1(ig2, jsp)
    1390           0 :             g(2) = gpt0(2, ig) - lapw_nkqpt%k2(ig2, jsp)
    1391           0 :             g(3) = gpt0(3, ig) - lapw_nkqpt%k3(ig2, jsp)
    1392           0 :             rdum = hybdat%stepfunc_r(g(1), g(2), g(3))/svol
    1393           0 :             DO n2 = bandoi, bandof
    1394           0 :                z0(n2, ig) = z0(n2, ig) + rarr1(n2)*rdum
    1395             :             END DO
    1396             :          END DO
    1397             :       END DO
    1398           0 :       call timestop("step function")
    1399             : 
    1400           0 :       call timestart("hybrid gptm")
    1401           0 :       ic = nbasm_mt
    1402           0 :       DO igptm = 1, hybrid%ngptm(iq)
    1403           0 :          rarr2 = 0
    1404           0 :          ic = ic + 1
    1405           0 :          iigptm = hybrid%pgptm(igptm, iq)
    1406             : 
    1407           0 :          DO ig1 = 1, lapw%nv(jsp)
    1408           0 :             g(1) = lapw%k1(ig1, jsp) + hybrid%gptm(1, iigptm) - g_t(1)
    1409           0 :             g(2) = lapw%k2(ig1, jsp) + hybrid%gptm(2, iigptm) - g_t(2)
    1410           0 :             g(3) = lapw%k3(ig1, jsp) + hybrid%gptm(3, iigptm) - g_t(3)
    1411             : 
    1412           0 :             ig2 = pointer(g(1), g(2), g(3))
    1413             : 
    1414           0 :             IF (ig2 == 0) THEN
    1415           0 :                STOP 'wavefproducts_inv5: pointer undefined'
    1416             :             END IF
    1417             : 
    1418           0 :             DO n1 = 1, bandf - bandi + 1
    1419           0 :                rdum1 = z_nk%data_r(ig1, n1)
    1420           0 :                DO n2 = bandoi, bandof
    1421           0 :                   rarr2(n2, n1) = rarr2(n2, n1) + rdum1*z0(n2, ig2)
    1422             :                END DO
    1423             :             END DO
    1424             : 
    1425             :          END DO
    1426           0 :          cprod(ic, :, :) = rarr2(:, :)
    1427             :       END DO
    1428           0 :       call timestop("hybrid gptm")
    1429           0 :       call timestop("calc convolution")
    1430             : 
    1431           0 :       WRITE (2005, *) 'Point B'
    1432           0 :       DO n2 = 1, 1
    1433           0 :          DO n1 = 1, 2
    1434           0 :             DO ic = 1, 20
    1435           0 :                WRITE (2010, '(3i7,f15.8)') ic, n1, n2, cprod(ic, n1, n2)
    1436             :             END DO
    1437             :          END DO
    1438             :       END DO
    1439             : 
    1440           0 :       DEALLOCATE (z0, pointer, gpt0)
    1441           0 :       CALL timestop("wavefproducts_inv5 IR")
    1442             : 
    1443             :       ! lmstart = lm start index for each l-quantum number and atom type (for cmt-coefficients)
    1444           0 :       DO itype = 1, atoms%ntype
    1445           0 :          DO l = 0, atoms%lmax(itype)
    1446           0 :             lmstart(l, itype) = sum((/(hybrid%nindx(ll, itype)*(2*ll + 1), ll=0, l - 1)/))
    1447             :          END DO
    1448             :       END DO
    1449             : 
    1450             :       ! read in cmt coefficient at k-point nk
    1451           0 :       ALLOCATE (ccmt_nk(dimension%neigd, hybrid%maxlmindx, atoms%nat), ccmt(dimension%neigd, hybrid%maxlmindx, atoms%nat), stat=ok)
    1452           0 :       IF (ok /= 0) STOP 'wavefproducts_inv5: error allocation ccmt_nk/ccmt'
    1453             : 
    1454           0 :       call read_cmt(ccmt_nk, nk)
    1455             :       !read in cmt coefficients at k+q point
    1456           0 :       call read_cmt(ccmt, nkqpt)
    1457             : 
    1458           0 :       iatom = 0
    1459           0 :       DO itype = 1, atoms%ntype
    1460           0 :          DO ieq = 1, atoms%neq(itype)
    1461           0 :             iatom = iatom + 1
    1462             : 
    1463           0 :             cexp(iatom) = exp((-img)*tpi_const*dotprod(kpts%bkf(:, iq) + kpts%bkf(:, nk), atoms%taual(:, iatom)))
    1464             : 
    1465           0 :             cexp_nk(iatom) = exp((-img)*tpi_const*dotprod(kpts%bkf(:, nk), atoms%taual(:, iatom)))
    1466             :          END DO
    1467             :       END DO
    1468             : 
    1469             :       rfac = 1./sr2
    1470             :       cfac = -img/sr2
    1471             :       iatom = 0
    1472           0 :       DO itype = 1, atoms%ntype
    1473           0 :          DO ieq = 1, atoms%neq(itype)
    1474           0 :             iatom = iatom + 1
    1475             :             ! determine the number of the inverse symmetric atom belonging to iatom
    1476           0 :             IF (sym%invsatnr(iatom) == 0) THEN
    1477             :                iiatom = iatom
    1478             :             ELSE
    1479           0 :                iiatom = sym%invsatnr(iatom)
    1480             :             END IF
    1481             :             ! the cmt coefficients at iatom and iiatom are made real in one step
    1482           0 :             IF (iiatom < iatom) CYCLE
    1483           0 :             lm1 = 0
    1484           0 :             DO l = 0, atoms%lmax(itype)
    1485           0 :                DO m = -l, l
    1486           0 :                   rdum = (-1)**(l + m)
    1487           0 :                   DO p = 1, hybrid%nindx(l, itype)
    1488           0 :                      lm1 = lm1 + 1
    1489             :                      ! lm index at l,-m
    1490           0 :                      lm2 = lm1 - 2*m*hybrid%nindx(l, itype)
    1491             : 
    1492           0 :                      IF (iatom == iiatom) THEN
    1493           0 :                         IF (m < 0) THEN
    1494           0 :                            cmt(:, lm1, iatom) = (ccmt(:, lm1, iatom) + rdum*ccmt(:, lm2, iiatom))*cexp(iatom)*rfac
    1495             : 
    1496           0 :                            cmt_nk(:, lm1, iatom) = (ccmt_nk(:, lm1, iatom) + rdum*ccmt_nk(:, lm2, iiatom))*cexp_nk(iatom)*rfac
    1497           0 :                         ELSE IF (m > 0) THEN
    1498             : 
    1499           0 :                            cmt(:, lm1, iatom) = (ccmt(:, lm1, iatom) - rdum*ccmt(:, lm2, iiatom))*cexp(iatom)*cfac
    1500             : 
    1501           0 :                            cmt_nk(:, lm1, iatom) = (ccmt_nk(:, lm1, iatom) - rdum*ccmt_nk(:, lm2, iiatom))*cexp_nk(iatom)*cfac
    1502             :                         ELSE
    1503           0 :                            IF (mod(l, 2) == 0) THEN
    1504           0 :                               cmt(:, lm1, iatom) = ccmt(:, lm1, iatom)*cexp(iatom)
    1505           0 :                               cmt_nk(:, lm1, iatom) = ccmt_nk(:, lm1, iatom)*cexp_nk(iatom)
    1506             :                            ELSE
    1507           0 :                               cmt(:, lm1, iatom) = ccmt(:, lm1, iatom)*(-img)*cexp(iatom)
    1508           0 :                               cmt_nk(:, lm1, iatom) = ccmt_nk(:, lm1, iatom)*(-img)*cexp_nk(iatom)
    1509             :                            END IF
    1510             :                         END IF
    1511             :                      ELSE
    1512           0 :                         cdum = rdum*cexp(iatom)*cexp(iiatom)
    1513           0 :                         cmt(:, lm1, iatom) = (ccmt(:, lm1, iatom) + cdum*ccmt(:, lm2, iiatom))*rfac
    1514             : 
    1515           0 :                         cmt(:, lm1, iiatom) = (ccmt(:, lm1, iatom) - cdum*ccmt(:, lm2, iiatom))*cfac
    1516             : 
    1517           0 :                         cdum = rdum*cexp_nk(iatom)*cexp_nk(iiatom)
    1518           0 :                         cmt_nk(:, lm1, iatom) = (ccmt_nk(:, lm1, iatom) + cdum*ccmt_nk(:, lm2, iiatom))*rfac
    1519             : 
    1520           0 :                         cmt_nk(:, lm1, iiatom) = (ccmt_nk(:, lm1, iatom) - cdum*ccmt_nk(:, lm2, iiatom))*cfac
    1521             :                      END IF
    1522             : 
    1523             :                   END DO
    1524             :                END DO
    1525             :             END DO
    1526             : 
    1527             :          END DO
    1528             :       END DO
    1529           0 :       DEALLOCATE (ccmt_nk, ccmt)
    1530             : 
    1531           0 :       lm_0 = 0
    1532           0 :       lm_00 = 0
    1533           0 :       iatom1 = 0
    1534           0 :       iiatom = 0
    1535             : 
    1536           0 :       DO itype = 1, atoms%ntype
    1537           0 :          ioffset = sum((/((2*ll + 1)*hybrid%nindxm1(ll, itype), ll=0, hybrid%lcutm1(itype))/))
    1538           0 :          lm_0 = lm_00
    1539           0 :          DO ieq = 1, atoms%neq(itype)
    1540           0 :             iatom1 = iatom1 + 1
    1541           0 :             IF (sym%invsatnr(iatom1) == 0) THEN
    1542             :                iatom2 = iatom1
    1543             :             ELSE
    1544           0 :                iatom2 = sym%invsatnr(iatom1)
    1545             :             END IF
    1546           0 :             IF (iatom1 > iatom2) CYCLE
    1547             : 
    1548           0 :             IF (iatom1 /= iatom2) THEN
    1549           0 :                call timestart("iatom1 neq iatom2")
    1550             :                ! loop over l of mixed basis
    1551           0 :                DO l = 0, hybrid%lcutm1(itype)
    1552             :                   ! loop over basis functions products, which belong to l
    1553           0 :                   DO n = 1, hybdat%nindxp1(l, itype)
    1554             : 
    1555             :                      ! determine l1,p1 and l2,p2 for the basis functions, which can generate l
    1556           0 :                      l1 = hybdat%prod(n, l, itype)%l1
    1557           0 :                      l2 = hybdat%prod(n, l, itype)%l2
    1558           0 :                      p1 = hybdat%prod(n, l, itype)%n1
    1559           0 :                      p2 = hybdat%prod(n, l, itype)%n2
    1560             : 
    1561             :                      ! condition for Gaunt coefficients
    1562           0 :                      IF (mod(l + l1 + l2, 2) /= 0) CYCLE
    1563             : 
    1564           0 :                      offdiag = l1 /= l2 .or. p1 /= p2 ! offdiag=true means that b1*b2 and b2*b1 are different combinations
    1565             :                      !(leading to the same basis-function product)
    1566             : 
    1567           0 :                      lm1_0 = lmstart(l1, itype) ! start at correct lm index of cmt-coefficients
    1568           0 :                      lm2_0 = lmstart(l2, itype) ! (corresponding to l1 and l2)
    1569             : 
    1570           0 :                      lm = lm_0
    1571           0 :                      lp1 = lm1_0 + p1
    1572           0 :                      lp2 = lm2_0 + p2
    1573             : 
    1574             :                      ! sum over different m of mixed basis functions with qn l
    1575           0 :                      DO m = -l, l
    1576           0 :                         rarr3 = 0.0
    1577             : 
    1578             :                         ! go to lm index for m1=-l1
    1579           0 :                         lmp1 = lm1_0 + p1
    1580             : 
    1581           0 :                         DO m1 = -l1, l1
    1582             :                            ! Gaunt condition -m1+m2-m=0
    1583           0 :                            m2 = m1 + m
    1584           0 :                            IF (abs(m2) <= l2) THEN
    1585           0 :                               lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1586             :                               ! precalculated Gaunt coefficient
    1587           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, m)
    1588           0 :                               IF (rdum /= 0) THEN
    1589           0 :                                  DO iband = bandi, bandf
    1590           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)
    1591           0 :                                     rdum2 = rdum*cmt_nk(iband, lmp1, iatom2)
    1592             :                                     ! loop over occupied bands
    1593           0 :                                     DO ibando = bandoi, bandof
    1594             : 
    1595             :                                        rarr3(1, ibando, iband) = rarr3(1, ibando, iband)&
    1596           0 :                 &                    + rdum1*cmt(ibando, lmp2, iatom1) + rdum2*cmt(ibando, lmp2, iatom2)
    1597             : 
    1598             :                                        rarr3(2, ibando, iband) = rarr3(2, ibando, iband)&
    1599           0 :                 &                    + rdum1*cmt(ibando, lmp2, iatom2) - rdum2*cmt(ibando, lmp2, iatom1)
    1600             : 
    1601             :                                     END DO  !ibando
    1602             :                                  END DO  !iband
    1603             :                               END IF  ! rdum
    1604             :                            END IF  ! abs(m2) .le. l2
    1605             : 
    1606           0 :                            m2 = m1 - m ! switch role of b1 and b2
    1607           0 :                            IF (abs(m2) <= l2 .and. offdiag) THEN
    1608           0 :                               lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1609           0 :                               rdum = hybdat%gauntarr(2, l1, l2, l, m1, m) ! precalculated Gaunt coefficient
    1610           0 :                               IF (rdum /= 0) THEN
    1611           0 :                                  DO iband = bandi, bandf
    1612           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)
    1613           0 :                                     rdum2 = rdum*cmt_nk(iband, lmp2, iatom2)
    1614             :                                     ! loop over occupied bands
    1615           0 :                                     DO ibando = bandoi, bandof
    1616             :                                        rarr3(1, ibando, iband) = rarr3(1, ibando, iband)&
    1617           0 :                 &                    + rdum1*cmt(ibando, lmp1, iatom1) + rdum2*cmt(ibando, lmp1, iatom2)
    1618             : 
    1619             :                                        rarr3(2, ibando, iband) = rarr3(2, ibando, iband)&
    1620           0 :                 &                    + rdum1*cmt(ibando, lmp1, iatom2) - rdum2*cmt(ibando, lmp1, iatom1)
    1621             :                                     END DO  !ibando
    1622             :                                  END DO  !iband
    1623             :                               END IF  ! rdum .ne. 0
    1624             :                            END IF  ! abs(m2) .le. l2 .and. offdiag
    1625             : 
    1626             :                            ! go to lmp start index for next m1-quantum number
    1627           0 :                            lmp1 = lmp1 + hybrid%nindx(l1, itype)
    1628             : 
    1629             :                         END DO  !m1
    1630             : 
    1631           0 :                         ishift = -2*m*hybrid%nindxm1(l, itype)
    1632             : 
    1633             :                         ! go to lm mixed basis startindx for l and m
    1634           0 :                         lm1 = lm + (iatom1 - 1 - iiatom)*ioffset
    1635           0 :                         lm2 = lm + (iatom2 - 1 - iiatom)*ioffset + ishift
    1636             : 
    1637           0 :                         rdum = tpi_const*dotprod(kpts%bkf(:, iq), atoms%taual(:, iatom1))
    1638           0 :                         rfac1 = sin(rdum)/sr2
    1639           0 :                         rfac2 = cos(rdum)/sr2
    1640           0 :                         DO iband = bandi, bandf
    1641           0 :                            DO ibando = bandoi, bandof
    1642           0 :                               rdum1 = rarr3(1, ibando, iband)
    1643           0 :                               rdum2 = rarr3(2, ibando, iband)
    1644             : !                       sin1  = rdum1*rfac1
    1645             : !                       cos1  = rdum1*rfac2
    1646             : !                       sin2  = rdum2*rfac1
    1647             : !                       cos2  = rdum2*rfac2
    1648           0 :                               add1 = rdum1*rfac2 + rdum2*rfac1
    1649           0 :                               add2 = rdum2*rfac2 - rdum1*rfac1
    1650           0 :                               DO i = 1, hybrid%nindxm1(l, itype)
    1651           0 :                                  j = lm1 + i
    1652           0 :                                  cprod(j, ibando, iband) = cprod(j, ibando, iband) + hybdat%prodm(i, n, l, itype)*add1!( cos1 + sin2 )
    1653           0 :                                  j = lm2 + i
    1654           0 :                                  cprod(j, ibando, iband) = cprod(j, ibando, iband) + hybdat%prodm(i, n, l, itype)*add2!( cos2 - sin1 )
    1655             : 
    1656             :                               END DO  !i -> loop over mixed basis functions
    1657             :                            END DO  !ibando
    1658             :                         END DO  !iband
    1659             : 
    1660             :                         ! go to lm start index for next m-quantum number
    1661           0 :                         lm = lm + hybrid%nindxm1(l, itype)
    1662             : 
    1663             :                      END DO  !m
    1664             : 
    1665             :                   END DO !n
    1666           0 :                   lm_0 = lm_0 + hybrid%nindxm1(l, itype)*(2*l + 1) ! go to the lm start index of the next l-quantum number
    1667           0 :                   IF (lm /= lm_0) STOP 'wavefproducts_inv5: counting of lm-index incorrect (bug?)'
    1668             :                END DO !l
    1669           0 :                call timestop("iatom1 neq iatom2")
    1670             :             ELSE !case: iatom1==iatom2
    1671           0 :                call timestart("iatom1 eq iatom2")
    1672             : 
    1673             :                ! loop over l of mixed basis
    1674           0 :                monepl = -1
    1675           0 :                DO l = 0, hybrid%lcutm1(itype)
    1676           0 :                   monepl = -monepl
    1677             :                   ! loop over basis functions products, which belong to l
    1678           0 :                   DO n = 1, hybdat%nindxp1(l, itype)
    1679             : 
    1680             :                      ! determine l1,p1 and l2,p2 for the basis functions, which can generate l
    1681           0 :                      l1 = hybdat%prod(n, l, itype)%l1
    1682           0 :                      l2 = hybdat%prod(n, l, itype)%l2
    1683           0 :                      p1 = hybdat%prod(n, l, itype)%n1
    1684           0 :                      p2 = hybdat%prod(n, l, itype)%n2
    1685             : 
    1686             :                      ! condition for Gaunt coefficients
    1687           0 :                      IF (mod(l + l1 + l2, 2) /= 0) CYCLE
    1688             : 
    1689           0 :                      offdiag = l1 /= l2 .or. p1 /= p2 ! offdiag=true means that b1*b2 and b2*b1 are different combinations
    1690             :                      !(leading to the same basis-function product)
    1691             : 
    1692           0 :                      lm1_0 = lmstart(l1, itype) ! start at correct lm index of cmt-coefficients
    1693           0 :                      lm2_0 = lmstart(l2, itype) ! (corresponding to l1 and l2)
    1694             : 
    1695           0 :                      lm = lm_0
    1696           0 :                      lp1 = lm1_0 + p1
    1697           0 :                      lp2 = lm2_0 + p2
    1698             : 
    1699             :                      ! calculate (-1)**l1 and (-1)**l2 (monepl = minus one power l)
    1700           0 :                      monepl1 = (-1)**l1
    1701           0 :                      monepl2 = (-1)**l2
    1702             : 
    1703             :                      ! sum over different m of mixed basis functions with qn l
    1704             : 
    1705             :                      !
    1706             :                      !case m<0
    1707             :                      !
    1708             : 
    1709           0 :                      monepm = -monepl
    1710           0 :                      DO m = -l, -1
    1711           0 :                         monepm = -monepm
    1712           0 :                         moneplm = monepl*monepm
    1713             : 
    1714             :                         ! calculate the contributions which are identical for m>0 and m <0
    1715           0 :                         rarr2 = 0.0
    1716           0 :                         IF (abs(m) <= l2) THEN
    1717           0 :                            lmp1 = lp1 + l1*hybrid%nindx(l1, itype)
    1718           0 :                            IF (mod(l1, 2) == 0) THEN
    1719           0 :                               lmp2 = lp2 + (m + l2)*hybrid%nindx(l2, itype)
    1720             :                            ELSE
    1721           0 :                               lmp2 = lp2 + (-m + l2)*hybrid%nindx(l2, itype)
    1722             :                            END IF
    1723             : 
    1724           0 :                            rdum = hybdat%gauntarr(1, l1, l2, l, 0, m)
    1725           0 :                            IF (rdum /= 0) THEN
    1726           0 :                               DO iband = bandi, bandf
    1727           0 :                                  rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)
    1728           0 :                                  IF (mod(l1, 2) /= 0) rdum1 = moneplm*rdum1
    1729           0 :                                  DO ibando = bandoi, bandof
    1730           0 :                                     rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    1731             :                                  END DO  ! ibando
    1732             :                               END DO  ! iband
    1733             :                            END IF  ! rdum .ne. 0
    1734             : 
    1735           0 :                            IF (offdiag) THEN
    1736           0 :                               rdum = hybdat%gauntarr(1, l2, l1, l, -m, m)
    1737           0 :                               IF (rdum /= 0) THEN
    1738           0 :                                  DO iband = bandi, bandf
    1739           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)
    1740           0 :                                     IF (mod(l1, 2) == 0) rdum1 = moneplm*rdum1
    1741           0 :                                     DO ibando = bandoi, bandof
    1742           0 :                                        rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp1, iatom1)
    1743             :                                     END DO  ! ibando
    1744             :                                  END DO  ! iband
    1745             :                               END IF  ! rdum .ne. 0
    1746             :                            END IF  ! offdiag
    1747             : 
    1748             :                         END IF  ! abs(m) .le. l2
    1749             : 
    1750           0 :                         IF (abs(m) <= l1) THEN
    1751           0 :                            IF (mod(l2, 2) == 0) THEN
    1752           0 :                               lmp3 = lp1 + (m + l1)*hybrid%nindx(l1, itype)
    1753             :                            ELSE
    1754           0 :                               lmp3 = lp1 + (-m + l1)*hybrid%nindx(l1, itype)
    1755             :                            END IF
    1756           0 :                            lmp2 = lp2 + l2*hybrid%nindx(l2, itype)
    1757             : 
    1758           0 :                            rdum = hybdat%gauntarr(1, l1, l2, l, -m, m)
    1759           0 :                            IF (rdum /= 0) THEN
    1760           0 :                               DO iband = bandi, bandf
    1761           0 :                                  rdum1 = rdum*cmt_nk(iband, lmp3, iatom1)
    1762           0 :                                  IF (mod(l2, 2) == 0) rdum1 = moneplm*rdum1
    1763           0 :                                  DO ibando = bandoi, bandof
    1764           0 :                                     rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    1765             :                                  END DO  ! ibando
    1766             :                               END DO  ! iband
    1767             :                            END IF  ! rdum .ne. 0
    1768             : 
    1769           0 :                            IF (offdiag) THEN
    1770           0 :                               rdum = hybdat%gauntarr(1, l2, l1, l, 0, m)
    1771           0 :                               IF (rdum /= 0) THEN
    1772           0 :                                  DO iband = bandi, bandf
    1773           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)
    1774           0 :                                     IF (mod(l2, 2) /= 0) rdum1 = moneplm*rdum1
    1775           0 :                                     DO ibando = bandoi, bandof
    1776           0 :                                        rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
    1777             :                                     END DO  ! ibando
    1778             :                                  END DO  ! iband
    1779             :                               END IF  ! rdum .ne. 0
    1780             :                            END IF  ! offdiag
    1781             : 
    1782             :                         END IF  ! abs(m) .le. l2
    1783             : 
    1784             :                         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1785             : 
    1786             :                         !go to lm index for m1=-l1
    1787           0 :                         lmp1 = lp1
    1788           0 :                         monepm1 = -monepl1
    1789           0 :                         DO m1 = -l1, l1
    1790           0 :                            monepm1 = -monepm1
    1791           0 :                            IF (m1 == 0) THEN
    1792           0 :                               lmp1 = lmp1 + hybrid%nindx(l1, itype)
    1793           0 :                               CYCLE
    1794             :                            END IF
    1795             :                            ! (-1)**(l1+m1)
    1796           0 :                            monepl1m1 = monepl1*monepm1
    1797           0 :                            m2 = m1 + m
    1798           0 :                            IF (abs(m2) <= l2 .and. m2 /= 0) THEN
    1799           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, m)
    1800           0 :                               IF (rdum /= 0) THEN
    1801           0 :                                  IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    1802           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1803             :                                  ELSE
    1804           0 :                                     lmp2 = lp2 + (-m2 + l2)*hybrid%nindx(l2, itype)
    1805           0 :                                     fac = 1/2.*moneplm*monepl1m1*(sign(1, m2) - sign(1, m1))
    1806             :                                  END IF
    1807           0 :                                  rdum = rdum/sr2
    1808           0 :                                  DO iband = bandi, bandf
    1809           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)!rdum*cmt_nk(iband,lmp1,iatom1)/sr2
    1810           0 :                                     IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
    1811           0 :                                     DO ibando = bandoi, bandof
    1812           0 :                                        rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    1813             :                                     END DO  ! ibando
    1814             :                                  END DO  ! iband
    1815             :                               END IF  ! rdum .ne. 0
    1816             : 
    1817           0 :                               IF (offdiag) THEN
    1818           0 :                                  rdum = hybdat%gauntarr(1, l2, l1, l, m2, -m)
    1819           0 :                                  IF (rdum /= 0) THEN
    1820           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1821           0 :                                     IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    1822             :                                        lmp3 = lmp1
    1823             :                                     ELSE
    1824           0 :                                        lmp3 = lmp1 - 2*m1*hybrid%nindx(l1, itype)
    1825           0 :                                        fac = 1/2.*monepl1m1*(sign(1, m1) - sign(1, m2))
    1826             :                                     END IF
    1827           0 :                                     rdum = moneplm*rdum/sr2
    1828           0 :                                     DO iband = bandi, bandf
    1829           0 :                                        rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)!moneplm*rdum*cmt_nk(iband,lmp2,iatom1)/sr2
    1830           0 :                                        IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
    1831           0 :                                        DO ibando = bandoi, bandof
    1832           0 :                                           rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
    1833             :                                        END DO  ! ibando
    1834             :                                     END DO  ! iband
    1835             :                                  END IF  ! rdum .ne. 0
    1836             :                               END IF  ! offdiag
    1837             : 
    1838             :                            END IF  ! abs(m2) .le. l2 .and. m2 .ne. 0
    1839             : 
    1840           0 :                            m2 = m1 - m
    1841           0 :                            IF (abs(m2) <= l2 .and. m2 /= 0) THEN
    1842             : 
    1843           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, -m)
    1844           0 :                               IF (rdum /= 0) THEN
    1845             : 
    1846           0 :                                  IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    1847           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1848             :                                  ELSE
    1849           0 :                                     lmp2 = lp2 + (-m2 + l2)*hybrid%nindx(l2, itype)
    1850           0 :                                     fac = 1/2.*moneplm*monepl1m1*(sign(1, m2) - sign(1, m1))
    1851             :                                  END IF
    1852           0 :                                  rdum = moneplm*rdum/sr2
    1853           0 :                                  DO iband = bandi, bandf
    1854           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)!moneplm*rdum*cmt_nk(iband,lmp1,iatom1)/sr2
    1855           0 :                                     IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
    1856           0 :                                     DO ibando = bandoi, bandof
    1857           0 :                                        rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    1858             :                                     END DO  ! ibando
    1859             :                                  END DO  ! iband
    1860             : 
    1861             :                               END IF  ! rdum .ne. 0
    1862             : 
    1863           0 :                               IF (offdiag) THEN
    1864           0 :                                  rdum = hybdat%gauntarr(1, l2, l1, l, m2, m)
    1865           0 :                                  IF (rdum /= 0) THEN
    1866           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1867           0 :                                     IF (sign(1, m1) + sign(1, m2) /= 0) THEN
    1868             :                                        lmp3 = lmp1
    1869             :                                     ELSE
    1870           0 :                                        lmp3 = lmp1 - 2*m1*hybrid%nindx(l1, itype)
    1871           0 :                                        fac = 1/2.*monepl1m1*(sign(1, m1) - sign(1, m2))
    1872             :                                     END IF
    1873           0 :                                     rdum = rdum/sr2
    1874           0 :                                     DO iband = bandi, bandf
    1875           0 :                                        rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)!rdum*cmt_nk(iband,lmp2,iatom1)/sr2
    1876           0 :                                        IF (sign(1, m1) + sign(1, m2) == 0) rdum1 = fac*rdum1
    1877           0 :                                        DO ibando = bandoi, bandof
    1878           0 :                                           rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
    1879             :                                        END DO  ! ibando
    1880             :                                     END DO  ! iband
    1881             :                                  END IF  ! rdum .ne. 0
    1882             :                               END IF  ! offdiag
    1883             : 
    1884             :                            END IF  ! abs(m2) .le. l2 .and. m1 .ne. 0
    1885             : 
    1886             :                            !go to lmp start index for next m1-quantum number
    1887           0 :                            lmp1 = lmp1 + hybrid%nindx(l1, itype)
    1888             :                         END DO  ! m1
    1889             : 
    1890             :                         ! go to lm mixed basis startindx for l and m
    1891           0 :                         lm1 = lm + (iatom1 - 1 - iiatom)*ioffset
    1892           0 :                         DO iband = bandi, bandf
    1893           0 :                            DO ibando = bandoi, bandof
    1894           0 :                               rdum = rarr2(ibando, iband)
    1895           0 :                               DO i = 1, hybrid%nindxm1(l, itype)
    1896           0 :                                  j = lm1 + i
    1897           0 :                                  cprod(j, ibando, iband) = cprod(j, ibando, iband) + hybdat%prodm(i, n, l, itype)*rdum
    1898             :                               END DO  !i -> loop over mixed basis functions
    1899             :                            END DO  !ibando
    1900             :                         END DO  !iband
    1901             : 
    1902             :                         ! go to lm start index for next m-quantum number
    1903           0 :                         lm = lm + hybrid%nindxm1(l, itype)
    1904             : 
    1905             :                      END DO  ! m=-l,-1
    1906             : 
    1907             :                      !
    1908             :                      !case m=0
    1909             :                      !
    1910             : 
    1911             :                      m = 0
    1912           0 :                      rarr2 = 0.0
    1913           0 :                      lmp1 = lp1
    1914             : 
    1915           0 :                      monepm1 = -monepl1
    1916             : 
    1917           0 :                      DO m1 = -l1, l1
    1918           0 :                         m2 = m1
    1919           0 :                         monepm1 = -monepm1
    1920           0 :                         IF (abs(m2) <= l2) THEN
    1921             : 
    1922           0 :                            IF (mod(l, 2) == 0) THEN
    1923           0 :                               lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1924             :                               !lmp3 and lmp4 are variables, which avoid an if clause in the loop
    1925           0 :                               lmp3 = lmp2
    1926           0 :                               lmp4 = lmp1
    1927             :                            ELSE
    1928           0 :                               lmp2 = lp2 + (-m2 + l2)*hybrid%nindx(l2, itype)
    1929             :                               !lmp3 and lmp3 are variables, which avoid an if clause in the loop
    1930           0 :                               lmp3 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    1931           0 :                               lmp4 = lmp1 - 2*m1*hybrid%nindx(l1, itype)
    1932             : 
    1933           0 :                               fac1 = monepl1*monepm1 ! (-1)**(l1+m1)
    1934           0 :                               fac2 = monepl2*monepm1 ! (-1)**(l2+m1)
    1935             :                            END IF
    1936             : 
    1937             :                            !precalculated Gaunt coefficient
    1938           0 :                            IF (mod(l, 2) == 0) THEN
    1939           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, m)
    1940             :                            ELSE
    1941           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, m)*fac1
    1942             :                            END IF
    1943           0 :                            IF (rdum /= 0) THEN
    1944           0 :                               DO iband = bandi, bandf
    1945           0 :                                  rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)
    1946           0 :                                  DO ibando = bandoi, bandof
    1947           0 :                                     rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    1948             :                                  END DO  ! ibando
    1949             :                               END DO  ! iband
    1950             :                            END IF  ! rdum.ne.0
    1951             : 
    1952             :                            !change role of b1 and b2
    1953           0 :                            IF (offdiag) THEN
    1954           0 :                               IF (mod(l, 2) == 0) THEN
    1955           0 :                                  rdum = hybdat%gauntarr(2, l1, l2, l, m1, m)
    1956             :                               ELSE
    1957           0 :                                  rdum = hybdat%gauntarr(2, l1, l2, l, m1, m)*fac2
    1958             :                               END IF
    1959           0 :                               IF (rdum /= 0) THEN
    1960           0 :                                  DO iband = bandi, bandf
    1961           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp3, iatom1)
    1962           0 :                                     DO ibando = bandoi, bandof
    1963           0 :                                        rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp4, iatom1)
    1964             :                                     END DO  ! ibando
    1965             :                                  END DO  ! iband
    1966             :                               END IF  ! rdum.ne.0
    1967             :                            END IF  ! offdiag
    1968             : 
    1969             :                         END IF  ! abs(m2).le.l2
    1970             : 
    1971             :                         ! go to lmp start index for next m1-quantum number
    1972           0 :                         lmp1 = lmp1 + hybrid%nindx(l1, itype)
    1973             :                      END DO  !m1
    1974             : 
    1975             :                      ! go to lm mixed basis startindx for l and m
    1976           0 :                      lm1 = lm + (iatom1 - 1 - iiatom)*ioffset
    1977           0 :                      DO iband = bandi, bandf
    1978           0 :                         DO ibando = bandoi, bandof
    1979           0 :                            rdum = rarr2(ibando, iband)
    1980           0 :                            DO i = 1, hybrid%nindxm1(l, itype)
    1981           0 :                               j = lm1 + i
    1982           0 :                               cprod(j, ibando, iband) = cprod(j, ibando, iband) + hybdat%prodm(i, n, l, itype)*rdum
    1983             :                            END DO  !i -> loop over mixed basis functions
    1984             :                         END DO  !ibando
    1985             :                      END DO  !iband
    1986             : 
    1987             :                      ! go to lm start index for next m-quantum number
    1988           0 :                      lm = lm + hybrid%nindxm1(l, itype)
    1989             : 
    1990             :                      !
    1991             :                      ! case: m>0
    1992             :                      !
    1993             : 
    1994           0 :                      rarr2 = 0.0
    1995             :                      monepm = 1
    1996           0 :                      DO m = 1, l
    1997           0 :                         monepm = -monepm
    1998           0 :                         moneplm = monepl*monepm
    1999             : 
    2000             :                         ! calculate the contributions which are identical for m>0 and m <0
    2001           0 :                         rarr2 = 0.0
    2002           0 :                         IF (abs(m) <= l2) THEN
    2003           0 :                            lmp1 = lp1 + l1*hybrid%nindx(l1, itype)
    2004           0 :                            IF (mod(l1, 2) == 0) THEN
    2005           0 :                               lmp2 = lp2 + (m + l2)*hybrid%nindx(l2, itype)
    2006             :                            ELSE
    2007           0 :                               lmp2 = lp2 + (-m + l2)*hybrid%nindx(l2, itype)
    2008             :                            END IF
    2009             : 
    2010           0 :                            rdum = hybdat%gauntarr(1, l1, l2, l, 0, m)
    2011           0 :                            IF (rdum /= 0) THEN
    2012           0 :                               DO iband = bandi, bandf
    2013           0 :                                  rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)
    2014           0 :                                  IF (mod(l1, 2) /= 0) rdum1 = moneplm*rdum1
    2015           0 :                                  DO ibando = bandoi, bandof
    2016           0 :                                     rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    2017             :                                  END DO  ! ibando
    2018             :                               END DO  ! iband
    2019             :                            END IF  ! rdum .ne. 0
    2020             : 
    2021           0 :                            IF (offdiag) THEN
    2022           0 :                               rdum = hybdat%gauntarr(1, l2, l1, l, -m, m)
    2023           0 :                               IF (rdum /= 0) THEN
    2024           0 :                                  DO iband = bandi, bandf
    2025           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)
    2026           0 :                                     IF (mod(l1, 2) == 0) rdum1 = moneplm*rdum1
    2027           0 :                                     DO ibando = bandoi, bandof
    2028           0 :                                        rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp1, iatom1)
    2029             :                                     END DO  ! ibando
    2030             :                                  END DO  ! iband
    2031             :                               END IF  ! rdum .ne. 0
    2032             :                            END IF  ! offdiag
    2033             : 
    2034             :                         END IF  ! abs(m) .le. l2
    2035             : 
    2036           0 :                         IF (abs(m) <= l1) THEN
    2037           0 :                            IF (mod(l2, 2) == 0) THEN
    2038           0 :                               lmp3 = lp1 + (m + l1)*hybrid%nindx(l1, itype)
    2039             :                            ELSE
    2040           0 :                               lmp3 = lp1 + (-m + l1)*hybrid%nindx(l1, itype)
    2041             :                            END IF
    2042           0 :                            lmp2 = lp2 + l2*hybrid%nindx(l2, itype)
    2043             : 
    2044           0 :                            rdum = hybdat%gauntarr(1, l1, l2, l, -m, m)
    2045           0 :                            IF (rdum /= 0) THEN
    2046           0 :                               DO iband = bandi, bandf
    2047           0 :                                  rdum1 = rdum*cmt_nk(iband, lmp3, iatom1)
    2048           0 :                                  IF (mod(l2, 2) == 0) rdum1 = moneplm*rdum1
    2049           0 :                                  DO ibando = bandoi, bandof
    2050           0 :                                     rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    2051             :                                  END DO  ! ibando
    2052             :                               END DO  ! iband
    2053             :                            END IF  ! rdum .ne. 0
    2054             : 
    2055           0 :                            IF (offdiag) THEN
    2056           0 :                               rdum = hybdat%gauntarr(1, l2, l1, l, 0, m)
    2057           0 :                               IF (rdum /= 0) THEN
    2058           0 :                                  DO iband = bandi, bandf
    2059           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)
    2060           0 :                                     IF (mod(l2, 2) /= 0) rdum1 = moneplm*rdum1
    2061           0 :                                     DO ibando = bandoi, bandof
    2062           0 :                                        rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
    2063             :                                     END DO  ! ibando
    2064             :                                  END DO  ! iband
    2065             :                               END IF  ! rdum .ne. 0
    2066             :                            END IF  ! offdiag
    2067             : 
    2068             :                         END IF  ! abs(m) .le. l2
    2069             : 
    2070             :                         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2071             : 
    2072             :                         !go to lm index for m1=-l1
    2073             :                         lmp1 = lp1
    2074             :                         monepm1 = -monepl1
    2075           0 :                         DO m1 = -l1, l1
    2076           0 :                            monepm1 = -monepm1
    2077           0 :                            IF (m1 == 0) THEN
    2078           0 :                               lmp1 = lmp1 + hybrid%nindx(l1, itype)
    2079           0 :                               CYCLE
    2080             :                            END IF
    2081           0 :                            m2 = m1 + m
    2082             :                            ! (-1)**(l1+m1)
    2083           0 :                            monepl1m1 = monepl1*monepm1
    2084           0 :                            IF (abs(m2) <= l2 .and. m2 /= 0) THEN
    2085           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, m)
    2086           0 :                               IF (rdum /= 0) THEN
    2087             : 
    2088           0 :                                  IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    2089           0 :                                     lmp2 = lp2 + (-m2 + l2)*hybrid%nindx(l2, itype)
    2090             :                                  ELSE
    2091           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    2092           0 :                                     fac = -moneplm*monepl1m1*(sign(1, m2) - sign(1, m1))/2
    2093             :                                  END IF
    2094             : 
    2095           0 :                                  rdum = -moneplm*monepl1m1*rdum/sr2
    2096           0 :                                  DO iband = bandi, bandf
    2097           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)!-moneplm*monepl1m1*rdum*cmt_nk(iband,lmp1,iatom1)/sr2
    2098           0 :                                     IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
    2099           0 :                                     DO ibando = bandoi, bandof
    2100           0 :                                        rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    2101             :                                     END DO  ! ibando
    2102             :                                  END DO  ! iband
    2103             : 
    2104             :                               END IF  ! rdum .ne. 0
    2105             : 
    2106           0 :                               IF (offdiag) THEN
    2107           0 :                                  rdum = hybdat%gauntarr(2, l1, l2, l, m1, -m)
    2108           0 :                                  IF (rdum /= 0) THEN
    2109           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    2110           0 :                                     IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    2111           0 :                                        lmp3 = lmp1 - 2*m1*hybrid%nindx(l1, itype)
    2112             :                                     ELSE
    2113           0 :                                        lmp3 = lmp1
    2114           0 :                                        fac = 1/2.*monepl1m1*(sign(1, m2) - sign(1, m1))
    2115             :                                     END IF
    2116           0 :                                     rdum = monepl1m1*moneplm*rdum/sr2
    2117           0 :                                     DO iband = bandi, bandf
    2118           0 :                                        rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)!monepl1m1*moneplm*rdum*cmt_nk(iband,lmp2,iatom1)/sr2
    2119           0 :                                        IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
    2120           0 :                                        DO ibando = bandoi, bandof
    2121           0 :                                           rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
    2122             :                                        END DO  ! ibando
    2123             :                                     END DO  ! iband
    2124             : 
    2125             :                                  END IF  ! rdum
    2126             : 
    2127             :                               END IF  ! offdiag
    2128             :                            END IF  ! abs(m2) .le. l2 .and. m2 .ne. 0
    2129             : 
    2130           0 :                            m2 = m1 - m
    2131           0 :                            IF (abs(m2) <= l2 .and. m2 /= 0) THEN
    2132             : 
    2133           0 :                               rdum = hybdat%gauntarr(1, l1, l2, l, m1, -m)
    2134           0 :                               IF (rdum /= 0) THEN
    2135             : 
    2136           0 :                                  IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    2137           0 :                                     lmp2 = lp2 + (-m2 + l2)*hybrid%nindx(l2, itype)
    2138             :                                  ELSE
    2139           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    2140           0 :                                     fac = 1/2.*moneplm*monepl1m1*(sign(1, m1) - sign(1, m2))
    2141             :                                  END IF
    2142           0 :                                  rdum = monepl1m1*rdum/sr2
    2143           0 :                                  DO iband = bandi, bandf
    2144           0 :                                     rdum1 = rdum*cmt_nk(iband, lmp1, iatom1)!monepl1m1*rdum*cmt_nk(iband,lmp1,iatom1)/sr2
    2145           0 :                                     IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = rdum1*fac
    2146           0 :                                     DO ibando = bandoi, bandof
    2147           0 :                                        rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp2, iatom1)
    2148             :                                     END DO  ! ibando
    2149             :                                  END DO  ! iband
    2150             :                               END IF  ! rdum .ne. 0
    2151             : 
    2152           0 :                               IF (offdiag) THEN
    2153           0 :                                  rdum = hybdat%gauntarr(2, l1, l2, l, m1, m)
    2154           0 :                                  IF (rdum /= 0) THEN
    2155           0 :                                     lmp2 = lp2 + (m2 + l2)*hybrid%nindx(l2, itype)
    2156           0 :                                     IF (sign(1, m2) + sign(1, m1) /= 0) THEN
    2157           0 :                                        lmp3 = lmp1 - 2*m1*hybrid%nindx(l1, itype)
    2158             :                                     ELSE
    2159           0 :                                        lmp3 = lmp1
    2160           0 :                                        fac = -monepl1m1*(sign(1, m1) - sign(1, m2))/2
    2161             :                                     END IF
    2162           0 :                                     rdum = -monepl1m1*rdum/sr2
    2163           0 :                                     DO iband = bandi, bandf
    2164           0 :                                        rdum1 = rdum*cmt_nk(iband, lmp2, iatom1)!-monepl1m1*rdum*cmt_nk(iband,lmp2,iatom1)/sr2
    2165           0 :                                        IF (sign(1, m2) + sign(1, m1) == 0) rdum1 = fac*rdum1
    2166           0 :                                        DO ibando = bandoi, bandof
    2167           0 :                                           rarr2(ibando, iband) = rarr2(ibando, iband) + rdum1*cmt(ibando, lmp3, iatom1)
    2168             :                                        END DO  ! ibando
    2169             :                                     END DO  ! iband
    2170             : 
    2171             :                                  END IF  ! rdum .ne. 0
    2172             :                               END IF  ! offdiag
    2173             : 
    2174             :                            END IF  !  abs(m2) .le. l2 .and. m2 .ne. 0
    2175             : 
    2176             :                            !go to lmp start index for next m1-quantum number
    2177           0 :                            lmp1 = lmp1 + hybrid%nindx(l1, itype)
    2178             :                         END DO  ! m1
    2179             : 
    2180             :                         ! multiply rarr2 by (-1)**(l+m+1)
    2181           0 :                         rarr2(:, :) = (-1)*moneplm*rarr2(:, :)
    2182             : 
    2183             :                         ! go to lm mixed basis startindx for l and m
    2184           0 :                         lm1 = lm + (iatom1 - 1 - iiatom)*ioffset
    2185             : 
    2186           0 :                         DO iband = bandi, bandf
    2187           0 :                            DO ibando = bandoi, bandof
    2188           0 :                               rdum = rarr2(ibando, iband)
    2189           0 :                               DO i = 1, hybrid%nindxm1(l, itype)
    2190           0 :                                  j = lm1 + i
    2191           0 :                                  cprod(j, ibando, iband) = cprod(j, ibando, iband) + hybdat%prodm(i, n, l, itype)*rdum
    2192             :                               END DO  !i -> loop over mixed basis functions
    2193             :                            END DO  !ibando
    2194             :                         END DO  !iband
    2195             : 
    2196             :                         ! go to lm start index for next m-quantum number
    2197           0 :                         lm = lm + hybrid%nindxm1(l, itype)
    2198             : 
    2199             :                      END DO  ! m=1,l
    2200             : 
    2201             :                   END DO !n
    2202           0 :                   lm_0 = lm_0 + hybrid%nindxm1(l, itype)*(2*l + 1) ! go to the m start index of the next l-quantum number
    2203           0 :                   IF (lm /= lm_0) STOP 'wavefproducts_inv5: counting of lm-index incorrect (bug?)'
    2204             :                END DO !l
    2205             : 
    2206           0 :                call timestop("iatom1 eq iatom2")
    2207             :             END IF  ! iatom1 .ne. iatom2
    2208             : 
    2209           0 :             lm_0 = lm_00
    2210             :          END DO !ieq
    2211           0 :          iiatom = iiatom + atoms%neq(itype)
    2212           0 :          lm_00 = lm_00 + atoms%neq(itype)*ioffset
    2213             :       END DO  !itype
    2214           0 :       CALL timestop("wavefproducts_inv5")
    2215             : 
    2216           0 :    END SUBROUTINE wavefproducts_inv5
    2217             : 
    2218           0 :    SUBROUTINE wavefproducts_noinv5(&
    2219             :   &                      bandi, bandf, bandoi, bandof,&
    2220             :   &                      nk, iq, dimension, input, jsp,&
    2221             :   &                      cell, atoms, hybrid,&
    2222             :   &                      hybdat,&
    2223             :   &                      kpts,&
    2224             :   &                      mnobd,&
    2225             :   &                      lapw, sym,&
    2226             :   &                      nbasm_mt,&
    2227             :   &                      noco,&
    2228           0 :   &                      nkqpt, cprod)
    2229             : 
    2230             :       USE m_constants
    2231             :       USE m_util, ONLY: modulo1
    2232             :       USE m_olap, ONLY: gptnorm
    2233             :       USE m_trafo
    2234             :       USE m_wrapper
    2235             :       USE m_types
    2236             :       USE m_io_hybrid
    2237             :       IMPLICIT NONE
    2238             :       TYPE(t_dimension), INTENT(IN)   :: dimension
    2239             :       TYPE(t_input), INTENT(IN)       :: input
    2240             :       TYPE(t_noco), INTENT(IN)        :: noco
    2241             :       TYPE(t_sym), INTENT(IN)         :: sym
    2242             :       TYPE(t_cell), INTENT(IN)        :: cell
    2243             :       TYPE(t_kpts), INTENT(IN)        :: kpts
    2244             :       TYPE(t_atoms), INTENT(IN)       :: atoms
    2245             :       TYPE(t_lapw), INTENT(IN)        :: lapw
    2246             :       TYPE(t_hybrid), INTENT(IN)      :: hybrid
    2247             :       TYPE(t_hybdat), INTENT(INOUT)   :: hybdat
    2248             : 
    2249             : !     - scalars -
    2250             :       INTEGER, INTENT(IN)      ::  bandi, bandf, bandoi, bandof
    2251             :       INTEGER, INTENT(IN)      ::  nk, iq, jsp
    2252             :       INTEGER, INTENT(IN)      :: mnobd, nbasm_mt
    2253             :       INTEGER, INTENT(OUT)     ::  nkqpt
    2254             : 
    2255             : !     - arrays -
    2256             : 
    2257             :       COMPLEX, INTENT(OUT)    ::  cprod(hybrid%maxbasm1, bandoi:bandof, bandf - bandi + 1)
    2258             : 
    2259             : !     - local scalars -
    2260             :       INTEGER                 ::  ic, l, n, l1, l2, n1, n2, lm_0, lm1_0, lm2_0, lm, lm1, lm2, m1, m2, i, j, ll
    2261             :       INTEGER                 ::  itype, ieq, ikpt, ikpt1, ikpt2, igpt, igptp, igpt1, igpt2, iband, iband1, iband2
    2262             :       INTEGER                 ::  k, ic1, ioffset, ibando, ig1, ig2, ig
    2263             :       INTEGER                 ::  igptm, iigptm
    2264             :       INTEGER                 ::  q, idum
    2265             :       INTEGER                 :: nbasm_ir, ngpt0
    2266             :       INTEGER                 ::  nbasmmt, nbasfcn
    2267             :       INTEGER                 ::  ok, m
    2268             : 
    2269             :       REAL                    ::  rdum, svol, s2
    2270             : 
    2271             :       COMPLEX                 ::  cdum, cdum0, cdum1
    2272             :       COMPLEX                 ::  cexp
    2273             :       COMPLEX, PARAMETER       ::  img = (0.0, 1.0)
    2274             : 
    2275             :       LOGICAL                 ::  offdiag
    2276           0 :       TYPE(t_lapw)            ::    lapw_nkqpt
    2277             : 
    2278             : !      - local arrays -
    2279             :       INTEGER                 ::  g(3), g_t(3)
    2280           0 :       INTEGER                 ::  lmstart(0:atoms%lmaxd, atoms%ntype)
    2281           0 :       INTEGER, ALLOCATABLE    ::  gpt0(:, :)
    2282           0 :       INTEGER, ALLOCATABLE    ::  pointer(:, :, :)
    2283             : 
    2284             :       REAL                    ::  bkpt(3)
    2285             :       REAL                    ::  kqpt(3), kqpthlp(3)
    2286             : 
    2287           0 :       COMPLEX                 ::  carr1(bandoi:bandof)
    2288           0 :       COMPLEX                 ::  carr2(bandoi:bandof, bandf - bandi + 1)
    2289           0 :       TYPE(t_mat)             ::  z_nk, z_kqpt
    2290           0 :       COMPLEX                 ::  cmt(dimension%neigd, hybrid%maxlmindx, atoms%nat)
    2291           0 :       COMPLEX                 ::  cmt_nk(dimension%neigd, hybrid%maxlmindx, atoms%nat)
    2292           0 :       COMPLEX, ALLOCATABLE     ::  z0(:, :)
    2293             : 
    2294           0 :       call timestart("wavefproducts_noinv5")
    2295           0 :       call timestart("wavefproducts_noinv5 IR")
    2296           0 :       cprod = 0
    2297           0 :       svol = sqrt(cell%omtil)
    2298           0 :       s2 = sqrt(2.0)
    2299             : 
    2300           0 :       nbasm_ir = maxval(hybrid%ngptm)
    2301             : 
    2302             :       !
    2303             :       ! compute k+q point for given q point in EIBZ(k)
    2304             :       !
    2305           0 :       kqpthlp = kpts%bkf(:, nk) + kpts%bkf(:, iq)
    2306             :       ! k+q can lie outside the first BZ, transfer
    2307             :       ! it back into the 1. BZ
    2308           0 :       kqpt = modulo1(kqpthlp, kpts%nkpt3)
    2309           0 :       g_t(:) = nint(kqpt - kqpthlp)
    2310             :       ! determine number of kqpt
    2311           0 :       nkqpt = 0
    2312           0 :       DO ikpt = 1, kpts%nkptf
    2313           0 :          IF (maxval(abs(kqpt - kpts%bkf(:, ikpt))) <= 1E-06) THEN
    2314           0 :             nkqpt = ikpt
    2315           0 :             EXIT
    2316             :          END IF
    2317             :       END DO
    2318           0 :       IF (nkqpt == 0) STOP 'wavefproducts: k-point not found'
    2319             : 
    2320             :       !
    2321             :       ! compute G's fulfilling |bk(:,nkqpt) + G| <= rkmax
    2322             :       !
    2323           0 :       CALL lapw_nkqpt%init(input, noco, kpts, atoms, sym, nkqpt, cell, sym%zrfs)
    2324           0 :       nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
    2325           0 :       call z_nk%alloc(.false., nbasfcn, dimension%neigd)
    2326           0 :       nbasfcn = MERGE(lapw_nkqpt%nv(1) + lapw_nkqpt%nv(2) + 2*atoms%nlotot, lapw_nkqpt%nv(1) + atoms%nlotot, noco%l_noco)
    2327           0 :       call z_kqpt%alloc(.false., nbasfcn, dimension%neigd)
    2328             : 
    2329             :       ! read in z at k-point nk and nkqpt
    2330           0 :       call timestart("read_z")
    2331           0 :       call read_z(z_nk, nk)
    2332           0 :       call read_z(z_kqpt, nkqpt)
    2333           0 :       call timestop("read_z")
    2334             : 
    2335             :       g(1) = maxval(abs(lapw%k1(:lapw%nv(jsp), jsp))) &
    2336             :      &     + maxval(abs(lapw_nkqpt%k1(:lapw_nkqpt%nv(jsp), jsp)))&
    2337           0 :      &     + maxval(abs(hybrid%gptm(1, hybrid%pgptm(:hybrid%ngptm(iq), iq)))) + 1
    2338             :       g(2) = maxval(abs(lapw%k2(:lapw%nv(jsp), jsp)))&
    2339             :      &     + maxval(abs(lapw_nkqpt%k2(:lapw_nkqpt%nv(jsp), jsp)))&
    2340           0 :      &     + maxval(abs(hybrid%gptm(2, hybrid%pgptm(:hybrid%ngptm(iq), iq)))) + 1
    2341             :       g(3) = maxval(abs(lapw%k3(:lapw%nv(jsp), jsp)))&
    2342             :      &     + maxval(abs(lapw_nkqpt%k3(:lapw_nkqpt%nv(jsp), jsp)))&
    2343           0 :      &     + maxval(abs(hybrid%gptm(3, hybrid%pgptm(:hybrid%ngptm(iq), iq)))) + 1
    2344             : 
    2345           0 :       ALLOCATE (pointer(-g(1):g(1), -g(2):g(2), -g(3):g(3)), stat=ok)
    2346           0 :       IF (ok /= 0) STOP 'wavefproducts_noinv2: error allocation pointer'
    2347           0 :       ALLOCATE (gpt0(3, size(pointer)), stat=ok)
    2348           0 :       IF (ok /= 0) STOP 'wavefproducts_noinv2: error allocation gpt0'
    2349             : 
    2350           0 :       if (.not. allocated(hybdat%stepfunc_c)) then
    2351           0 :          call timestart("setup stepfunc")
    2352           0 :          ALLOCATE (hybdat%stepfunc_c(-g(1):g(1), -g(2):g(2), -g(3):g(3)), stat=ok)
    2353           0 :          IF (ok /= 0) then
    2354           0 :             call juDFT_error('wavefproducts_noinv2: error allocation stepfunc')
    2355             :          endif
    2356           0 :          DO i = -g(1), g(1)
    2357           0 :             DO j = -g(2), g(2)
    2358           0 :                DO k = -g(3), g(3)
    2359           0 :                   hybdat%stepfunc_c(i, j, k) = stepfunction(cell, atoms, (/i, j, k/))
    2360             :                END DO
    2361             :             END DO
    2362             :          END DO
    2363           0 :          call timestop("setup stepfunc")
    2364             :       endif
    2365             : 
    2366             :       !
    2367             :       ! convolute phi(n,k) with the step function and store in cpw0
    2368             :       !
    2369             : 
    2370             :       !(1) prepare list of G vectors
    2371           0 :       call timestart("prep list of Gvec")
    2372           0 :       pointer = 0
    2373           0 :       ic = 0
    2374           0 :       DO ig1 = 1, lapw%nv(jsp)
    2375           0 :          DO igptm = 1, hybrid%ngptm(iq)
    2376           0 :             iigptm = hybrid%pgptm(igptm, iq)
    2377           0 :             g(1) = lapw%k1(ig1, jsp) + hybrid%gptm(1, iigptm) - g_t(1)
    2378           0 :             g(2) = lapw%k2(ig1, jsp) + hybrid%gptm(2, iigptm) - g_t(2)
    2379           0 :             g(3) = lapw%k3(ig1, jsp) + hybrid%gptm(3, iigptm) - g_t(3)
    2380           0 :             IF (pointer(g(1), g(2), g(3)) == 0) THEN
    2381           0 :                ic = ic + 1
    2382           0 :                gpt0(:, ic) = g
    2383           0 :                pointer(g(1), g(2), g(3)) = ic
    2384             :             END IF
    2385             :          END DO
    2386             :       END DO
    2387           0 :       ngpt0 = ic
    2388           0 :       call timestop("prep list of Gvec")
    2389             : 
    2390             :       !(2) calculate convolution
    2391           0 :       call timestart("calc convolution")
    2392           0 :       call timestart("step function")
    2393           0 :       ALLOCATE (z0(bandoi:bandof, ngpt0))
    2394           0 :       z0 = 0
    2395           0 :       DO ig2 = 1, lapw_nkqpt%nv(jsp)
    2396           0 :          carr1 = z_kqpt%data_c(ig2, bandoi:bandof)
    2397           0 :          DO ig = 1, ngpt0
    2398           0 :             g(1) = gpt0(1, ig) - lapw_nkqpt%k1(ig2, jsp)
    2399           0 :             g(2) = gpt0(2, ig) - lapw_nkqpt%k2(ig2, jsp)
    2400           0 :             g(3) = gpt0(3, ig) - lapw_nkqpt%k3(ig2, jsp)
    2401           0 :             cdum = hybdat%stepfunc_c(g(1), g(2), g(3))/svol
    2402           0 :             DO n2 = bandoi, bandof
    2403           0 :                z0(n2, ig) = z0(n2, ig) + carr1(n2)*cdum
    2404             :             END DO
    2405             :          END DO
    2406             :       END DO
    2407           0 :       call timestop("step function")
    2408             : 
    2409           0 :       call timestart("hybrid gptm")
    2410           0 :       ic = nbasm_mt
    2411           0 :       DO igptm = 1, hybrid%ngptm(iq)
    2412           0 :          carr2 = 0
    2413           0 :          ic = ic + 1
    2414           0 :          iigptm = hybrid%pgptm(igptm, iq)
    2415             : 
    2416           0 :          DO ig1 = 1, lapw%nv(jsp)
    2417           0 :             g(1) = lapw%k1(ig1, jsp) + hybrid%gptm(1, iigptm) - g_t(1)
    2418           0 :             g(2) = lapw%k2(ig1, jsp) + hybrid%gptm(2, iigptm) - g_t(2)
    2419           0 :             g(3) = lapw%k3(ig1, jsp) + hybrid%gptm(3, iigptm) - g_t(3)
    2420             : 
    2421           0 :             ig2 = pointer(g(1), g(2), g(3))
    2422             : 
    2423           0 :             IF (ig2 == 0) THEN
    2424           0 :                STOP 'wavefproducts_noinv2: pointer undefined'
    2425             :             END IF
    2426             : 
    2427           0 :             DO n1 = 1, bandf - bandi + 1
    2428           0 :                cdum1 = conjg(z_nk%data_c(ig1, n1))
    2429           0 :                DO n2 = bandoi, bandof
    2430           0 :                   carr2(n2, n1) = carr2(n2, n1) + cdum1*z0(n2, ig2)
    2431             :                END DO
    2432             :             END DO
    2433             : 
    2434             :          END DO
    2435           0 :          cprod(ic, :, :) = carr2(:, :)
    2436             :       END DO
    2437           0 :       call timestop("hybrid gptm")
    2438           0 :       DEALLOCATE (z0, pointer, gpt0)
    2439           0 :       call timestop("calc convolution")
    2440             : 
    2441           0 :       call timestop("wavefproducts_noinv5 IR")
    2442             : 
    2443             : !       RETURN
    2444             : 
    2445             :       !
    2446             :       ! MT contribution
    2447             :       !
    2448             : 
    2449             :       ! lmstart = lm start index for each l-quantum number and atom type (for cmt-coefficients)
    2450           0 :       DO itype = 1, atoms%ntype
    2451           0 :          DO l = 0, atoms%lmax(itype)
    2452           0 :             lmstart(l, itype) = sum((/(hybrid%nindx(ll, itype)*(2*ll + 1), ll=0, l - 1)/))
    2453             :          END DO
    2454             :       END DO
    2455             : 
    2456             :       ! read in cmt coefficients from direct access file cmt
    2457           0 :       call read_cmt(cmt_nk(:, :, :), nk)
    2458           0 :       call read_cmt(cmt(:, :, :), nkqpt)
    2459             : 
    2460           0 :       lm_0 = 0
    2461           0 :       ic = 0
    2462             : 
    2463           0 :       DO itype = 1, atoms%ntype
    2464           0 :          DO ieq = 1, atoms%neq(itype)
    2465           0 :             ic = ic + 1
    2466           0 :             ic1 = 0
    2467             : 
    2468           0 :             cexp = exp(-img*tpi_const*dot_product(kpts%bkf(:, iq), atoms%taual(:, ic)))
    2469             : 
    2470           0 :             DO l = 0, hybrid%lcutm1(itype)
    2471           0 :                DO n = 1, hybdat%nindxp1(l, itype) ! loop over basis-function products
    2472             : 
    2473           0 :                   l1 = hybdat%prod(n, l, itype)%l1 !
    2474           0 :                   l2 = hybdat%prod(n, l, itype)%l2 ! current basis-function product
    2475           0 :                   n1 = hybdat%prod(n, l, itype)%n1 ! = bas(:,n1,l1,itype)*bas(:,n2,l2,itype) = b1*b2
    2476           0 :                   n2 = hybdat%prod(n, l, itype)%n2 !
    2477             : 
    2478           0 :                   IF (mod(l1 + l2 + l, 2) /= 0) cycle
    2479             : 
    2480           0 :                   offdiag = l1 /= l2 .or. n1 /= n2 ! offdiag=true means that b1*b2 and b2*b1 are different combinations
    2481             :                   !(leading to the same basis-function product)
    2482             : 
    2483           0 :                   lm1_0 = lmstart(l1, itype) ! start at correct lm index of cmt-coefficients
    2484           0 :                   lm2_0 = lmstart(l2, itype) ! (corresponding to l1 and l2)
    2485             : 
    2486           0 :                   lm = lm_0
    2487           0 :                   DO m = -l, l
    2488             : 
    2489           0 :                      carr2 = 0.0
    2490             : 
    2491           0 :                      lm1 = lm1_0 + n1 ! go to lm index for m1=-l1
    2492           0 :                      DO m1 = -l1, l1
    2493           0 :                         m2 = m1 + m ! Gaunt condition -m1+m2-m=0
    2494           0 :                         IF (abs(m2) <= l2) THEN
    2495           0 :                            lm2 = lm2_0 + n2 + (m2 + l2)*hybrid%nindx(l2, itype)
    2496           0 :                            rdum = hybdat%gauntarr(1, l1, l2, l, m1, m) ! precalculated Gaunt coefficient
    2497           0 :                            IF (rdum /= 0) THEN
    2498           0 :                               DO iband = bandi, bandf
    2499           0 :                                  cdum = rdum*conjg(cmt_nk(iband, lm1, ic)) !nk
    2500           0 :                                  DO iband1 = bandoi, bandof
    2501           0 :                                     carr2(iband1, iband) = carr2(iband1, iband) + cdum*cmt(iband1, lm2, ic) !ikpt
    2502             : 
    2503             :                                  END DO
    2504             :                               END DO
    2505             :                            END IF
    2506             :                         END IF
    2507             : 
    2508           0 :                         m2 = m1 - m ! switch role of b1 and b2
    2509           0 :                         IF (abs(m2) <= l2 .and. offdiag) THEN
    2510           0 :                            lm2 = lm2_0 + n2 + (m2 + l2)*hybrid%nindx(l2, itype)
    2511           0 :                            rdum = hybdat%gauntarr(2, l1, l2, l, m1, m) ! precalculated Gaunt coefficient
    2512           0 :                            IF (rdum /= 0) THEN
    2513           0 :                               DO iband = bandi, bandf
    2514           0 :                                  cdum = rdum*conjg(cmt_nk(iband, lm2, ic)) !nk
    2515           0 :                                  DO iband1 = bandoi, bandof
    2516           0 :                                     carr2(iband1, iband) = carr2(iband1, iband) + cdum*cmt(iband1, lm1, ic)
    2517             :                                  END DO
    2518             :                               END DO
    2519             :                            END IF
    2520             :                         END IF
    2521             : 
    2522           0 :                         lm1 = lm1 + hybrid%nindx(l1, itype) ! go to lm start index for next m1-quantum number
    2523             : 
    2524             :                      END DO  !m1
    2525             : 
    2526           0 :                      DO iband = bandi, bandf
    2527           0 :                         DO iband1 = bandoi, bandof
    2528           0 :                            cdum = carr2(iband1, iband)*cexp
    2529           0 :                            DO i = 1, hybrid%nindxm1(l, itype)
    2530           0 :                               j = lm + i
    2531           0 :                               cprod(j, iband1, iband) = cprod(j, iband1, iband) + hybdat%prodm(i, n, l, itype)*cdum
    2532             :                            END DO
    2533             : 
    2534             :                         END DO
    2535             :                      END DO
    2536             : 
    2537           0 :                      lm = lm + hybrid%nindxm1(l, itype) ! go to lm start index for next m-quantum number
    2538             : 
    2539             :                   END DO
    2540             : 
    2541             :                END DO
    2542           0 :                lm_0 = lm_0 + hybrid%nindxm1(l, itype)*(2*l + 1) ! go to the lm start index of the next l-quantum number
    2543           0 :                IF (lm /= lm_0) STOP 'wavefproducts_noinv2: counting of lm-index incorrect (bug?)'
    2544             :             END DO
    2545             :          END DO
    2546             :       END DO
    2547             : 
    2548           0 :       call timestop("wavefproducts_noinv5")
    2549             : 
    2550           0 :    END SUBROUTINE wavefproducts_noinv5
    2551             : 
    2552             :    !private subroutine
    2553           0 :    FUNCTION stepfunction(cell, atoms, g)
    2554             :       USE m_types
    2555             :       USE m_constants
    2556             :       USE m_olap
    2557             :       IMPLICIT NONE
    2558             : 
    2559             :       TYPE(t_cell), INTENT(IN)    :: cell
    2560             :       TYPE(t_atoms), INTENT(IN)   :: atoms
    2561             : 
    2562             :       INTEGER, INTENT(IN) :: g(3)
    2563             :       COMPLEX             :: stepfunction  !Is real in inversion case
    2564             :       REAL                :: gnorm, gnorm3, r, fgr
    2565             :       INTEGER             :: itype, ieq, icent
    2566             : 
    2567           0 :       gnorm = gptnorm(g, cell%bmat)
    2568           0 :       gnorm3 = gnorm**3
    2569           0 :       IF (gnorm == 0) THEN
    2570           0 :          stepfunction = 1
    2571           0 :          DO itype = 1, atoms%ntype
    2572           0 :             stepfunction = stepfunction - atoms%neq(itype)*atoms%volmts(itype)/cell%omtil
    2573             :          END DO
    2574             :       ELSE
    2575           0 :          stepfunction = 0
    2576           0 :          icent = 0
    2577           0 :          DO itype = 1, atoms%ntype
    2578           0 :             r = gnorm*atoms%rmt(itype)
    2579           0 :             fgr = fpi_const*(sin(r) - r*cos(r))/gnorm3/cell%omtil
    2580           0 :             DO ieq = 1, atoms%neq(itype)
    2581           0 :                icent = icent + 1
    2582           0 :                stepfunction = stepfunction - fgr*exp(-cmplx(0., tpi_const*dot_product(atoms%taual(:, icent), g)))
    2583             :             ENDDO
    2584             :          ENDDO
    2585             :       ENDIF
    2586             : 
    2587           0 :    END FUNCTION stepfunction
    2588             : 
    2589             : END MODULE m_wavefproducts

Generated by: LCOV version 1.13