LCOV - code coverage report
Current view: top level - hybrid - coulombmatrix.F90 (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 539 545 98.9 %
Date: 2024-03-28 04:22:06 Functions: 8 8 100.0 %

          Line data    Source code
       1             : !
       2             : !     Calculates the Coulomb matrix
       3             : !
       4             : !     v      =  < M    | v | M    >
       5             : !      k,IJ        k,I        k,J
       6             : !
       7             : !     with the mixed-basis functions M (indices I and J).
       8             : !
       9             : !     Note that
      10             : !                 *
      11             : !     v      =  v     .
      12             : !      k,JI      k,IJ
      13             : !
      14             : !     In the code: coulomb(IJ,k) = v     where only the upper triangle (I<=J) is stored.
      15             : !                                   k,IJ
      16             : !
      17             : !     The Coulomb matrix v(IJ,k) diverges at the Gamma-point. Here, we apply the decomposition
      18             : !
      19             : !              (0)        (1)   *        2-l              (0)*   (0)    (1)*        m  (1)
      20             : !     v     = v    + SUM v   * Y  (k) / k        with    v    = v   ,  v      = (-1)  v
      21             : !      k,IJ    IJ     lm  IJ    lm                        JI     IJ     JI,lm          IJ,l,-m
      22             : !
      23             : !     where a = atom index, R  = position vector, T  = Wigner-Seitz radius (scalar).
      24             : !                            a                     0
      25             : !                                    (0)
      26             : !     In the code: coulomb(IJ,1)  = v    where only the upper triangle (I<=J) is stored,
      27             : !                                    IJ
      28             : !                                    (1)
      29             : !                  coulfac(IJ,lm) = v                                    IJ,lm
      30             : !
      31             : !     For the PW contribution we have to construct plane waves within the MT spheres with the help
      32             : !     of spherical Bessel functions. The value lexp (LEXP in gwinp) is the corresponding cutoff.
      33             : !
      34             : MODULE m_coulombmatrix
      35             : #ifdef CPP_MPI
      36             :    use mpi
      37             : #endif
      38             :    use m_types
      39             :    USE m_intgrf, ONLY: intgrf, intgrf_init
      40             :    use m_sphbes, only: sphbes
      41             :    use m_glob_tofrom_loc
      42             :    USE m_trafo, ONLY: symmetrize_mpimat, symmetrize, bramat_trafo
      43             :    use m_gamma_double_gpt_loop
      44             : CONTAINS
      45             : 
      46          12 :    SUBROUTINE coulombmatrix(fmpi, fi, mpdata, hybdat, xcpot)
      47             :       use m_work_package
      48             :       use m_structureconstant
      49             :       USE m_types
      50             :       USE m_types_mpimat
      51             :       use m_types_mat
      52             :       USE m_types_hybdat
      53             :       USE m_juDFT
      54             :       USE m_constants
      55             :       use m_util, only: primitivef
      56             :       USE m_hsefunctional, ONLY: change_coulombmatrix
      57             :       USE m_wrapper
      58             :       USE m_io_hybrid
      59             :       use m_ylm
      60             :       use m_calc_l_m_from_lm
      61             :       use m_calc_mpsmat
      62             :       use m_copy_coul
      63             :       use m_apply_inverse_olap
      64             :       IMPLICIT NONE
      65             : 
      66             :       TYPE(t_xcpot_inbuild), INTENT(IN) :: xcpot
      67             :       TYPE(t_mpi), INTENT(IN)           :: fmpi
      68             :       type(t_fleurinput), intent(in)    :: fi
      69             :       TYPE(t_mpdata), intent(in)        :: mpdata
      70             :       TYPE(t_hybdat), INTENT(INOUT)     :: hybdat
      71             : 
      72             :       ! - local scalars -
      73             :       INTEGER                    :: inviop
      74             :       INTEGER                    :: nqnrm, iqnrm, iqnrm1, iqnrm2
      75             :       INTEGER                    :: itype, l, ix, iy, iy0, i, j, lm, l1, l2, m1, m2, ineq, ikpt
      76             :       INTEGER                    :: lm1, lm2, itype1, itype2, ineq1, ineq2, n1, n2, iat2
      77             :       INTEGER                    :: ic, ic1, ic2
      78             :       INTEGER                    :: igpt, igpt1, igpt2, igptp, igptp1, igptp2
      79             :       INTEGER                    :: isym, isym1, isym2, igpt0
      80             :       INTEGER                    :: iatm1, iatm2
      81             :       INTEGER                    :: m, im
      82             :       INTEGER                    :: maxfac, ix_loc, pe, pe_ix
      83             : 
      84             :       LOGICAL                    :: lsym
      85             : 
      86             :       REAL                       :: rdum, rdum1, rdum2
      87             :       REAL                       :: svol, qnorm1, qnorm2
      88             :       REAL                       :: fcoulfac
      89             : 
      90             :       COMPLEX                    :: cdum, cexp, csum
      91             : 
      92             :       ! - local arrays -
      93             :       INTEGER                    :: g(3)
      94          12 :       INTEGER, ALLOCATABLE   :: pqnrm(:, :)
      95          24 :       INTEGER                    :: rrot(3, 3, fi%sym%nsym), invrrot(3, 3, fi%sym%nsym)
      96          24 :       INTEGER, ALLOCATABLE   :: iarr(:), POINTER(:, :, :, :)!,pointer(:,:,:)
      97          24 :       INTEGER, ALLOCATABLE   :: nsym_gpt(:, :), sym_gpt(:, :, :)
      98          24 :       INTEGER                    :: nsym1(fi%kpts%nkpt + 1), sym1(fi%sym%nsym, fi%kpts%nkpt + 1)
      99             : 
     100          12 :       INTEGER, ALLOCATABLE   ::  ngptm1(:)
     101          12 :       INTEGER, ALLOCATABLE   ::  pgptm1(:, :)
     102             : 
     103             :       REAL                       :: q(3), q1(3), q2(3), mtmt_term
     104          24 :       REAL                       :: integrand(fi%atoms%jmtd), primf1(fi%atoms%jmtd), primf2(fi%atoms%jmtd)
     105         152 :       REAL                       :: moment(maxval(mpdata%num_radbasfn), 0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype), &
     106         132 :                                     moment2(maxval(mpdata%num_radbasfn), fi%atoms%ntype)
     107          24 :       REAL, ALLOCATABLE   :: gmat(:, :), qnrm(:)
     108          12 :       REAL, ALLOCATABLE   :: sphbesmoment(:, :, :)
     109          12 :       REAL, ALLOCATABLE   :: sphbes0(:, :, :)
     110          24 :       REAL, ALLOCATABLE   :: olap(:, :, :, :), integral(:, :, :, :)
     111          12 :       REAL, ALLOCATABLE   :: gridf(:, :)
     112          64 :       REAL                       :: facA(0:MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1))
     113          64 :       REAL                       :: facB(0:MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1))
     114          64 :       REAL                       :: facC(-1:MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1))
     115          32 :       REAL    :: sphbes_var(fi%atoms%jmtd, 0:maxval(fi%hybinp%lcutm1))
     116          32 :       REAL    :: sphbesmoment1(fi%atoms%jmtd, 0:maxval(fi%hybinp%lcutm1))
     117             : 
     118          12 :       COMPLEX     :: y((fi%hybinp%lexp + 1)**2), smat
     119         112 :       COMPLEX     :: dwgn(-maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), -maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), 0:maxval(fi%hybinp%lcutm1), fi%sym%nsym)
     120          12 :       COMPLEX, ALLOCATABLE   :: carr2(:, :), carr2a(:, :), carr2b(:, :)
     121          12 :       COMPLEX, ALLOCATABLE   :: structconst(:,:,:,:)
     122             : 
     123             :       INTEGER                    :: ierr
     124             :       INTEGER                    :: iatom, mtmt_idx
     125          12 :       TYPE(t_mat)                :: mat
     126          12 :       type(t_mat), allocatable   :: mtmt_repl(:)
     127          12 :       class(t_mat), allocatable  :: coul(:)
     128             : 
     129          12 :       CALL timestart("Coulomb matrix setup")
     130          12 :       call timestart("prep in coulomb")
     131          12 :       if (fmpi%is_root()) write (*, *) "start of coulomb calculation"
     132             : 
     133          72 :       allocate(structconst((2*fi%hybinp%lexp + 1)**2, fi%atoms%nat, fi%atoms%nat, fi%kpts%nkpt), stat=ierr)
     134          12 :       if(ierr /= 0) call judft_error("can't allocate structconst. error: " // int2str(ierr))
     135             : 
     136          12 :       svol = SQRT(fi%cell%vol)
     137          12 :       fcoulfac = fpi_const/fi%cell%vol
     138          32 :       maxfac = MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1)
     139             : 
     140          12 :       facA(0) = 1                    !
     141          12 :       facB(0) = 1                    ! Define:
     142          36 :       facC(-1:0) = 1                    ! facA(i)    = i!
     143         792 :       DO i = 1, maxfac                       ! facB(i)   = sqrt(i!)
     144         780 :          facA(i) = facA(i - 1)*i            ! facC(i) = (2i+1)!!
     145         780 :          facB(i) = facB(i - 1)*SQRT(i*1.0) !
     146         792 :          facC(i) = facC(i - 1)*(2*i + 1)   !
     147             :       END DO
     148             : 
     149          12 :       CALL intgrf_init(fi%atoms%ntype, fi%atoms%jmtd, fi%atoms%jri, fi%atoms%dx, fi%atoms%rmsh, gridf)
     150             : 
     151             :       !     Calculate the structure constant
     152          12 :       CALL structureconstant(structconst, fi%cell, fi%hybinp, fi%atoms, fi%kpts, fmpi)
     153             : 
     154          12 :       IF (fmpi%irank == 0) WRITE (oUnit, '(//A)') '### subroutine: coulombmatrix ###'
     155             : 
     156             :       !
     157             :       !     Matrix allocation
     158             :       !
     159             : 
     160          12 :       call timestart("coulomb allocation")
     161             : 
     162          12 :       if(fmpi%n_size == 1) then
     163           0 :          allocate(t_mat::coul(fi%kpts%nkpt))
     164             :       else
     165          72 :          allocate(t_mpimat::coul(fi%kpts%nkpt))
     166             :       endif
     167          48 :       do ikpt = 1, fi%kpts%nkpt 
     168          84 :          if(any(ikpt == fmpi%k_list))then 
     169          36 :             call coul(ikpt)%init(.False., hybdat%nbasm(ikpt), hybdat%nbasm(ikpt), fmpi%sub_comm, .false.)
     170             :          else
     171           0 :             call coul(ikpt)%init(.False., 1, 1, fmpi%sub_comm, .false.)
     172             :          endif 
     173             :       enddo
     174          12 :       call timestop("coulomb allocation") 
     175             : 
     176          12 :       IF (fmpi%irank == 0) then
     177             :          write (oUnit,*) "Size of coulomb matrix: " //&
     178          42 :                             float2str(sum([(coul(fmpi%k_list(i))%size_mb(), i=1,size(fmpi%k_list))])) // " MB"
     179             :       endif
     180             : 
     181             :       !     Generate Symmetry:
     182             :       !     Reduce list of g-Points so that only one of each symm-equivalent is calculated
     183             :       ! calculate rotations in reciprocal space
     184         588 :       DO isym = 1, fi%sym%nsym
     185         588 :          IF (isym <= fi%sym%nop) THEN
     186         480 :             inviop = fi%sym%invtab(isym)
     187        6240 :             rrot(:, :, isym) = TRANSPOSE(fi%sym%mrot(:, :, inviop))
     188        4128 :             DO l = 0, maxval(fi%hybinp%lcutm1)
     189             :                dwgn(:, :, l, isym) = TRANSPOSE(fi%hybinp%d_wgn2(-maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), &
     190      228960 :                                                                 -maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), l, isym))
     191             :             END DO
     192             :          ELSE
     193          96 :             inviop = isym - fi%sym%nop
     194        1248 :             rrot(:, :, isym) = -rrot(:, :, inviop)
     195       43776 :             dwgn(:, :, :, isym) = dwgn(:, :, :, inviop)
     196         864 :             DO l = 0, maxval(fi%hybinp%lcutm1)
     197        2976 :                DO m1 = -l, l
     198        9600 :                   DO m2 = -l, -1
     199        6720 :                      cdum = dwgn(m1, m2, l, isym)
     200        6720 :                      dwgn(m1, m2, l, isym) = dwgn(m1, -m2, l, isym)*(-1)**m2
     201        9120 :                      dwgn(m1, -m2, l, isym) = cdum*(-1)**m2
     202             :                   END DO
     203             :                END DO
     204             :             END DO
     205             :          END IF
     206             :       END DO
     207        6744 :       invrrot(:, :, :fi%sym%nop) = rrot(:, :, fi%sym%invtab)
     208          12 :       IF (fi%sym%nsym > fi%sym%nop) THEN
     209        1348 :          invrrot(:, :, fi%sym%nop + 1:) = rrot(:, :, fi%sym%invtab + fi%sym%nop)
     210             :       END IF
     211             : 
     212             :       ! Get symmetry operations that leave bk(:,ikpt) invariant -> sym1
     213          60 :       nsym1 = 0
     214          48 :       DO ikpt = 1, fi%kpts%nkpt
     215          36 :          isym1 = 0
     216        1764 :          DO isym = 1, fi%sym%nsym
     217             :             ! temporary fix until bramat_trafo is correct
     218             :             ! for systems with symmetries including translations
     219        1728 :             IF (isym > fi%sym%nop) THEN
     220         288 :                isym2 = isym - fi%sym%nop
     221             :             ELSE
     222             :                isym2 = isym
     223             :             END IF
     224        6048 :             IF (ANY(abs(fi%sym%tau(:, isym2)) > 1e-12)) CYCLE
     225             : 
     226       34404 :             IF (ALL(ABS(MATMUL(rrot(:, :, isym), fi%kpts%bk(:, ikpt)) - fi%kpts%bk(:, ikpt)) < 1e-12)) THEN
     227         544 :                isym1 = isym1 + 1
     228         544 :                sym1(isym1, ikpt) = isym
     229             :             END IF
     230             :          END DO
     231          48 :          nsym1(ikpt) = isym1
     232             :       END DO
     233             :       ! Define reduced lists of G points -> pgptm1(:,ikpt), ikpt=1,..,nkpt
     234             :       !if(allocated(pgptm1)) deallocate(fi%hybinp%pgptm1)
     235       14512 :       allocate (pgptm1(maxval(mpdata%n_g), fi%kpts%nkptf), source=0) !in mixedbasis
     236        1916 :       allocate (iarr(maxval(mpdata%n_g)), source=0)
     237             :       allocate (POINTER(fi%kpts%nkpt, &
     238             :                         MINVAL(mpdata%g(1, :)) - 1:MAXVAL(mpdata%g(1, :)) + 1, &
     239             :                         MINVAL(mpdata%g(2, :)) - 1:MAXVAL(mpdata%g(2, :)) + 1, &
     240             :                         MINVAL(mpdata%g(3, :)) - 1:MAXVAL(mpdata%g(3, :)) + 1), &
     241       58292 :                 source=0)
     242          36 :       allocate (ngptm1, mold=mpdata%n_g)
     243         108 :       ngptm1 = 0
     244             : 
     245          48 :       DO ikpt = 1, fi%kpts%nkpt
     246        4952 :          DO igpt = 1, mpdata%n_g(ikpt)
     247       19664 :             g = mpdata%g(:, mpdata%gptm_ptr(igpt, ikpt))
     248        4952 :             POINTER(ikpt, g(1), g(2), g(3)) = igpt
     249             :          END DO
     250        5388 :          iarr = 0
     251             :          j = 0
     252        4952 :          DO igpt = mpdata%n_g(ikpt), 1, -1
     253        4952 :             IF (iarr(igpt) == 0) THEN
     254        1432 :                j = j + 1
     255        1432 :                pgptm1(j, ikpt) = igpt
     256       10936 :                DO isym1 = 1, nsym1(ikpt)
     257      123552 :                   g = MATMUL(rrot(:, :, sym1(isym1, ikpt)), mpdata%g(:, mpdata%gptm_ptr(igpt, ikpt)))
     258        9504 :                   i = POINTER(ikpt, g(1), g(2), g(3))
     259        9504 :                   IF (i == 0) call judft_error('coulombmatrix: zero pointer (bug?)')
     260       10936 :                   iarr(i) = 1
     261             :                END DO
     262             :             END IF
     263             :          END DO
     264          48 :          ngptm1(ikpt) = j
     265             :       END DO
     266          12 :       deallocate (iarr)
     267             : 
     268             :       ! Distribute the work as equally as possible over the processes
     269          12 :       call timestop("prep in coulomb")
     270             : 
     271          12 :       call timestart("define gmat")
     272             :       ! Define gmat (symmetric)
     273          48 :       allocate (gmat((fi%hybinp%lexp + 1)**2, (fi%hybinp%lexp + 1)**2))
     274     1005732 :       gmat = 0
     275             :       lm1 = 0
     276         216 :       DO l1 = 0, fi%hybinp%lexp
     277        3684 :          DO m1 = -l1, l1
     278        3468 :             lm1 = lm1 + 1
     279        3468 :             lm2 = 0
     280       41412 :             lp1: DO l2 = 0, l1
     281      547332 :                DO m2 = -l2, l2
     282      506124 :                   lm2 = lm2 + 1
     283      506124 :                   IF (lm2 > lm1) EXIT lp1 ! Don't cross the diagonal!
     284             :                   gmat(lm1, lm2) = facB(l1 + l2 + m2 - m1)*facB(l1 + l2 + m1 - m2)/ &
     285             :                                    (facB(l1 + m1)*facB(l1 - m1)*facB(l2 + m2)*facB(l2 - m2))/ &
     286      502860 :                                    SQRT(1.0*(2*l1 + 1)*(2*l2 + 1)*(2*(l1 + l2) + 1))*(fpi_const)**1.5
     287      540600 :                   gmat(lm2, lm1) = gmat(lm1, lm2)
     288             :                END DO
     289             :             END DO LP1
     290             :          END DO
     291             :       END DO
     292          12 :       call timestop("define gmat")
     293             : 
     294             :       ! Calculate moments of MT functions
     295          12 :       call timestart("calc moments of MT")
     296          32 :       DO itype = 1, fi%atoms%ntype
     297         120 :          DO l = 0, fi%hybinp%lcutm1(itype)
     298         880 :             DO i = 1, mpdata%num_radbasfn(l, itype)
     299             :                ! note that mpdata%radbasfn_mt already contains the factor rgrid
     300             :                moment(i, l, itype) = intgrf(fi%atoms%rmsh(:, itype)**(l + 1)*mpdata%radbasfn_mt(:, i, l, itype), &
     301      618020 :                                             fi%atoms, itype, gridf)
     302             :             END DO
     303             :          END DO
     304         212 :          DO i = 1, mpdata%num_radbasfn(0, itype)
     305             :             moment2(i, itype) = intgrf(fi%atoms%rmsh(:, itype)**3*mpdata%radbasfn_mt(:, i, 0, itype), &
     306      146512 :                                        fi%atoms, itype, gridf)
     307             :          END DO
     308             :       END DO
     309          12 :       call timestop("calc moments of MT")
     310             : 
     311          12 :       call timestart("getnorm")
     312             :       ! Look for different qnorm = |k+G|, definition of qnrm and pqnrm.
     313          12 :       CALL getnorm(fi%kpts, mpdata%g, mpdata%n_g, mpdata%gptm_ptr, qnrm, nqnrm, pqnrm, fi%cell)
     314           0 :       allocate (sphbesmoment(0:fi%hybinp%lexp, fi%atoms%ntype, nqnrm), &
     315           0 :                 olap(maxval(mpdata%num_radbasfn), 0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype, nqnrm), &
     316         460 :                 integral(maxval(mpdata%num_radbasfn), 0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype, nqnrm))
     317       48332 :       sphbes_var = 0
     318       13956 :       sphbesmoment = 0
     319       48332 :       sphbesmoment1 = 0
     320       39092 :       olap = 0
     321       39092 :       integral = 0
     322             : 
     323             :       ! Calculate moments of spherical Bessel functions (for (2) and (3))              (->sphbesmoment)
     324             :       ! Calculate overlap of spherical Bessel functions with basis functions (for (2)) (->olap)
     325             :       ! Calculate overlap of sphbesmoment1(r,l)         with basis functions (for (2)) (->integral)
     326             :       ! We use           sphbes(r,l) = j_l(qr)
     327             :       ! and       sphbesmoment1(r,l) = 1/r**(l-1) * INT(0..r) r'**(l+2) * j_l(qr') dr'
     328             :       !                                + r**(l+2) * INT(r..S) r'**(1-l) * j_l(qr') dr' .
     329             : 
     330          12 :       call timestop("getnorm")
     331             : 
     332          12 :       call bessel_calculation(fi, fmpi, mpdata, nqnrm, gridf, qnrm, sphbesmoment, olap, integral)
     333             : 
     334             :       !
     335             :       !     (1) Case < MT | v | MT >
     336             :       !
     337             : 
     338             : 
     339             :       !       (1a) r,r' in same MT
     340          12 :       call timestart("loop 1")
     341          12 :       ix = 0
     342          12 :       iy = 0
     343          12 :       iy0 = 0
     344             : 
     345             : 
     346         252 :       call mat%alloc(.True., maxval(mpdata%num_radbasfn), maxval(mpdata%num_radbasfn))
     347             : 
     348         536 :       allocate(mtmt_repl(calc_num_mtmts(fi)))
     349             : 
     350          12 :       mtmt_idx = 0
     351          32 :       DO itype = 1, fi%atoms%ntype
     352          52 :          DO ineq = 1, fi%atoms%neq(itype)
     353             :             ! Here the diagonal block matrices do not depend on ineq. In (1b) they do depend on ineq, though,
     354         140 :             DO l = 0, fi%hybinp%lcutm1(itype)
     355         100 :                mat%matsize1=mpdata%num_radbasfn(l, itype)
     356         100 :                mat%matsize2=mpdata%num_radbasfn(l, itype)
     357         860 :                DO n2 = 1, mpdata%num_radbasfn(l, itype)
     358             :                   ! note that mpdata%radbasfn_mt already contains the factor rgrid
     359             :                   CALL primitivef(primf1, mpdata%radbasfn_mt(:, n2, l, itype) &
     360             :                                     *fi%atoms%rmsh(:, itype)**(l + 1), fi%atoms%rmsh, fi%atoms%dx, &
     361      617920 :                                     fi%atoms%jri, fi%atoms%jmtd, itype, fi%atoms%ntype)
     362             :                   ! -itype is to enforce inward integration
     363             :                   CALL primitivef(primf2, mpdata%radbasfn_mt(:fi%atoms%jri(itype), n2, l, itype) &
     364             :                                     /fi%atoms%rmsh(:fi%atoms%jri(itype), itype)**l, fi%atoms%rmsh, fi%atoms%dx, &
     365      606832 :                                     fi%atoms%jri, fi%atoms%jmtd, -itype, fi%atoms%ntype)
     366             : 
     367      606072 :                   primf1(:fi%atoms%jri(itype)) = primf1(:fi%atoms%jri(itype))/fi%atoms%rmsh(:fi%atoms%jri(itype), itype)**l
     368      617160 :                   primf2 = primf2*fi%atoms%rmsh(:, itype)**(l + 1)
     369             : 
     370        4212 :                   DO n1 = 1, n2
     371     2708224 :                      integrand = mpdata%radbasfn_mt(:, n1, l, itype)*(primf1 + primf2)
     372        4112 :                      mat%data_r(n1, n2) = fpi_const/(2*l + 1) * intgrf(integrand, fi%atoms, itype, gridf)
     373             :                   END DO
     374             :                END DO
     375             : 
     376             :                ! distribute mat for m=-l,l on coulomb in block-matrix form
     377         620 :                DO M = -l, l
     378         500 :                   mtmt_idx = mtmt_idx + 1
     379         500 :                   call mtmt_repl(mtmt_idx)%alloc(.True., mpdata%num_radbasfn(l, itype), mpdata%num_radbasfn(l, itype))
     380             : 
     381        3980 :                   DO n2 = 1, mpdata%num_radbasfn(l, itype)
     382       18180 :                      DO n1 = 1, n2
     383       17680 :                         mtmt_repl(mtmt_idx)%data_r(n1, n2) = mat%data_r(n1, n2)
     384             :                      END DO
     385             :                   END DO
     386         600 :                   call mtmt_repl(mtmt_idx)%u2l()
     387             :                END DO
     388             : 
     389             :             END DO
     390             :          END DO
     391             :       END DO
     392          12 :       call mat%free()
     393          12 :       call timestop("loop 1")
     394             : 
     395          48 :       DO im = 1, size(fmpi%k_list)
     396          36 :          ikpt = fmpi%k_list(im)
     397             : 
     398             :          ! only the first rank handles the MT-MT part
     399          36 :          call timestart("MT-MT part")
     400          36 :          ix = 0
     401          36 :          ic2 = 0
     402          36 :          mtmt_idx = 0
     403          96 :          DO itype2 = 1, fi%atoms%ntype
     404         156 :             DO ineq2 = 1, fi%atoms%neq(itype2)
     405          60 :                ic2 = ic2 + 1
     406          60 :                lm2 = 0
     407         420 :                DO l2 = 0, fi%hybinp%lcutm1(itype2)
     408        1860 :                   DO m2 = -l2, l2
     409        1500 :                      mtmt_idx = mtmt_idx + 1
     410        1500 :                      lm2 = lm2 + 1
     411       12240 :                      DO n2 = 1, mpdata%num_radbasfn(l2, itype2)
     412       10440 :                         ix = ix + 1
     413       10440 :                         call glob_to_loc(fmpi, ix, pe, ix_loc)
     414       11940 :                         if(fmpi%n_rank == pe) then
     415             :                            iy = 0
     416             :                            ic1 = 0
     417       12420 :                            lp2: DO itype1 = 1, itype2
     418       19620 :                               DO ineq1 = 1, fi%atoms%neq(itype1)
     419        7200 :                                  ic1 = ic1 + 1
     420        7200 :                                  lm1 = 0
     421       50400 :                                  DO l1 = 0, fi%hybinp%lcutm1(itype1)
     422      223200 :                                     DO m1 = -l1, l1
     423      180000 :                                        lm1 = lm1 + 1
     424     1480206 :                                        DO n1 = 1, mpdata%num_radbasfn(l1, itype1)
     425     1264206 :                                           iy = iy + 1
     426     1264206 :                                           rdum = (-1)**(l2 + m2)*moment(n1, l1, itype1)*moment(n2, l2, itype2)*gmat(lm1, lm2)
     427     1264206 :                                           l = l1 + l2
     428     1264206 :                                           lm = l**2 + l + m1 - m2 + 1
     429             : 
     430     1264206 :                                           if(itype2 /= itype1 .or. ineq2 /= ineq1 .or. l2 /= l1 .or. m2 /= m1) then 
     431             :                                              mtmt_term = 0.0
     432             :                                           else 
     433       37380 :                                              mtmt_term = mtmt_repl(mtmt_idx)%data_r(n1, n2)
     434             :                                           endif
     435             :                                        
     436             :                                           coul(ikpt)%data_c(iy, ix_loc) = mtmt_term + EXP(CMPLX(0.0, 1.0)*tpi_const* &
     437             :                                                                      dot_PRODUCT(fi%kpts%bk(:, ikpt), &
     438             :                                                                                  fi%atoms%taual(:, ic2) - fi%atoms%taual(:, ic1))) &
     439     5236824 :                                                                *rdum*structconst(lm, ic1, ic2, ikpt)
     440             :                                        END DO
     441             :                                     END DO
     442             :                                  END DO
     443             :                               END DO
     444             :                            END DO lp2
     445             :                         endif ! pe = n_rank
     446             :                      END DO
     447             :                   END DO
     448             :                END DO
     449             :             END DO
     450             :          END DO
     451             : 
     452          36 :          IF (fi%sym%invs) THEN
     453             :             !symmetrize makes the Coulomb matrix real symmetric     
     454             :             CALL symmetrize_mpimat(fi, fmpi, coul(ikpt)%data_c, [1,1],[hybdat%nbasp, hybdat%nbasp], &
     455          72 :                                    3, .FALSE., mpdata%num_radbasfn)
     456             :          ENDIF
     457          48 :          call timestop("MT-MT part")
     458             :       END DO
     459             : 
     460         108 :       IF (maxval(mpdata%n_g) /= 0) THEN ! skip calculation of plane-wave contribution if mixed basis does not contain plane waves
     461             : 
     462             :          !
     463             :          !     (2) Case < MT | v | PW >
     464             :          !
     465             : 
     466             :          !     (2a) r in MT, r' everywhere
     467             :          !     (2b) r,r' in same MT
     468             :          !     (2c) r,r' in different MT
     469             : 
     470          12 :          call timestart("loop over interst.")
     471          48 :          DO im = 1, size(fmpi%k_list)
     472          36 :             ikpt = fmpi%k_list(im)
     473             :             call loop_over_interst(fi, hybdat, mpdata, fmpi, structconst, sphbesmoment, moment, moment2, &
     474          48 :                                    qnrm, facc, gmat, integral, olap, pqnrm, pgptm1, ngptm1, ikpt, coul(ikpt))
     475             : 
     476             :          END DO
     477             : 
     478          12 :          call timestop("loop over interst.")
     479          12 :          deallocate (olap, integral)
     480             : 
     481             :          !
     482             :          !     (3) Case < PW | v | PW >
     483             :          !
     484             :          !     (3a) r,r' everywhere; r everywhere, r' in MT; r in MT, r' everywhere
     485             : 
     486             :          ! Coulomb matrix, contribution (3a)
     487          12 :          call timestart("coulomb matrix 3a")
     488          48 :          DO im = 1, size(fmpi%k_list)
     489          36 :             ikpt = fmpi%k_list(im)
     490             : 
     491        1480 :             DO igpt0 = 1, ngptm1(ikpt)
     492        1432 :                igpt2 = pgptm1(igpt0, ikpt)
     493        1432 :                igptp2 = mpdata%gptm_ptr(igpt2, ikpt)
     494        1432 :                ix = hybdat%nbasp + igpt2
     495        1432 :                call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
     496        1468 :                if(fmpi%n_rank == pe_ix) then
     497       11456 :                   q2 = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp2), fi%cell%bmat)
     498        2864 :                   rdum2 = SUM(q2**2)
     499         716 :                   IF (abs(rdum2) > 1e-12) rdum2 = fpi_const/rdum2
     500             : 
     501       57126 :                   DO igpt1 = 1, igpt2
     502       56410 :                      igptp1 = mpdata%gptm_ptr(igpt1, ikpt)
     503       56410 :                      iy = hybdat%nbasp + igpt1
     504      902560 :                      q1 = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp1), fi%cell%bmat)
     505      225640 :                      rdum1 = SUM(q1**2)
     506       56410 :                      IF (abs(rdum1) > 1e-12) rdum1 = fpi_const/rdum1
     507       56410 :                      smat = calc_smat_elem(fi, mpdata, igptp1, igptp2)
     508             : 
     509       57126 :                      IF (ikpt == 1) THEN
     510        6542 :                         IF (igpt1 /= 1) THEN
     511        6456 :                            coul(1)%data_c(iy,ix_loc) = -smat*rdum1/fi%cell%vol 
     512             :                         END IF
     513        6542 :                         IF (igpt2 /= 1) THEN
     514             :                            coul(1)%data_c(iy,ix_loc) &
     515        6536 :                               = coul(1)%data_c(iy,ix_loc) - smat*rdum2/fi%cell%vol
     516             :                         END IF
     517             :                      ELSE
     518       49868 :                         coul(ikpt)%data_c(iy,ix_loc) = -smat*(rdum1 + rdum2)/fi%cell%vol
     519             :                      END IF
     520             :                   END DO
     521         716 :                   IF (ikpt /= 1 .OR. igpt2 /= 1) THEN
     522         710 :                      coul(ikpt)%data_c(iy,ix_loc) = coul(ikpt)%data_c(iy,ix_loc)  + rdum2
     523             :                   END IF
     524             :                endif
     525             :             END DO
     526             :          END DO
     527          12 :          call timestop("coulomb matrix 3a")
     528             :          !     (3b) r,r' in different MT
     529             : 
     530          12 :          call timestart("coulomb matrix 3b")
     531          48 :          DO im = 1, size(fmpi%k_list)
     532          36 :             ikpt = fmpi%k_list(im)
     533          36 :             if (fmpi%is_root()) write (*, *) "coulomb pw-loop nk: ("//int2str(ikpt)//"/"//int2str(fi%kpts%nkpt)//")"
     534             :             ! group together quantities which depend only on l,m and igpt -> carr2a
     535         828 :             allocate (carr2a((fi%hybinp%lexp + 1)**2, maxval(mpdata%n_g)), carr2b(fi%atoms%nat, maxval(mpdata%n_g)))
     536     1567512 :             carr2a = 0; carr2b = 0
     537        4952 :             DO igpt = 1, mpdata%n_g(ikpt)
     538        4916 :                igptp = mpdata%gptm_ptr(igpt, ikpt)
     539        4916 :                iqnrm = pqnrm(igpt, ikpt)
     540       78656 :                q = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%cell%bmat)
     541             : 
     542        4916 :                call ylm4(fi%hybinp%lexp, q, y)
     543             : 
     544     1425640 :                y = CONJG(y)
     545             :                lm = 0
     546       88488 :                DO l = 0, fi%hybinp%lexp
     547     1509212 :                   DO M = -l, l
     548     1420724 :                      lm = lm + 1
     549     1504296 :                      carr2a(lm, igpt) = fpi_const*CMPLX(0.0, 1.0)**(l)*y(lm)
     550             :                   END DO
     551             :                END DO
     552       14220 :                DO ic = 1, fi%atoms%nat
     553             :                   carr2b(ic, igpt) = EXP(-CMPLX(0.0, 1.0)*tpi_const* &
     554       41988 :                                          dot_PRODUCT(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%atoms%taual(:, ic)))
     555             :                END DO
     556             :             END DO
     557             : 
     558             :             !finally we can loop over the plane waves (G: igpt1,igpt2)
     559          36 :             call timestart("loop over plane waves")
     560             : 
     561             :             !$OMP parallel default(none) private(carr2, igpt0, igpt1, igpt2, ix, iy, ix_loc, pe_ix, iatom, lm1, l1, m1)&
     562             :             !$OMP private(lm2, l2, m2, cdum, ic, lm, l, m, itype, itype2, igptp1, csum, igptp2, iqnrm1, iqnrm2, cexp)&
     563             :             !$OMP shared(fmpi, fi, ikpt, ngptm1, pgptm1, pqnrm, coul, gmat, structconst, sphbesmoment, hybdat, mpdata)&
     564          36 :             !$OMP shared(carr2a, carr2b)
     565             :             allocate(carr2(fi%atoms%nat, (fi%hybinp%lexp + 1)**2))               
     566             :             !$OMP do schedule(dynamic)
     567             :             DO igpt0 = 1, ngptm1(ikpt)
     568             :                igpt2 = pgptm1(igpt0, ikpt)
     569             :                ix = hybdat%nbasp + igpt2
     570             :                call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
     571             :                if(fmpi%n_rank == pe_ix) then
     572             :                   igptp2 = mpdata%gptm_ptr(igpt2, ikpt)
     573             :                   iqnrm2 = pqnrm(igpt2, ikpt)
     574             : 
     575             :                   carr2 = 0
     576             : 
     577             : !                  call timestart("itype loops")
     578             :                   do iatom = 1,fi%atoms%nat
     579             :                      itype2 = fi%atoms%itype(iatom)
     580             :                      cexp = CONJG(carr2b(iatom, igpt2))
     581             :                      
     582             :                      DO lm1 = 1, (fi%hybinp%lexp+1)**2
     583             :                         call calc_l_m_from_lm(lm1, l1, m1)
     584             :                         do lm2 = 1, (fi%hybinp%lexp+1)**2
     585             :                            call calc_l_m_from_lm(lm2, l2, m2)
     586             :                            cdum = (-1)**(l2 + m2)*sphbesmoment(l2, itype2, iqnrm2)*cexp*carr2a(lm2, igpt2)*gmat(lm1, lm2)
     587             :                            l = l1 + l2
     588             :                            lm = l**2 + l - l1 - m2 + (m1 + l1) + 1
     589             :                            do iat2 =1,fi%atoms%nat
     590             :                               carr2(iat2, lm1) = carr2(iat2,lm1) + cdum*structconst(lm, iat2, iatom, ikpt)
     591             :                            enddo
     592             :                         enddo
     593             :                      enddo
     594             :                   end do ! iatom
     595             : 
     596             : !                  call timestop("itype loops")
     597             : 
     598             : !                  call timestart("igpt1")
     599             :                   DO igpt1 = 1, igpt2
     600             :                      iy = hybdat%nbasp + igpt1
     601             :                      igptp1 = mpdata%gptm_ptr(igpt1, ikpt)
     602             :                      iqnrm1 = pqnrm(igpt1, ikpt)
     603             :                      csum = 0
     604             :                      do ic = 1, fi%atoms%nat
     605             :                         do lm = 1, (fi%hybinp%lexp+1)**2
     606             :                            itype = fi%atoms%itype(ic)
     607             :                            call calc_l_m_from_lm(lm, l, m)
     608             :                            cdum = carr2b(ic, igpt1)*sphbesmoment(l, itype, iqnrm1)
     609             :                            csum = csum + cdum*carr2(ic, lm)*CONJG(carr2a(lm, igpt1)) ! for coulomb
     610             :                         END DO
     611             :                      END DO
     612             :                      coul(ikpt)%data_c(iy,ix_loc) = coul(ikpt)%data_c(iy,ix_loc) + csum/fi%cell%vol
     613             :                   END DO
     614             : !                  call timestop("igpt1")
     615             :                endif ! pe_ix
     616             :             END DO !igpt0
     617             :             !$omp end do
     618             :             deallocate (carr2) 
     619             :             !$OMP end parallel
     620          36 :             deallocate(carr2a, carr2b)
     621             : 
     622          48 :             call timestop("loop over plane waves")
     623             :          END DO !ikpt
     624          12 :          call timestop("coulomb matrix 3b")
     625             : 
     626             :          ! check if I own the gamma point
     627          12 :          if(any(fmpi%k_list == 1)) then
     628             :             !     Add corrections from higher orders in (3b) to coulomb(:,1)
     629             :             ! (1) igpt1 > 1 , igpt2 > 1  (finite G vectors)
     630          12 :             call timestart("add corrections from higher orders")
     631          12 :             call gamma_double_gpt_loop(fi, fmpi, hybdat, mpdata, sphbesmoment, gmat,  ngptm1, pgptm1, pqnrm, coul(1)%data_c) 
     632             : 
     633          12 :             rdum = (fpi_const)**(1.5)/fi%cell%vol**2*gmat(1, 1)
     634             : 
     635             :             ! (2) igpt1 = 1 , igpt2 > 1  (first G vector vanishes, second finite)
     636          12 :             call timestart("igpt1=1 loop")
     637          12 :             iy = hybdat%nbasp + 1
     638         184 :             DO igpt0 = 1, ngptm1(1)
     639         172 :                igpt2 = pgptm1(igpt0, 1)
     640         184 :                IF (igpt2 /= 1) then
     641         160 :                   ix = hybdat%nbasp + igpt2
     642         160 :                   call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
     643         160 :                   if(fmpi%n_rank == pe_ix) then
     644          80 :                      iqnrm2 = pqnrm(igpt2, 1)
     645          80 :                      igptp2 = mpdata%gptm_ptr(igpt2, 1)
     646          80 :                      qnorm2 = qnrm(iqnrm2)
     647         232 :                      DO itype1 = 1, fi%atoms%ntype
     648         384 :                         DO ineq1 = 1, fi%atoms%neq(itype1)
     649             :                            ic2 = 0
     650         600 :                            DO itype2 = 1, fi%atoms%ntype
     651         744 :                               DO ineq2 = 1, fi%atoms%neq(itype2)
     652         296 :                                  ic2 = ic2 + 1
     653        1184 :                                  cdum = EXP(CMPLX(0.0, 1.0)*tpi_const*dot_PRODUCT(mpdata%g(:, igptp2), fi%atoms%taual(:, ic2)))
     654             :                                  coul(1)%data_c(iy, ix_loc) = coul(1)%data_c(iy, ix_loc) &
     655             :                                                    + rdum*cdum*fi%atoms%rmt(itype1)**3*( &
     656             :                                                    +sphbesmoment(0, itype2, iqnrm2)/30*fi%atoms%rmt(itype1)**2 &
     657             :                                                    - sphbesmoment(2, itype2, iqnrm2)/18 &
     658         592 :                                                    + sphbesmoment(1, itype2, iqnrm2)/6/qnorm2)
     659             :                               END DO
     660             :                            END DO
     661             :                         END DO
     662             :                      END DO
     663             :                   endif !pe_ix
     664             :                endif
     665             :             END DO
     666          12 :             call timestop("igpt1=1 loop")
     667             : 
     668             :             ! (2) igpt1 = 1 , igpt2 = 1  (vanishing G vectors)
     669          12 :             call timestart("igpt1=igpt2=1 loop")
     670          12 :             iy = hybdat%nbasp + 1
     671          12 :             ix = hybdat%nbasp + 1
     672          12 :             call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
     673          12 :             if(pe_ix == fmpi%n_rank) then
     674          16 :                DO itype1 = 1, fi%atoms%ntype
     675          26 :                   DO ineq1 = 1, fi%atoms%neq(itype1)
     676          38 :                      DO itype2 = 1, fi%atoms%ntype
     677          46 :                         DO ineq2 = 1, fi%atoms%neq(itype2)
     678             :                            coul(1)%data_c(iy, ix_loc) = coul(1)%data_c(iy, ix_loc) &
     679             :                                              + rdum*fi%atoms%rmt(itype1)**3*fi%atoms%rmt(itype2)**3* &
     680          36 :                                              (fi%atoms%rmt(itype1)**2 + fi%atoms%rmt(itype2)**2)/90
     681             :                         END DO
     682             :                      END DO
     683             :                   END DO
     684             :                END DO
     685             :             endif ! pe_ix
     686          12 :             call timestop("igpt1=igpt2=1 loop")
     687          12 :             call timestop("add corrections from higher orders")
     688             :          endif
     689             : 
     690             :          !     (3c) r,r' in same MT
     691             : 
     692             :          ! Calculate sphbesintegral
     693          12 :          call timestart("sphbesintegral")
     694           0 :          allocate (sphbes0(-1:fi%hybinp%lexp + 2, fi%atoms%ntype, nqnrm),&
     695         192 :               &           carr2((fi%hybinp%lexp + 1)**2, maxval(mpdata%n_g)))
     696      533572 :          sphbes0 = 0; carr2 = 0
     697         420 :          DO iqnrm = 1, nqnrm
     698        1172 :             DO itype = 1, fi%atoms%ntype
     699         752 :                rdum = qnrm(iqnrm)*fi%atoms%rmt(itype)
     700         752 :                call sphbes(fi%hybinp%lexp + 2, rdum, sphbes0(0, itype, iqnrm))
     701        1160 :                IF (abs(rdum) > 1e-12) sphbes0(-1, itype, iqnrm) = COS(rdum)/rdum
     702             :             END DO
     703             :          END DO
     704          12 :          call timestop("sphbesintegral")
     705             : 
     706          12 :          call timestart("loop 2")
     707          48 :          DO im = 1, size(fmpi%k_list)
     708          36 :             ikpt = fmpi%k_list(im)
     709          36 :             call timestart("harmonics setup")
     710        4952 :             DO igpt = 1, mpdata%n_g(ikpt)
     711        4916 :                igptp = mpdata%gptm_ptr(igpt, ikpt)
     712       78656 :                q = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%cell%bmat)
     713        4952 :                call ylm4(fi%hybinp%lexp, q, carr2(:, igpt))
     714             :             END DO
     715          36 :             call timestop("harmonics setup")
     716             :             call perform_double_g_loop(fi, hybdat, fmpi, mpdata, sphbes0, carr2, ngptm1,pgptm1,&
     717          36 :                                        pqnrm,qnrm, nqnrm, ikpt, coul(ikpt))
     718             :             ! this one is needed
     719             : #ifdef CPP_MPI
     720          36 :             call timestart("post dblgloop barrier")
     721          36 :             call MPI_Barrier(fmpi%sub_comm, ierr)
     722          36 :             call timestop("post dblgloop barrier")
     723             : #endif
     724          84 :             call coul(ikpt)%u2l()
     725             :          END DO
     726          12 :          call timestop("loop 2")
     727          12 :          deallocate (carr2)
     728             : 
     729             :          !
     730             :          !     Symmetry-equivalent G vectors
     731             :          !
     732             :          ! All elements are needed so send all data to all processes treating the
     733             :          ! respective k-points
     734             : 
     735         252 :          allocate (carr2(maxval(hybdat%nbasm), 2), iarr(maxval(mpdata%n_g)))
     736             :          allocate (nsym_gpt(mpdata%num_gpts(), fi%kpts%nkpt), &
     737         144 :                    sym_gpt(MAXVAL(nsym1), mpdata%num_gpts(), fi%kpts%nkpt))
     738      258948 :          nsym_gpt = 0; sym_gpt = 0
     739          12 :          call timestart("loop 3")
     740          48 :          DO im = 1, size(fmpi%k_list)
     741          36 :             ikpt = fmpi%k_list(im)
     742       37044 :             carr2 = 0; iarr = 0
     743        1468 :             iarr(pgptm1(:ngptm1(ikpt), ikpt)) = 1
     744        1480 :             DO igpt0 = 1, ngptm1(ikpt)
     745        1432 :                lsym = (1 <= igpt0) .AND. (ngptm1(ikpt) >= igpt0)
     746        1432 :                igpt2 = pgptm1(igpt0, ikpt)
     747        1432 :                ix = hybdat%nbasp + igpt2
     748        1432 :                call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
     749        1432 :                if(pe_ix == fmpi%n_rank) then 
     750      355344 :                   carr2(:hybdat%nbasm(ikpt),2) = coul(ikpt)%data_c(:hybdat%nbasm(ikpt),ix_loc)
     751             :                endif
     752             : #ifdef CPP_MPI
     753        1432 :                call timestart("bcast carr2")
     754        1432 :                call MPI_Bcast(carr2(1,2), hybdat%nbasm(ikpt), MPI_DOUBLE_COMPLEX, pe_ix, fmpi%sub_comm, ierr)
     755        1432 :                call timestop("bcast carr2")
     756             : #endif
     757             : 
     758        1432 :                IF (lsym) THEN
     759        1432 :                   ic = 1
     760        1432 :                   sym_gpt(ic, igpt0, ikpt) = igpt2
     761             :                END IF
     762        9504 :                DO isym1 = 2, nsym1(ikpt)
     763        8072 :                   isym = sym1(isym1, ikpt)
     764             :                   CALL bramat_trafo(carr2(:, 2), igpt2, ikpt, isym, .FALSE., POINTER(ikpt, :, :, :), &
     765             :                                     fi%sym, rrot(:, :, isym), invrrot(:, :, isym), mpdata, fi%hybinp, &
     766             :                                     fi%kpts, maxval(fi%hybinp%lcutm1), fi%atoms, fi%hybinp%lcutm1, &
     767             :                                     mpdata%num_radbasfn, maxval(mpdata%num_radbasfn), dwgn(:, :, :, isym), &
     768     9830032 :                                     hybdat%nbasp, hybdat%nbasm, carr2(:, 1), igpt1)
     769        9504 :                   IF (iarr(igpt1) == 0) THEN
     770             :                      CALL bramat_trafo(carr2(:, 2), igpt2, ikpt, isym, .TRUE., POINTER(ikpt, :, :, :), &
     771             :                                        fi%sym, rrot(:, :, isym), invrrot(:, :, isym), mpdata, fi%hybinp, &
     772             :                                        fi%kpts, maxval(fi%hybinp%lcutm1), fi%atoms, fi%hybinp%lcutm1, &
     773             :                                        mpdata%num_radbasfn, maxval(mpdata%num_radbasfn), &
     774     4311088 :                                        dwgn(:, :, :, isym), hybdat%nbasp, hybdat%nbasm, carr2(:, 1), igpt1)
     775        3484 :                      l = (hybdat%nbasp + igpt1 - 1)*(hybdat%nbasp + igpt1)/2
     776        3484 :                      ix = hybdat%nbasp + igpt1
     777        3484 :                      call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
     778        3484 :                      if(pe_ix == fmpi%n_rank) then
     779      710330 :                         coul(ikpt)%data_c(:hybdat%nbasp + igpt1,ix_loc) = carr2(:hybdat%nbasp + igpt1, 1)
     780             :                      endif 
     781             : 
     782     1420660 :                      do ix = 1,hybdat%nbasp + igpt1
     783     1417176 :                         call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
     784     1420660 :                         if(pe_ix == fmpi%n_rank) coul(ikpt)%data_c(hybdat%nbasp + igpt1, ix_loc) = conjg(carr2(ix, 1))
     785             :                      enddo
     786             : 
     787        3484 :                      iarr(igpt1) = 1
     788        3484 :                      IF (lsym) THEN
     789        3484 :                         ic = ic + 1
     790        3484 :                         sym_gpt(ic, igpt0, ikpt) = igpt1
     791             :                      END IF
     792             :                   END IF
     793             :                END DO
     794        1468 :                nsym_gpt(igpt0, ikpt) = ic
     795             :             END DO ! igpt0
     796             :          END DO ! ikpt
     797          12 :          call timestop("loop 3")
     798          12 :          call timestart("gap 1:")
     799          12 :          deallocate (carr2, iarr, pgptm1)
     800             :       END IF
     801          12 :       deallocate (qnrm, pqnrm)
     802             : 
     803          12 :       IF (xcpot%is_name("hse") .OR. xcpot%is_name("vhse")) THEN
     804             :          !
     805             :          ! The HSE functional is realized subtracting erf/r from
     806             :          ! the normal Coulomb matrix
     807             :          !
     808             :       ELSE
     809             :          ! check for gamma
     810          12 :          if(any(fmpi%k_list == 1)) then
     811             :             CALL subtract_sphaverage(fi%sym, fi%cell, fi%atoms, mpdata, &
     812          12 :                                     fi%hybinp, hybdat, fmpi, hybdat%nbasm, gridf, coul(1))
     813             :          endif
     814             :       END IF
     815             :       
     816             :       ! transform Coulomb matrix to the biorthogonal set
     817          12 :       call timestop("gap 1:")
     818          48 :       DO im = 1, size(fmpi%k_list)
     819          36 :          ikpt = fmpi%k_list(im)
     820          36 :          call apply_inverse_olaps(mpdata, fi%atoms, fi%cell, hybdat, fmpi, fi%sym, ikpt, coul(ikpt))
     821             :          ! lower to upper, because the lower half is better in memory
     822             : #ifdef CPP_MPI
     823          36 :          call timestart("post inverse barrier")
     824          36 :          call MPI_BARRIER(fmpi%sub_comm, ierr)
     825          36 :          call timestop("post inverse barrier")
     826             : #endif
     827          84 :          call coul(ikpt)%l2u()
     828             :       enddo
     829             : 
     830             :       !call plot_coulombmatrix() -> code was shifted to plot_coulombmatrix.F90
     831             :       !
     832             :       ! rearrange coulomb matrix
     833             :       !
     834          12 :       if(.not. allocated(hybdat%coul)) allocate(hybdat%coul(fi%kpts%nkpt))
     835             : 
     836          48 :       DO im = 1, size(fmpi%k_list)
     837          36 :          ikpt = fmpi%k_list(im)
     838             :          ! unpack coulomb into coulomb(ikpt)
     839          48 :          call copy_from_dense_to_sparse(fi, fmpi, mpdata, coul, ikpt, hybdat)
     840             :       END DO ! ikpt
     841          12 :       CALL timestop("Coulomb matrix setup")
     842             : 
     843         524 :    END SUBROUTINE coulombmatrix
     844             : 
     845             :    !     Calculate body of Coulomb matrix at Gamma point: v_IJ = SUM(G) c^*_IG c_JG 4*pi/G**2 .
     846             :    !     For this we must subtract from coulomb(:,1) the spherical average of a term that comes
     847             :    !     from the fact that MT functions have k-dependent Fourier coefficients (see script).
     848          12 :    SUBROUTINE subtract_sphaverage(sym, cell, atoms, mpdata, hybinp, hybdat, fmpi, nbasm1, gridf, coulomb)
     849             : 
     850             :       USE m_types
     851             :       USE m_constants
     852             :       USE m_wrapper
     853             :       USE m_trafo
     854             :       USE m_util
     855             :       use m_intgrf
     856             :       USE m_olap
     857             :       IMPLICIT NONE
     858             : 
     859             :       TYPE(t_sym), INTENT(IN)    :: sym
     860             :       TYPE(t_cell), INTENT(IN)    :: cell
     861             :       TYPE(t_atoms), INTENT(IN)    :: atoms
     862             :       TYPE(t_mpdata), intent(in)  :: mpdata
     863             :       TYPE(t_hybinp), INTENT(IN)    :: hybinp
     864             :       TYPE(t_hybdat), INTENT(IN)    :: hybdat
     865             :       type(t_mpi), intent(in)    :: fmpi
     866             : 
     867             :       INTEGER, INTENT(IN)    :: nbasm1(:)
     868             :       REAL, INTENT(IN)    :: gridf(:, :)
     869             :       class(t_mat), intent(inout) :: coulomb
     870             : 
     871             :       ! - local scalars -
     872             :       INTEGER               :: l, ix, iy, ix_loc, pe_ix, i, j, n, nn, itype, ieq, M, ierr
     873             : 
     874             :       ! - local arrays -
     875          12 :       TYPE(t_mat) :: olap
     876             :       !COMPLEX , ALLOCATABLE :: constfunc(:)  !can also be real in inversion case
     877          12 :       COMPLEX      :: coeff(1,nbasm1(1)), cderiv(-1:1, nbasm1(1)), claplace(1,nbasm1(1))
     878             : 
     879          12 :       call timestart("subtract_sphaverage")
     880          12 :       CALL olap%alloc(sym%invs, mpdata%n_g(1), mpdata%n_g(1), 0.)
     881             : 
     882          12 :       n = nbasm1(1)
     883          12 :       nn = n*(n + 1)/2
     884        8124 :       CALL olap_pw(olap, mpdata%g(:, mpdata%gptm_ptr(:mpdata%n_g(1), 1)), mpdata%n_g(1), atoms, cell, fmpi)
     885             : 
     886             :       ! Define coefficients (coeff) and their derivatives (cderiv,claplace)
     887       10212 :       coeff = 0
     888       20412 :       cderiv = 0
     889       10212 :       claplace = 0
     890          12 :       j = 0
     891          32 :       DO itype = 1, atoms%ntype
     892          52 :          DO ieq = 1, atoms%neq(itype)
     893         140 :             DO l = 0, hybinp%lcutm1(itype)
     894         620 :                DO M = -l, l
     895        4080 :                   DO i = 1, mpdata%num_radbasfn(l, itype)
     896        3480 :                      j = j + 1
     897        3980 :                      IF (l == 0) THEN
     898             :                         coeff(1,j) = SQRT(fpi_const) &
     899             :                                    *intgrf(atoms%rmsh(:, itype)*mpdata%radbasfn_mt(:, i, 0, itype), &
     900             :                                            atoms, itype, gridf) &
     901      146492 :                                    /SQRT(cell%vol)
     902             : 
     903             :                         claplace(1,j) = -SQRT(fpi_const) &
     904             :                                       *intgrf(atoms%rmsh(:, itype)**3*mpdata%radbasfn_mt(:, i, 0, itype), &
     905             :                                               atoms, itype, gridf) &
     906      146492 :                                       /SQRT(cell%vol)
     907             : 
     908        3300 :                      ELSE IF (l == 1) THEN
     909             :                         cderiv(M,j) = -SQRT(fpi_const/3)*CMPLX(0.0, 1.0) &
     910             :                                        *intgrf(atoms%rmsh(:, itype)**2*mpdata%radbasfn_mt(:, i, 1, itype), &
     911             :                                                atoms, itype, gridf) &
     912      408240 :                                        /SQRT(cell%vol)
     913             :                      END IF
     914             :                   END DO
     915             :                END DO
     916             :             END DO
     917             :          END DO
     918             :       END DO
     919          12 :       IF (olap%l_real) THEN
     920         952 :          coeff(1,hybdat%nbasp + 1:n) = olap%data_r(1, 1:n - hybdat%nbasp)
     921             :       else
     922         680 :          coeff(1,hybdat%nbasp + 1:n) = olap%data_c(1, 1:n - hybdat%nbasp)
     923             :       END IF
     924          12 :       IF (sym%invs) THEN
     925             :          CALL symmetrize(coeff, 1, nbasm1(1), 2, &
     926             :                          atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
     927          20 :                          mpdata%num_radbasfn, sym)
     928             :          CALL symmetrize(claplace, 1, nbasm1(1), 2, &
     929             :                          atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
     930          20 :                          mpdata%num_radbasfn, sym)
     931             :          CALL symmetrize(cderiv(-1:-1,:), 1, nbasm1(1), 2, &
     932             :                          atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
     933          20 :                          mpdata%num_radbasfn, sym)
     934             :          CALL symmetrize(cderiv(0:0,:), 1, nbasm1(1), 2, &
     935             :                          atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
     936          20 :                          mpdata%num_radbasfn, sym)
     937             :          CALL symmetrize(cderiv(1:1,:), 1, nbasm1(1), 2, &
     938             :                          atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
     939          20 :                          mpdata%num_radbasfn, sym)
     940             :       ENDIF
     941             :       ! Subtract head contributions from coulomb(:nn,1) to obtain the body
     942          12 :       l = 0
     943        5112 :       DO ix = 1, n
     944        5100 :          call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
     945        5112 :          if(fmpi%n_rank == pe_ix) then
     946      597218 :             DO iy = 1, ix
     947      594668 :                l = l + 1
     948             :                coulomb%data_c(iy,ix_loc) = coulomb%data_c(iy,ix_loc) - fpi_const/3 &
     949             :                                           *(dot_PRODUCT(cderiv(:,iy), cderiv(:,ix)) &
     950             :                                           + (CONJG(coeff(1,iy))*claplace(1,ix) &
     951     2381222 :                                              + CONJG(claplace(1,iy))*coeff(1,ix))/2)
     952             :             END DO
     953             :          endif
     954             :       END DO
     955             : 
     956             :       !needed bc apply inverse uses lower half 
     957             : #ifdef CPP_MPI
     958          12 :       call timestart("post subtr avg barrier")
     959          12 :       call MPI_Barrier(fmpi%sub_comm, ierr)
     960          12 :       call timestop("post subtr avg barrier")
     961             : #endif
     962          12 :       call coulomb%u2l()
     963          12 :       call timestop("subtract_sphaverage")
     964          12 :    END SUBROUTINE subtract_sphaverage
     965             : 
     966             : 
     967             : 
     968             :    !     ---------
     969             : 
     970             :    !     Returns a list of (k+G) vector lengths in qnrm(1:nqnrm) and the corresponding pointer pqnrm(1:ngpt(ikpt),ikpt)
     971          12 :    SUBROUTINE getnorm(kpts, gpt, ngpt, pgpt, qnrm, nqnrm, pqnrm, cell)
     972             :       USE m_types
     973             :       USE m_juDFT
     974             :       IMPLICIT NONE
     975             :       TYPE(t_cell), INTENT(IN)   :: cell
     976             :       TYPE(t_kpts), INTENT(IN)   :: kpts
     977             : 
     978             :       INTEGER, INTENT(IN)  :: ngpt(:), gpt(:, :), pgpt(:, :)!(dim,kpts%nkpt)
     979          12 :       REAL, ALLOCATABLE :: qnrm(:), help(:)
     980             :       INTEGER, INTENT(INOUT) :: nqnrm
     981             :       INTEGER, ALLOCATABLE :: pqnrm(:, :)
     982             :       INTEGER               :: i, j, ikpt, igpt, igptp
     983             :       REAL                  :: q(3), qnorm
     984             : 
     985        5484 :       allocate (qnrm(MAXVAL(ngpt)*kpts%nkpt), source=0.0)
     986        5532 :       allocate (pqnrm(MAXVAL(ngpt), kpts%nkpt), source=0)
     987             :       i = 0
     988          48 :       DO ikpt = 1, kpts%nkpt
     989        4964 :          igptloop: DO igpt = 1, ngpt(ikpt)
     990        4916 :             igptp = pgpt(igpt, ikpt)
     991        4916 :             IF (igptp == 0) call judft_error('getnorm: zero pointer (bug?)')
     992       78656 :             q = MATMUL(kpts%bk(:, ikpt) + gpt(:, igptp), cell%bmat)
     993       19664 :             qnorm = norm2(q)
     994      104060 :             DO j = 1, i
     995      104060 :                IF (ABS(qnrm(j) - qnorm) < 1e-12) THEN
     996        4508 :                   pqnrm(igpt, ikpt) = j
     997        4508 :                   CYCLE igptloop
     998             :                END IF
     999             :             END DO
    1000         408 :             i = i + 1
    1001         408 :             qnrm(i) = qnorm
    1002         444 :             pqnrm(igpt, ikpt) = i
    1003             :          END DO igptloop
    1004             :       END DO
    1005          12 :       nqnrm = i
    1006             : 
    1007          36 :       allocate (help(nqnrm))
    1008         420 :       help(1:nqnrm) = qnrm(1:nqnrm)
    1009          12 :       deallocate (qnrm)
    1010          36 :       allocate (qnrm(1:nqnrm))
    1011         432 :       qnrm = help
    1012             : 
    1013          12 :    END SUBROUTINE getnorm
    1014             : 
    1015         324 :    subroutine loop_over_interst(fi, hybdat, mpdata, fmpi, structconst, sphbesmoment, moment, moment2, &
    1016          36 :                                 qnrm, facc, gmat, integral, olap, pqnrm, pgptm1, ngptm1, ikpt, coul)
    1017             :       use m_types
    1018             :       use m_juDFT
    1019             :       use m_ylm, only: ylm4
    1020             :       use m_constants, only: fpi_const, tpi_const
    1021             :       USE m_trafo, ONLY: symmetrize
    1022             :       use m_calc_l_m_from_lm
    1023             :       implicit none
    1024             : 
    1025             :       type(t_fleurinput), intent(in)    :: fi
    1026             :       type(t_hybdat), intent(in)        :: hybdat
    1027             :       type(t_mpdata), intent(in)        :: mpdata
    1028             :       type(t_mpi), intent(in)           :: fmpi
    1029             :       REAL, intent(in)                  :: sphbesmoment(0:, :, :), qnrm(:), facC(-1:), gmat(:, :), moment(:, 0:, :), moment2(:, :)
    1030             :       real, intent(in)                  :: integral(:, 0:, :, :), olap(:, 0:, :, :)
    1031             :       integer, intent(in)               :: ikpt, ngptm1(:), pqnrm(:, :), pgptm1(:, :)
    1032             :       complex, intent(in)               :: structconst(:, :, :, :)
    1033             :       class(t_mat), intent(inout)       :: coul
    1034             : 
    1035             :       integer  :: igpt0, igpt, igptp, iqnrm, niter
    1036             :       integer  :: ix, iy, ic, itype, lm, l, m, itype1, ic1, l1, m1, lm1, loc_from
    1037             :       integer  :: l2, m2, lm2, n, i, iatm, j_type, j_l, iy_start, j_m, j_lm, pe_ix, ix_loc
    1038             :       real     :: q(3), qnorm, svol, tmp_vec(3)
    1039          36 :       COMPLEX  :: y((fi%hybinp%lexp + 1)**2), y1((fi%hybinp%lexp + 1)**2), y2((fi%hybinp%lexp + 1)**2)
    1040             :       complex  :: csum, csumf(9), cdum, cexp
    1041          36 :       integer, allocatable :: lm_arr(:), ic_arr(:)
    1042             : 
    1043             : 
    1044          36 :       call range_from_glob_to_loc(fmpi, hybdat%nbasp+1, loc_from)
    1045      795666 :       coul%data_c(:hybdat%nbasp,loc_from:) = 0 
    1046             : 
    1047          36 :       svol = SQRT(fi%cell%vol)
    1048             :       ! start to loop over interstitial plane waves
    1049             :       !DO igpt0 = 1, ngptm1(ikpt)
    1050        1468 :       do igpt0 = 1, ngptm1(ikpt)
    1051        1432 :          igpt = pgptm1(igpt0, ikpt)
    1052        1432 :          igptp = mpdata%gptm_ptr(igpt, ikpt)
    1053        1432 :          ix = hybdat%nbasp + igpt
    1054        1432 :          call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
    1055        1468 :          if(pe_ix == fmpi%n_rank) then 
    1056       11456 :             q = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%cell%bmat)
    1057        2864 :             qnorm = norm2(q)
    1058         716 :             iqnrm = pqnrm(igpt, ikpt)
    1059         716 :             IF (ABS(qnrm(iqnrm) - qnorm) > 1e-12) then
    1060           0 :                call judft_error('coulombmatrix: qnorm does not equal corresponding & element in qnrm (bug?)') ! We shouldn't st op here!
    1061             :             endif
    1062             : 
    1063        9308 :             tmp_vec = MATMUL(fi%kpts%bk(:, fi%kpts%nkpt), fi%cell%bmat)
    1064         716 :             call ylm4(2, tmp_vec, y1)
    1065       11456 :             tmp_vec = MATMUL(mpdata%g(:, igptp), fi%cell%bmat)
    1066         716 :             call ylm4(2, tmp_vec, y2)
    1067         716 :             call ylm4(fi%hybinp%lexp, q, y)
    1068      621488 :             y1 = CONJG(y1); y2 = CONJG(y2); y = CONJG(y)
    1069             : 
    1070             :             ! this unrolls the do ic=1,atoms%nat{do lm=1,..{}} 
    1071         716 :             call collapse_ic_and_lm_loop(fi%atoms, fi%hybinp%lcutm1, niter, ic_arr, lm_arr)
    1072             : 
    1073             :             !$OMP PARALLEL DO default(none) &
    1074             :             !$OMP private(ic, lm, itype, l, m, csum, csumf, ic1, itype1, cexp, lm1, l2, cdum, m2, lm2, iy) &
    1075             :             !$OMP private(j_m, j_type, iy_start, l1, m1) &
    1076             :             !$OMP shared(ic_arr, lm_arr, fi, mpdata, olap, qnorm, moment, integral, hybdat, svol) &
    1077             :             !$OMP shared(moment2, ix, igpt, facc, structconst, y, y1, y2, gmat, iqnrm, sphbesmoment, ikpt) &
    1078             :             !$OMP shared(igptp, niter, fmpi, pe_ix, coul, ix_loc) &
    1079         716 :             !$OMP schedule(dynamic)
    1080             :             do i = 1,niter 
    1081             :                ic = ic_arr(i)
    1082             :                lm = lm_arr(i)
    1083             : 
    1084             :                itype = fi%atoms%itype(ic)
    1085             :                call calc_l_m_from_lm(lm, l, m)
    1086             : 
    1087             :                ! calculate sum over lm and centers for (2c) -> csum, csumf
    1088             :                csum = 0
    1089             :                csumf = 0
    1090             :                do ic1 = 1, fi%atoms%nat
    1091             :                   itype1 = fi%atoms%itype(ic1)
    1092             :                   cexp = fpi_const*EXP(CMPLX(0.0, 1.0)*tpi_const &
    1093             :                                        *(dot_PRODUCT(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%atoms%taual(:, ic1)) &
    1094             :                                           - dot_PRODUCT(fi%kpts%bk(:, ikpt), fi%atoms%taual(:, ic))))
    1095             : 
    1096             :                   do lm1 = 1, (fi%hybinp%lexp+1)**2
    1097             :                      call calc_l_m_from_lm(lm1, l1, m1)
    1098             :                      l2 = l + l1 ! for structconst
    1099             :                      cdum = sphbesmoment(l1, itype1, iqnrm)*CMPLX(0.0, 1.0)**(l1)*cexp
    1100             :                      m2 = M - m1              ! for structconst
    1101             :                      lm2 = l2**2 + l2 + m2 + 1 !
    1102             :                      csum = csum - (-1)**(m1 + l1)*gmat(lm1, lm)*y(lm1)*cdum*structconst(lm2, ic, ic1, ikpt)
    1103             :                   END DO
    1104             : 
    1105             :                   ! add contribution of (2c) to csum and csumf coming from linear and quadratic orders of Y_lm*(G) / G * j_(l+1)(GS)
    1106             :                   IF (ikpt == 1 .AND. l <= 2) THEN
    1107             :                      cexp = EXP(CMPLX(0.0, 1.0)*tpi_const*dot_PRODUCT(mpdata%g(:, igptp), fi%atoms%taual(:, ic1))) &
    1108             :                               *gmat(lm, 1)*fpi_const/fi%cell%vol
    1109             :                      csumf(lm) = csumf(lm) - cexp*SQRT(fpi_const)* &
    1110             :                                  CMPLX(0.0, 1.0)**l*sphbesmoment(0, itype1, iqnrm)/facC(l - 1)
    1111             :                      IF (l == 0) THEN
    1112             :                         IF (igpt /= 1) THEN
    1113             :                            csum = csum - cexp*(sphbesmoment(0, itype1, iqnrm)*fi%atoms%rmt(itype1)**2 - &
    1114             :                                                 sphbesmoment(2, itype1, iqnrm)*2.0/3)/10
    1115             :                         ELSE
    1116             :                            csum = csum - cexp*fi%atoms%rmt(itype1)**5/30
    1117             :                         END IF
    1118             :                      ELSE IF (l == 1) THEN
    1119             :                         csum = csum + cexp*CMPLX(0.0, 1.0)*SQRT(fpi_const) &
    1120             :                                  *sphbesmoment(1, itype1, iqnrm)*y(lm)/3
    1121             :                      END IF
    1122             :                   END IF
    1123             :                END DO
    1124             : 
    1125             :                ! add contribution of (2a) to csumf
    1126             :                IF (ikpt == 1 .AND. igpt == 1 .AND. l <= 2) THEN
    1127             :                   csumf(lm) = csumf(lm) + (fpi_const)**2*CMPLX(0.0, 1.0)**l/facC(l)
    1128             :                END IF
    1129             : 
    1130             :                ! finally define coulomb
    1131             :                cdum = (fpi_const)**2*CMPLX(0.0, 1.0)**(l)*y(lm) &
    1132             :                         *EXP(CMPLX(0.0, 1.0)*tpi_const &
    1133             :                            *dot_PRODUCT(mpdata%g(:, igptp), fi%atoms%taual(:, ic)))
    1134             : 
    1135             :                !calclate iy_start on the fly for OpenMP
    1136             :                iy_start = 0
    1137             :                do iatm = 1, ic-1
    1138             :                   j_type = fi%atoms%itype(iatm)
    1139             :                   do j_l = 0,fi%hybinp%lcutm1(j_type)
    1140             :                      iy_start = iy_start + mpdata%num_radbasfn(j_l, j_type) * (2*j_l+1)
    1141             :                   end do
    1142             :                end do
    1143             :                do j_lm = 1,lm-1
    1144             :                   call calc_l_m_from_lm(j_lm, j_l, j_m)
    1145             :                   iy_start = iy_start + mpdata%num_radbasfn(j_l, itype) 
    1146             :                enddo
    1147             : 
    1148             :                DO n = 1, mpdata%num_radbasfn(l, itype)
    1149             :                   iy = iy_start + n
    1150             : 
    1151             :                   IF (ikpt == 1 .AND. igpt == 1) THEN
    1152             :                      IF (l == 0) coul%data_c(iy, ix_loc) = -cdum*moment2(n, itype)/6/svol 
    1153             :                      coul%data_c(iy, ix_loc) = coul%data_c(iy, ix_loc) &
    1154             :                                                       + (-cdum/(2*l + 1)*integral(n, l, itype, iqnrm) & ! (2b)&
    1155             :                                                          + csum*moment(n, l, itype))/svol          ! (2c)
    1156             :                   ELSE
    1157             :                      coul%data_c(iy, ix_loc) = (cdum*olap(n, l, itype, iqnrm)/qnorm**2 &  ! (2a)&
    1158             :                                                 - cdum/(2*l + 1)*integral(n, l, itype, iqnrm) & ! (2b)&
    1159             :                                                 + csum*moment(n, l, itype))/svol          ! (2c)
    1160             :                   endif
    1161             :                END DO
    1162             :             END DO ! collapsed atom & lm loop (ic)
    1163             :             !$OMP END PARALLEL DO
    1164             :          endif !pe_ix
    1165             :       END DO
    1166             : 
    1167          36 :       IF (fi%sym%invs) THEN
    1168             :          call symmetrize_mpimat(fi, fmpi, coul%data_c, [1,hybdat%nbasp+1], [hybdat%nbasp, hybdat%nbasp+mpdata%n_g(ikpt)],&
    1169         120 :                                 1, .false., mpdata%num_radbasfn)
    1170             :       ENDIF
    1171          36 :    endsubroutine loop_over_interst
    1172             : 
    1173          36 :    subroutine perform_double_g_loop(fi, hybdat, fmpi, mpdata, sphbes0, carr2, ngptm1,pgptm1,pqnrm,qnrm, nqnrm, ikpt, coulomb)
    1174             :       use m_juDFT
    1175             :       use m_constants, only: tpi_const,fpi_const
    1176             :       use m_sphbessel_integral
    1177             :       implicit none
    1178             :       type(t_fleurinput), intent(in)    :: fi
    1179             :       TYPE(t_mpdata), intent(in)        :: mpdata
    1180             :       TYPE(t_hybdat), INTENT(IN)        :: hybdat
    1181             :       type(t_mpi), intent(in)           :: fmpi
    1182             :       integer, intent(in)               :: ikpt, ngptm1(:), pqnrm(:,:),pgptm1(:, :), nqnrm
    1183             :       real, intent(in)                  :: qnrm(:), sphbes0(:, :, :)
    1184             :       complex, intent(in)               :: carr2(:, :)
    1185             :       !complex, intent(inout)            :: coulomb(:) ! only at ikpt
    1186             :       class(t_mat), intent(inout)        :: coulomb
    1187             : 
    1188             :       integer :: igpt0, igpt1, igpt2, ix, iy, igptp1, igptp2, iqnrm1, iqnrm2
    1189             :       integer :: ic, itype, lm, m, l, pe_ix, ix_loc
    1190             :       real    :: q1(3), q2(3)
    1191          36 :       complex :: y1((fi%hybinp%lexp + 1)**2), y2((fi%hybinp%lexp + 1)**2)
    1192          36 :       COMPLEX :: cexp1(fi%atoms%ntype)
    1193             :       complex :: cdum, cdum1 
    1194             :       logical :: ldum
    1195             : 
    1196          36 :       call timestart("double g-loop")
    1197             : 
    1198        1468 :       DO igpt0 = 1, ngptm1(ikpt)
    1199        1432 :          igpt2 = pgptm1(igpt0, ikpt)
    1200        1432 :          ix = hybdat%nbasp + igpt2
    1201        1432 :          call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
    1202        1468 :          if(pe_ix == fmpi%n_rank) then
    1203         716 :             igptp2 = mpdata%gptm_ptr(igpt2, ikpt)
    1204         716 :             iqnrm2 = pqnrm(igpt2, ikpt)
    1205       11456 :             q2 = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp2), fi%cell%bmat)
    1206      207640 :             y2 = CONJG(carr2(:, igpt2))
    1207             :             
    1208             :             !$OMP PARALLEL DO default(none) &
    1209             :             !$OMP private(igpt1, iy, igptp1, iqnrm1, q1, y1, cexp1, ic, itype, lm) &
    1210             :             !$OMP private(cdum, l, cdum1, m, ldum) &
    1211             :             !$OMP shared(igpt2, coulomb, hybdat, mpdata, ikpt, fi, carr2, pqnrm, igptp2)&
    1212         716 :             !$OMP shared(qnrm, sphbes0, iqnrm2, nqnrm, y2, ix_loc)
    1213             :             DO igpt1 = 1, igpt2
    1214             :                iy = hybdat%nbasp + igpt1
    1215             :                igptp1 = mpdata%gptm_ptr(igpt1, ikpt)
    1216             :                iqnrm1 = pqnrm(igpt1, ikpt)
    1217             :                q1 = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp1), fi%cell%bmat)
    1218             :                y1 = carr2(:, igpt1)
    1219             :                cexp1 = 0
    1220             :                do ic = 1,fi%atoms%nat 
    1221             :                   itype = fi%atoms%itype(ic)
    1222             :                   cexp1(itype) = cexp1(itype) + &
    1223             :                                  EXP(CMPLX(0.0, 1.0)*tpi_const*dot_PRODUCT( &
    1224             :                                        (mpdata%g(:, igptp2) - mpdata%g(:, igptp1)), fi%atoms%taual(:, ic)))
    1225             :                ENDDO
    1226             :                lm = 0
    1227             :                cdum = 0
    1228             :                DO l = 0, fi%hybinp%lexp
    1229             :                   cdum1 = 0
    1230             :                   DO itype = 1, fi%atoms%ntype
    1231             :                      cdum1 = cdum1 + cexp1(itype)*sphbessel_integral( &
    1232             :                               fi%atoms, itype, qnrm, nqnrm, &
    1233             :                               iqnrm1, iqnrm2, l, fi%hybinp, &
    1234             :                               sphbes0, .False., ldum) &
    1235             :                               /(2*l + 1)
    1236             :                   END DO
    1237             :                   DO M = -l, l
    1238             :                      lm = lm + 1
    1239             :                      cdum = cdum + cdum1*y1(lm)*y2(lm)
    1240             :                   ENDDO
    1241             :                ENDDO
    1242             :                coulomb%data_c(iy,ix_loc) = coulomb%data_c(iy,ix_loc) + (fpi_const)**3*cdum/fi%cell%vol
    1243             :             END DO
    1244             :             !$OMP end parallel do
    1245             :          endif !pe_ix
    1246             :       END DO
    1247          36 :       call timestop("double g-loop")
    1248          36 :    end subroutine perform_double_g_loop
    1249             : 
    1250         716 :    subroutine collapse_ic_and_lm_loop(atoms, lcutm1, niter, ic_arr, lm_arr)
    1251             :       use m_types
    1252             :       implicit none 
    1253             :       type(t_atoms), intent(in) :: atoms 
    1254             :       integer, intent(in)       :: lcutm1(:)
    1255             :       integer, intent(out)      :: niter 
    1256             :       integer, intent(inout), allocatable :: ic_arr(:),  lm_arr(:) 
    1257             : 
    1258             :       integer :: ic, lm, itype 
    1259             : 
    1260         716 :       if(allocated(ic_arr)) deallocate(ic_arr)
    1261         716 :       if(allocated(lm_arr)) deallocate(lm_arr)
    1262             : 
    1263         716 :       niter = 0
    1264        2082 :       do ic = 1, atoms%nat
    1265        1366 :          itype = atoms%itype(ic)
    1266       36232 :          do lm = 1,(lcutm1(itype)+1)**2
    1267       35516 :             niter = niter + 1
    1268             :          enddo 
    1269             :       enddo
    1270             : 
    1271        2864 :       allocate( lm_arr(niter), ic_arr(niter))
    1272         716 :       niter = 0
    1273        2082 :       do ic = 1, atoms%nat
    1274        1366 :          itype = atoms%itype(ic)
    1275       36232 :          do lm = 1,(lcutm1(itype)+1)**2
    1276       34150 :             niter = niter + 1
    1277       34150 :             lm_arr(niter)    = lm
    1278       35516 :             ic_arr(niter)    = ic
    1279             :          enddo 
    1280             :       enddo
    1281         716 :    end subroutine collapse_ic_and_lm_loop
    1282             : 
    1283          12 :    subroutine bessel_calculation(fi, fmpi, mpdata, nqnrm, gridf, qnrm, sphbesmoment, olap, integral)
    1284             :       implicit NONE 
    1285             :       type(t_fleurinput), intent(in)    :: fi
    1286             :       type(t_mpi), intent(in)           :: fmpi
    1287             :       type(t_mpdata), intent(in)        :: mpdata
    1288             :       integer, intent(in)               :: nqnrm
    1289             :       real, intent(in)                  :: gridf(:,:), qnrm(:)
    1290             :       real, intent(inout)               :: sphbesmoment(0:,:,:), olap(:,0:,:,:), integral(:,0:,:,:)
    1291             : 
    1292             :       integer :: iqnrm, itype, i, l, n, ng, buf_sz, root, ierr
    1293          32 :       REAL    :: sphbes_var(fi%atoms%jmtd, 0:maxval(fi%hybinp%lcutm1))
    1294          32 :       REAL    :: sphbesmoment1(fi%atoms%jmtd, 0:maxval(fi%hybinp%lcutm1))
    1295          32 :       REAL    :: rarr(0:fi%hybinp%lexp + 1), rarr1(0:maxval(fi%hybinp%lcutm1))
    1296             :       real    :: rdum, qnorm, rdum1
    1297             : 
    1298          12 :       call timestart("Bessel calculation")
    1299             :       
    1300          12 :       do iqnrm = 1+fmpi%irank, nqnrm, fmpi%isize
    1301         204 :          qnorm = qnrm(iqnrm)
    1302             :          !$OMP parallel do default(none) &
    1303             :          !$OMP shared(olap, integral, sphbesmoment, fi,qnorm, iqnrm, mpdata, gridf) &
    1304         204 :          !$OMP private(itype, rdum, sphbes_var, sphbesmoment1, ng, rarr, rarr1, rdum1, i, l)
    1305             :          DO itype = 1, fi%atoms%ntype
    1306             :             ng = fi%atoms%jri(itype)
    1307             :             rdum = fi%atoms%rmt(itype)
    1308             :             sphbes_var = 0
    1309             :             sphbesmoment1 = 0
    1310             :             IF (abs(qnorm) < 1e-12) THEN
    1311             :                sphbesmoment(0, itype, iqnrm) = rdum**3/3
    1312             :                DO i = 1, ng
    1313             :                   sphbes_var(i, 0) = 1
    1314             :                   sphbesmoment1(i, 0) = fi%atoms%rmsh(i, itype)**2/3 + (rdum**2 - fi%atoms%rmsh(i, itype)**2)/2
    1315             :                END DO
    1316             :             ELSE
    1317             :                call sphbes(fi%hybinp%lexp + 1, qnorm*rdum, rarr)
    1318             :                DO l = 0, fi%hybinp%lexp
    1319             :                   sphbesmoment(l, itype, iqnrm) = rdum**(l + 2)*rarr(l + 1)/qnorm
    1320             :                END DO
    1321             :                DO i = ng, 1, -1
    1322             :                   rdum = fi%atoms%rmsh(i, itype)
    1323             :                   call sphbes(fi%hybinp%lcutm1(itype) + 1, qnorm*rdum, rarr)
    1324             :                   DO l = 0, fi%hybinp%lcutm1(itype)
    1325             :                      sphbes_var(i, l) = rarr(l)
    1326             :                      IF (l /= 0) THEN; rdum1 = -rdum**(1 - l)*rarr(l - 1)
    1327             :                      ELSE; rdum1 = -COS(qnorm*rdum)/qnorm
    1328             :                      ENDIF
    1329             :                      IF (i == ng) rarr1(l) = rdum1
    1330             :                      sphbesmoment1(i, l) = (rdum**(l + 2)*rarr(l + 1)/rdum**(l + 1) &
    1331             :                                             + (rarr1(l) - rdum1)*rdum**l)/qnorm
    1332             :                   END DO
    1333             :                END DO
    1334             :             END IF
    1335             :             DO l = 0, fi%hybinp%lcutm1(itype)
    1336             :                DO n = 1, mpdata%num_radbasfn(l, itype)
    1337             :                   ! note that mpdata%radbasfn_mt already contains one factor rgrid
    1338             :                   olap(n, l, itype, iqnrm) = &
    1339             :                      intgrf(fi%atoms%rmsh(:, itype)*mpdata%radbasfn_mt(:, n, l, itype)*sphbes_var(:, l), &
    1340             :                             fi%atoms, itype, gridf)
    1341             : 
    1342             :                   integral(n, l, itype, iqnrm) = &
    1343             :                      intgrf(fi%atoms%rmsh(:, itype)*mpdata%radbasfn_mt(:, n, l, itype)*sphbesmoment1(:, l), &
    1344             :                             fi%atoms, itype, gridf)
    1345             : 
    1346             :                END DO
    1347             :             END DO
    1348             :          END DO
    1349             :       END DO
    1350             : 
    1351             : #ifdef CPP_MPI 
    1352          12 :       call timestart("bcast bessel")
    1353         420 :       do iqnrm = 1, nqnrm
    1354         408 :          root = mod(iqnrm - 1,fmpi%isize)
    1355         408 :          buf_sz = size(olap,1) * size(olap,2) * size(olap,3)
    1356         408 :          call MPI_Bcast(olap(1,0,1,iqnrm), buf_sz, MPI_DOUBLE_PRECISION, root, fmpi%mpi_comm, ierr)
    1357             : 
    1358         408 :          buf_sz = size(integral,1) * size(integral,2) * size(integral,3)
    1359         408 :          call MPI_Bcast(integral(1,0,1,iqnrm), buf_sz, MPI_DOUBLE_PRECISION, root, fmpi%mpi_comm, ierr)
    1360             : 
    1361         408 :          buf_sz = size(sphbesmoment,1) * size(sphbesmoment,2)
    1362         420 :          call MPI_Bcast(sphbesmoment(0,1,iqnrm), buf_sz, MPI_DOUBLE_PRECISION, root, fmpi%mpi_comm, ierr)
    1363             :       enddo
    1364          12 :       call timestop("bcast bessel")
    1365             : #endif
    1366          12 :       call timestop("Bessel calculation")
    1367          12 :    end subroutine bessel_calculation
    1368             : 
    1369          12 :    function calc_num_mtmts(fi) result(num_mtmt)
    1370             :       implicit none 
    1371             :       type(t_fleurinput), intent(in) :: fi
    1372             :       integer                        :: num_mtmt 
    1373             : 
    1374             :       integer :: itype, ineq, l, m
    1375             : 
    1376          12 :       num_mtmt = 0
    1377          32 :       DO itype = 1, fi%atoms%ntype
    1378          52 :          DO ineq = 1, fi%atoms%neq(itype)
    1379         140 :             DO l = 0, fi%hybinp%lcutm1(itype)
    1380         620 :                DO M = -l, l
    1381         600 :                   num_mtmt = num_mtmt + 1
    1382             :                enddo
    1383             :             enddo 
    1384             :          enddo 
    1385             :       enddo 
    1386          12 :    end function calc_num_mtmts
    1387             :    
    1388             : END MODULE m_coulombmatrix

Generated by: LCOV version 1.14