LCOV - code coverage report
Current view: top level - hybrid - coulombmatrix.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 1183 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 6 0.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             : 
      36             : CONTAINS
      37             : 
      38           0 :    SUBROUTINE coulombmatrix(mpi, atoms, kpts, cell, sym, hybrid, xcpot, l_restart)
      39             : 
      40             :       USE m_types
      41             :       USE m_juDFT
      42             :       USE m_constants, ONLY: pi_const
      43             :       USE m_olap, ONLY: olap_pw, gptnorm
      44             :       USE m_trafo, ONLY: symmetrize, bramat_trafo
      45             :       USE m_util, ONLY: intgrf, intgrf_init, primitivef
      46             :       USE m_hsefunctional, ONLY: change_coulombmatrix
      47             :       USE m_wrapper
      48             :       USE m_io_hybrid
      49             :       use m_ylm
      50             :       use m_sphbes, only: sphbes
      51             :       IMPLICIT NONE
      52             : 
      53             :       TYPE(t_xcpot_inbuild), INTENT(IN)     :: xcpot
      54             :       TYPE(t_mpi), INTENT(IN)       :: mpi
      55             :       TYPE(t_hybrid), INTENT(INOUT) :: hybrid
      56             :       TYPE(t_sym), INTENT(IN)       :: sym
      57             :       TYPE(t_cell), INTENT(IN)      :: cell
      58             :       TYPE(t_kpts), INTENT(IN)      :: kpts
      59             :       TYPE(t_atoms), INTENT(IN)     :: atoms
      60             : 
      61             :       ! - scalars -
      62             :       LOGICAL, INTENT(IN)    :: l_restart
      63             : 
      64             :       ! - local scalars -
      65             :       INTEGER                    :: inviop
      66             :       INTEGER                    :: nqnrm, iqnrm, iqnrm1, iqnrm2, iqnrmstart, iqnrmstep
      67             :       INTEGER                    :: itype, l, ix, iy, iy0, i, j, lm, l1, l2, m1, m2, ineq, idum, ikpt, ikpt0, ikpt1
      68             :       INTEGER                    :: lm1, lm2, itype1, itype2, ineq1, ineq2, n, n1, n2, ng
      69             :       INTEGER                    :: ic, ic1, ic2, ic3, ic4, ic5, ic6, ic7, ic8
      70             :       INTEGER                    :: igpt, igpt1, igpt2, igptp, igptp1, igptp2
      71             :       INTEGER                    :: isym, isym1, isym2, igpt0
      72             :       INTEGER                    :: ok
      73             :       INTEGER                    :: m
      74             :       INTEGER                    :: ikptmin, ikptmax, nkminmax
      75             :       INTEGER                    :: maxfac
      76             : 
      77             :       LOGICAL                    :: lsym
      78             : 
      79             :       REAL                       :: rdum, rdum1, rdum2
      80             :       REAL                       :: svol, qnorm, qnorm1, qnorm2, gnorm
      81             :       REAL                       :: fcoulfac
      82             :       REAL                       :: time1, time2
      83             : 
      84             :       COMPLEX                    :: cdum, cdum1, cexp, csum
      85             : 
      86             :       ! - local arrays -
      87             :       INTEGER                    :: g(3)
      88           0 :       INTEGER                    :: nbasm1(kpts%nkptf)
      89           0 :       INTEGER, ALLOCATABLE   :: pqnrm(:, :)
      90           0 :       INTEGER                    :: rrot(3, 3, sym%nsym), invrrot(3, 3, sym%nsym)
      91           0 :       INTEGER, ALLOCATABLE   :: iarr(:), POINTER(:, :, :, :)!,pointer(:,:,:)
      92           0 :       INTEGER                    :: igptmin(kpts%nkpt), igptmax(kpts%nkpt)
      93           0 :       INTEGER, ALLOCATABLE   :: nsym_gpt(:, :), sym_gpt(:, :, :)
      94           0 :       INTEGER                    :: nsym1(kpts%nkpt + 1), sym1(sym%nsym, kpts%nkpt + 1)
      95             : 
      96           0 :       LOGICAL                    :: calc_mt(kpts%nkpt)
      97             : 
      98             :       REAL                       :: q(3), q1(3), q2(3)
      99           0 :       REAL                       :: integrand(atoms%jmtd), primf1(atoms%jmtd), primf2(atoms%jmtd)
     100           0 :       REAL                       :: mat(hybrid%maxindxm1*(hybrid%maxindxm1 + 1)/2)
     101           0 :       REAL                       :: moment(hybrid%maxindxm1, 0:hybrid%maxlcutm1, atoms%ntype), &
     102           0 :                                     moment2(hybrid%maxindxm1, atoms%ntype)
     103           0 :       REAL                       :: sphbes_var(atoms%jmtd, 0:hybrid%maxlcutm1)
     104           0 :       REAL                       :: sphbesmoment1(atoms%jmtd, 0:hybrid%maxlcutm1)
     105           0 :       REAL                       :: rarr(0:hybrid%lexp + 1), rarr1(0:hybrid%maxlcutm1)
     106           0 :       REAL, ALLOCATABLE   :: gmat(:, :), qnrm(:)
     107           0 :       REAL, ALLOCATABLE   :: sphbesmoment(:, :, :)
     108           0 :       REAL, ALLOCATABLE   :: sphbes0(:, :, :)
     109           0 :       REAL, ALLOCATABLE   :: olap(:, :, :, :), integral(:, :, :, :)
     110           0 :       REAL, ALLOCATABLE   :: gridf(:, :)
     111           0 :       REAL                       :: facA(0:MAX(2*atoms%lmaxd + hybrid%maxlcutm1 + 1, 4*MAX(hybrid%maxlcutm1, hybrid%lexp) + 1))
     112           0 :       REAL                       :: facB(0:MAX(2*atoms%lmaxd + hybrid%maxlcutm1 + 1, 4*MAX(hybrid%maxlcutm1, hybrid%lexp) + 1))
     113           0 :       REAL                       :: facC(-1:MAX(2*atoms%lmaxd + hybrid%maxlcutm1 + 1, 4*MAX(hybrid%maxlcutm1, hybrid%lexp) + 1))
     114             : 
     115           0 :       COMPLEX     :: cexp1(atoms%ntype), csumf(9)
     116           0 :       COMPLEX     :: structconst((2*hybrid%lexp + 1)**2, atoms%nat, atoms%nat, kpts%nkpt)             ! nw = 1
     117           0 :       COMPLEX     :: y((hybrid%lexp + 1)**2), y1((hybrid%lexp + 1)**2), y2((hybrid%lexp + 1)**2)
     118           0 :       COMPLEX     :: dwgn(-hybrid%maxlcutm1:hybrid%maxlcutm1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1, sym%nsym)
     119           0 :       COMPLEX, ALLOCATABLE   :: smat(:, :)
     120           0 :       COMPLEX, ALLOCATABLE   :: coulmat(:, :)
     121           0 :       COMPLEX, ALLOCATABLE   :: carr2(:, :), carr2a(:, :), carr2b(:, :)
     122           0 :       COMPLEX, ALLOCATABLE   :: structconst1(:, :)
     123           0 :       REAL, ALLOCATABLE   :: coulomb_mt1(:, :, :, :, :)
     124             : 
     125             :       !REAL       , ALLOCATABLE   :: coulomb(:,:) !At the moment we always calculate a complex coulomb matrix
     126           0 :       REAL, ALLOCATABLE   :: coulomb_mt2_r(:, :, :, :, :), coulomb_mt3_r(:, :, :, :)
     127           0 :       REAL, ALLOCATABLE   :: coulomb_mtir_r(:, :, :), coulombp_mtir_r(:, :)
     128           0 :       COMPLEX, ALLOCATABLE   :: coulomb(:, :)
     129           0 :       COMPLEX, ALLOCATABLE   :: coulomb_mt2_c(:, :, :, :, :), coulomb_mt3_c(:, :, :, :)
     130           0 :       COMPLEX, ALLOCATABLE   :: coulomb_mtir_c(:, :, :), coulombp_mtir_c(:, :)
     131             : 
     132             :       INTEGER                    :: ishift, ishift1
     133             :       INTEGER                    :: iatom, iatom1
     134             :       INTEGER                    :: indx1, indx2, indx3, indx4
     135             :       LOGICAL                    :: l_found, l_warn, l_warned, l_plot = .FALSE.!.true.!.false.
     136           0 :       TYPE(t_mat)                :: olapm, coulhlp
     137             : 
     138           0 :       CALL timestart("Coulomb matrix setup")
     139             : 
     140           0 :       svol = SQRT(cell%vol)
     141           0 :       fcoulfac = 4*pi_const/cell%vol
     142           0 :       maxfac = MAX(2*atoms%lmaxd + hybrid%maxlcutm1 + 1, 4*MAX(hybrid%maxlcutm1, hybrid%lexp) + 1)
     143             : 
     144           0 :       facA(0) = 1                    !
     145           0 :       facB(0) = 1                    ! Define:
     146           0 :       facC(-1:0) = 1                    ! facA(i)    = i!
     147           0 :       DO i = 1, maxfac                       ! facB(i)   = sqrt(i!)
     148           0 :          facA(i) = facA(i - 1)*i            ! facC(i) = (2i+1)!!
     149           0 :          facB(i) = facB(i - 1)*SQRT(i*1.0) !
     150           0 :          facC(i) = facC(i - 1)*(2*i + 1)   !
     151             :       END DO
     152             : 
     153           0 :       CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, gridf)
     154             : 
     155           0 :       nbasm1 = hybrid%nbasp + hybrid%ngptm(:)
     156             : 
     157             :       !     Calculate the structure constant
     158           0 :       CALL structureconstant(structconst, cell, hybrid, atoms, kpts, mpi)
     159             : 
     160           0 :       IF (mpi%irank == 0) WRITE (6, '(//A)') '### subroutine: coulombmatrix ###'
     161             : 
     162             :       !
     163             :       !     Matrix allocation
     164             :       !
     165             : 
     166           0 :       call timestart("coulomb allocation")
     167             :       IF (ALLOCATED(coulomb)) DEALLOCATE (coulomb)
     168             : 
     169           0 :       ALLOCATE (coulomb(hybrid%maxbasm1*(hybrid%maxbasm1 + 1)/2, kpts%nkpt), stat=ok)
     170           0 :       IF (ok /= 0) STOP 'coulombmatrix: failure allocation coulomb matrix'
     171           0 :       coulomb = 0
     172           0 :       call timestop("coulomb allocation")
     173             : 
     174           0 :       IF (mpi%irank == 0) WRITE (6, '(/A,F6.1," MB")') 'Size of coulomb matrix:', 16.0/1048576*SIZE(coulomb)
     175             : 
     176             :       !     Generate Symmetry:
     177             :       !     Reduce list of g-Points so that only one of each symm-equivalent is calculated
     178             : 
     179             : #ifndef CPP_NOCOULSYM
     180           0 :       IF (mpi%irank == 0) WRITE (6, '(/A)', advance='no') 'Setup for symmetry...'
     181           0 :       CALL cpu_TIME(time1)
     182             :       ! calculate rotations in reciprocal space
     183           0 :       DO isym = 1, sym%nsym
     184           0 :          IF (isym <= sym%nop) THEN
     185           0 :             inviop = sym%invtab(isym)
     186           0 :             rrot(:, :, isym) = TRANSPOSE(sym%mrot(:, :, inviop))
     187           0 :             DO l = 0, hybrid%maxlcutm1
     188             :                dwgn(:, :, l, isym) = TRANSPOSE(hybrid%d_wgn2(-hybrid%maxlcutm1:hybrid%maxlcutm1, &
     189           0 :                                                              -hybrid%maxlcutm1:hybrid%maxlcutm1, l, isym))
     190             :             END DO
     191             :          ELSE
     192           0 :             inviop = isym - sym%nop
     193           0 :             rrot(:, :, isym) = -rrot(:, :, inviop)
     194           0 :             dwgn(:, :, :, isym) = dwgn(:, :, :, inviop)
     195           0 :             DO l = 0, hybrid%maxlcutm1
     196           0 :                DO m1 = -l, l
     197           0 :                   DO m2 = -l, -1
     198           0 :                      cdum = dwgn(m1, m2, l, isym)
     199           0 :                      dwgn(m1, m2, l, isym) = dwgn(m1, -m2, l, isym)*(-1)**m2
     200           0 :                      dwgn(m1, -m2, l, isym) = cdum*(-1)**m2
     201             :                   END DO
     202             :                END DO
     203             :             END DO
     204             :          END IF
     205             :       END DO
     206           0 :       invrrot(:, :, :sym%nop) = rrot(:, :, sym%invtab)
     207           0 :       IF (sym%nsym > sym%nop) THEN
     208           0 :          invrrot(:, :, sym%nop + 1:) = rrot(:, :, sym%invtab + sym%nop)
     209             :       END IF
     210             : 
     211             :       ! Get symmetry operations that leave bk(:,ikpt) invariant -> sym1
     212           0 :       nsym1 = 0
     213           0 :       DO ikpt = 1, kpts%nkpt
     214           0 :          isym1 = 0
     215           0 :          DO isym = 1, sym%nsym
     216             :             ! temporary fix until bramat_trafo is correct
     217             :             ! for systems with symmetries including translations
     218           0 :             IF (isym > sym%nop) THEN
     219           0 :                isym2 = isym - sym%nop
     220             :             ELSE
     221             :                isym2 = isym
     222             :             END IF
     223           0 :             IF (ANY(sym%tau(:, isym2) /= 0)) CYCLE
     224             : 
     225           0 :             IF (ALL(ABS(MATMUL(rrot(:, :, isym), kpts%bk(:, ikpt)) - kpts%bk(:, ikpt)) < 1e-12)) THEN
     226           0 :                isym1 = isym1 + 1
     227           0 :                sym1(isym1, ikpt) = isym
     228             :             END IF
     229             :          END DO
     230           0 :          nsym1(ikpt) = isym1
     231             :       END DO
     232             :       ! Define reduced lists of G points -> pgptm1(:,ikpt), ikpt=1,..,nkpt
     233             :       !ALLOCATE ( hybrid%pgptm1(hybrid%maxgptm,kpts%nkpt)) !in mixedbasis
     234             :       ALLOCATE (iarr(hybrid%maxgptm), POINTER(kpts%nkpt, &
     235             :                                               MINVAL(hybrid%gptm(1, :)) - 1:MAXVAL(hybrid%gptm(1, :)) + 1, &
     236             :                                               MINVAL(hybrid%gptm(2, :)) - 1:MAXVAL(hybrid%gptm(2, :)) + 1, &
     237           0 :                                               MINVAL(hybrid%gptm(3, :)) - 1:MAXVAL(hybrid%gptm(3, :)) + 1))
     238           0 :       hybrid%pgptm1 = 0; iarr = 0; POINTER = 0
     239           0 :       DO ikpt = 1, kpts%nkpt
     240           0 :          DO igpt = 1, hybrid%ngptm(ikpt)
     241           0 :             g = hybrid%gptm(:, hybrid%pgptm(igpt, ikpt))
     242           0 :             POINTER(ikpt, g(1), g(2), g(3)) = igpt
     243             :          END DO
     244           0 :          iarr = 0
     245           0 :          j = 0
     246           0 :          DO igpt = hybrid%ngptm(ikpt), 1, -1
     247           0 :             IF (iarr(igpt) == 0) THEN
     248           0 :                j = j + 1
     249           0 :                hybrid%pgptm1(j, ikpt) = igpt
     250           0 :                DO isym1 = 1, nsym1(ikpt)
     251           0 :                   g = MATMUL(rrot(:, :, sym1(isym1, ikpt)), hybrid%gptm(:, hybrid%pgptm(igpt, ikpt)))
     252           0 :                   i = POINTER(ikpt, g(1), g(2), g(3))
     253           0 :                   IF (i == 0) STOP 'coulombmatrix: zero pointer (bug?)'
     254           0 :                   iarr(i) = 1
     255             :                END DO
     256             :             END IF
     257             :          END DO
     258           0 :          hybrid%ngptm1(ikpt) = j
     259             :       END DO
     260           0 :       DEALLOCATE (iarr)
     261             : 
     262           0 :       IF (mpi%irank == 0) WRITE (6, '(12X,A)', advance='no') 'done'
     263           0 :       CALL cpu_TIME(time2)
     264           0 :       IF (mpi%irank == 0) WRITE (6, '(2X,A,F8.2,A)') '( Timing:', time2 - time1, ' )'
     265             : 
     266             :       ! no symmetry used
     267             : #else
     268             :       ALLOCATE (hybrid%pgptm1(hybrid%maxgptm, kpts%nkpt))
     269             :       DO ikpt = 1, kpts%nkpt
     270             :          hybrid%pgptm1(:, ikpt) = (/(igpt0, igpt0=1, hybrid%maxgptm)/)
     271             :          hybrid%ngptm1(ikpt) = hybrid%ngptm(ikpt)
     272             :       END DO
     273             : #endif
     274             : 
     275             :       ! Distribute the work as equally as possible over the processes
     276             :       ikptmin = 1
     277             :       ikptmax = kpts%nkpt
     278           0 :       igptmin = 1
     279           0 :       igptmax = hybrid%ngptm1(:kpts%nkpt)
     280           0 :       calc_mt = .TRUE.
     281           0 :       nkminmax = kpts%nkpt
     282             : 
     283           0 :       IF (mpi%irank == 0) WRITE (6, '(A)', advance='no') 'Preparations...'
     284           0 :       CALL cpu_TIME(time1)
     285             : 
     286           0 :       call timestart("define gmat")
     287             :       ! Define gmat (symmetric)
     288           0 :       i = (hybrid%lexp + 1)**2
     289           0 :       ALLOCATE (gmat(i, i))
     290           0 :       gmat = 0
     291             :       lm1 = 0
     292           0 :       DO l1 = 0, hybrid%lexp
     293           0 :          DO m1 = -l1, l1
     294           0 :             lm1 = lm1 + 1
     295           0 :             lm2 = 0
     296           0 :             lp1: DO l2 = 0, l1
     297           0 :                DO m2 = -l2, l2
     298           0 :                   lm2 = lm2 + 1
     299           0 :                   IF (lm2 > lm1) EXIT lp1 ! Don't cross the diagonal!
     300             :                   gmat(lm1, lm2) = facB(l1 + l2 + m2 - m1)*facB(l1 + l2 + m1 - m2)/ &
     301             :                                    (facB(l1 + m1)*facB(l1 - m1)*facB(l2 + m2)*facB(l2 - m2))/ &
     302           0 :                                    SQRT(1.0*(2*l1 + 1)*(2*l2 + 1)*(2*(l1 + l2) + 1))*(4*pi_const)**1.5
     303           0 :                   gmat(lm2, lm1) = gmat(lm1, lm2)
     304             :                END DO
     305             :             END DO LP1
     306             :          END DO
     307             :       END DO
     308           0 :       call timestop("define gmat")
     309             : 
     310             :       ! Calculate moments of MT functions
     311           0 :       call timestart("calc moments of MT")
     312           0 :       DO itype = 1, atoms%ntype
     313           0 :          DO l = 0, hybrid%lcutm1(itype)
     314           0 :             DO i = 1, hybrid%nindxm1(l, itype)
     315             :                ! note that hybrid%basm1 already contains the factor rgrid
     316             :                moment(i, l, itype) = intgrf(atoms%rmsh(:, itype)**(l + 1)*hybrid%basm1(:, i, l, itype), &
     317           0 :                                             atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     318             :             END DO
     319             :          END DO
     320           0 :          DO i = 1, hybrid%nindxm1(0, itype)
     321             :             moment2(i, itype) = intgrf(atoms%rmsh(:, itype)**3*hybrid%basm1(:, i, 0, itype), &
     322           0 :                                        atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     323             :          END DO
     324             :       END DO
     325           0 :       call timestop("calc moments of MT")
     326             : 
     327           0 :       call timestart("getnorm")
     328             :       ! Look for different qnorm = |k+G|, definition of qnrm and pqnrm.
     329           0 :       CALL getnorm(kpts, hybrid%gptm, hybrid%ngptm, hybrid%pgptm, qnrm, nqnrm, pqnrm, cell)
     330             :       ALLOCATE (sphbesmoment(0:hybrid%lexp, atoms%ntype, nqnrm), &
     331             :                 olap(hybrid%maxindxm1, 0:hybrid%maxlcutm1, atoms%ntype, nqnrm), &
     332           0 :                 integral(hybrid%maxindxm1, 0:hybrid%maxlcutm1, atoms%ntype, nqnrm))
     333           0 :       sphbes_var = 0
     334           0 :       sphbesmoment = 0
     335           0 :       sphbesmoment1 = 0
     336           0 :       olap = 0
     337           0 :       integral = 0
     338             : 
     339             :       ! Calculate moments of spherical Bessel functions (for (2) and (3))              (->sphbesmoment)
     340             :       ! Calculate overlap of spherical Bessel functions with basis functions (for (2)) (->olap)
     341             :       ! Calculate overlap of sphbesmoment1(r,l)         with basis functions (for (2)) (->integral)
     342             :       ! We use           sphbes(r,l) = j_l(qr)
     343             :       ! and       sphbesmoment1(r,l) = 1/r**(l-1) * INT(0..r) r'**(l+2) * j_l(qr') dr'
     344             :       !                                + r**(l+2) * INT(r..S) r'**(1-l) * j_l(qr') dr' .
     345             : 
     346           0 :       iqnrmstart = mpi%irank + 1
     347           0 :       iqnrmstep = mpi%isize
     348           0 :       call timestop("getnorm")
     349             : 
     350           0 :       call timestart("Bessel calculation")
     351           0 :       DO iqnrm = iqnrmstart, nqnrm, iqnrmstep
     352           0 :          qnorm = qnrm(iqnrm)
     353           0 :          DO itype = 1, atoms%ntype
     354           0 :             ng = atoms%jri(itype)
     355           0 :             rdum = atoms%rmt(itype)
     356           0 :             sphbes_var = 0
     357           0 :             sphbesmoment1 = 0
     358           0 :             IF (qnorm == 0) THEN
     359           0 :                sphbesmoment(0, itype, iqnrm) = rdum**3/3
     360           0 :                DO i = 1, ng
     361           0 :                   sphbes_var(i, 0) = 1
     362           0 :                   sphbesmoment1(i, 0) = atoms%rmsh(i, itype)**2/3 + (rdum**2 - atoms%rmsh(i, itype)**2)/2
     363             :                END DO
     364             :             ELSE
     365           0 :                call sphbes(hybrid%lexp + 1, qnorm*rdum, rarr)
     366           0 :                DO l = 0, hybrid%lexp
     367           0 :                   sphbesmoment(l, itype, iqnrm) = rdum**(l + 2)*rarr(l + 1)/qnorm
     368             :                END DO
     369           0 :                DO i = ng, 1, -1
     370           0 :                   rdum = atoms%rmsh(i, itype)
     371           0 :                   call sphbes(hybrid%lcutm1(itype) + 1, qnorm*rdum, rarr)
     372           0 :                   DO l = 0, hybrid%lcutm1(itype)
     373           0 :                      sphbes_var(i, l) = rarr(l)
     374           0 :                      IF (l /= 0) THEN; rdum1 = -rdum**(1 - l)*rarr(l - 1)
     375           0 :                      ELSE; rdum1 = -COS(qnorm*rdum)/qnorm
     376             :                      ENDIF
     377           0 :                      IF (i == ng) rarr1(l) = rdum1
     378             :                      sphbesmoment1(i, l) = (rdum**(l + 2)*rarr(l + 1)/rdum**(l + 1) &
     379           0 :                                             + (rarr1(l) - rdum1)*rdum**l)/qnorm
     380             :                   END DO
     381             :                END DO
     382             :             END IF
     383           0 :             DO l = 0, hybrid%lcutm1(itype)
     384           0 :                DO n = 1, hybrid%nindxm1(l, itype)
     385             :                   ! note that hybrid%basm1 already contains one factor rgrid
     386             :                   olap(n, l, itype, iqnrm) = &
     387             :                      intgrf(atoms%rmsh(:, itype)*hybrid%basm1(:, n, l, itype)*sphbes_var(:, l), &
     388           0 :                             atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     389             : 
     390             :                   integral(n, l, itype, iqnrm) = &
     391             :                      intgrf(atoms%rmsh(:, itype)*hybrid%basm1(:, n, l, itype)*sphbesmoment1(:, l), &
     392           0 :                             atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     393             : 
     394             :                END DO
     395             :             END DO
     396             :          END DO
     397             :       END DO
     398           0 :       call timestop("Bessel calculation")
     399             : 
     400           0 :       IF (mpi%irank == 0) THEN
     401           0 :          WRITE (6, '(18X,A)', advance='no') 'done'
     402           0 :          CALL cpu_TIME(time2)
     403           0 :          WRITE (6, '(2X,A,F8.2,A)', advance='no') '( Timing:', time2 - time1, ' )'
     404           0 :          WRITE (6, *)
     405             :       END IF
     406             : 
     407             :       !
     408             :       !     (1) Case < MT | v | MT >
     409             :       !
     410             : 
     411           0 :       IF (mpi%irank == 0) WRITE (6, '(A)', advance='no') '< MT | v | MT > contribution...'
     412             : 
     413           0 :       CALL cpu_TIME(time1)
     414             : 
     415           0 :       IF (ANY(calc_mt)) THEN
     416             : 
     417             :          !       (1a) r,r' in same MT
     418           0 :          call timestart("loop 1")
     419           0 :          ix = 0
     420           0 :          iy = 0
     421           0 :          iy0 = 0
     422           0 :          DO itype = 1, atoms%ntype
     423           0 :             DO ineq = 1, atoms%neq(itype)
     424             :                ! Here the diagonal block matrices do not depend on ineq. In (1b) they do depend on ineq, though,
     425           0 :                DO l = 0, hybrid%lcutm1(itype)
     426           0 :                   DO n2 = 1, hybrid%nindxm1(l, itype)
     427             :                      ! note that hybrid%basm1 already contains the factor rgrid
     428             :                      CALL primitivef(primf1, hybrid%basm1(:, n2, l, itype) &
     429             :                                      *atoms%rmsh(:, itype)**(l + 1), atoms%rmsh, atoms%dx, &
     430           0 :                                      atoms%jri, atoms%jmtd, itype, atoms%ntype)
     431             :                      ! -itype is to enforce inward integration
     432             :                      CALL primitivef(primf2, hybrid%basm1(:atoms%jri(itype), n2, l, itype) &
     433             :                                      /atoms%rmsh(:atoms%jri(itype), itype)**l, atoms%rmsh, atoms%dx, &
     434           0 :                                      atoms%jri, atoms%jmtd, -itype, atoms%ntype)
     435             : 
     436           0 :                      primf1(:atoms%jri(itype)) = primf1(:atoms%jri(itype))/atoms%rmsh(:atoms%jri(itype), itype)**l
     437           0 :                      primf2 = primf2*atoms%rmsh(:, itype)**(l + 1)
     438             : 
     439           0 :                      DO n1 = 1, n2
     440           0 :                         integrand = hybrid%basm1(:, n1, l, itype)*(primf1 + primf2)
     441             :                         !                 call intgr0( (4*pimach())/(2*l+1)*integrand,rmsh(1,itype),dx(itype),jri(itype),mat(n2*(n2-1)/2+n1) )
     442             :                         mat(n2*(n2 - 1)/2 + n1) = (4*pi_const)/(2*l + 1) &
     443             :                                                   *intgrf(integrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, &
     444           0 :                                                           atoms%ntype, itype, gridf)
     445             :                      END DO
     446             :                   END DO
     447             : 
     448             :                   ! distribute mat for m=-l,l on coulomb in block-matrix form
     449           0 :                   DO M = -l, l
     450           0 :                      DO n2 = 1, hybrid%nindxm1(l, itype)
     451           0 :                         ix = ix + 1
     452           0 :                         iy = iy0
     453           0 :                         DO n1 = 1, n2
     454           0 :                            iy = iy + 1
     455           0 :                            i = ix*(ix - 1)/2 + iy
     456           0 :                            j = n2*(n2 - 1)/2 + n1
     457           0 :                            coulomb(i, kpts%nkpt) = mat(j)
     458             :                         END DO
     459             :                      END DO
     460           0 :                      iy0 = ix
     461             :                   END DO
     462             : 
     463             :                END DO
     464             :             END DO
     465             :          END DO
     466           0 :          call timestop("loop 1")
     467             : 
     468             :          !       (1b) r,r' in different MT
     469             : 
     470           0 :          ALLOCATE (coulmat(hybrid%nbasp, hybrid%nbasp), stat=ok)
     471           0 :          IF (ok /= 0) STOP 'coulombmatrix: failure allocation coulmat'
     472           0 :          coulmat = 0
     473             : 
     474             :       END IF
     475             : 
     476           0 :       DO ikpt = ikptmin, ikptmax
     477             : 
     478             :          ! only the first rank handles the MT-MT part
     479           0 :          call timestart("MT-MT part")
     480           0 :          IF (calc_mt(ikpt)) THEN
     481             : 
     482           0 :             ix = 0
     483           0 :             ic2 = 0
     484           0 :             DO itype2 = 1, atoms%ntype
     485           0 :                DO ineq2 = 1, atoms%neq(itype2)
     486           0 :                   ic2 = ic2 + 1
     487           0 :                   lm2 = 0
     488           0 :                   DO l2 = 0, hybrid%lcutm1(itype2)
     489           0 :                      DO m2 = -l2, l2
     490           0 :                         lm2 = lm2 + 1
     491           0 :                         DO n2 = 1, hybrid%nindxm1(l2, itype2)
     492           0 :                            ix = ix + 1
     493             : 
     494           0 :                            iy = 0
     495           0 :                            ic1 = 0
     496           0 :                            lp2: DO itype1 = 1, itype2
     497           0 :                               DO ineq1 = 1, atoms%neq(itype1)
     498           0 :                                  ic1 = ic1 + 1
     499           0 :                                  lm1 = 0
     500           0 :                                  DO l1 = 0, hybrid%lcutm1(itype1)
     501           0 :                                     DO m1 = -l1, l1
     502           0 :                                        lm1 = lm1 + 1
     503           0 :                                        DO n1 = 1, hybrid%nindxm1(l1, itype1)
     504           0 :                                           iy = iy + 1
     505           0 :                                           IF (iy > ix) EXIT lp2 ! Don't cross the diagonal!
     506           0 :                                           rdum = (-1)**(l2 + m2)*moment(n1, l1, itype1)*moment(n2, l2, itype2)*gmat(lm1, lm2)
     507           0 :                                           l = l1 + l2
     508           0 :                                           lm = l**2 + l + m1 - m2 + 1
     509           0 :                                           idum = ix*(ix - 1)/2 + iy
     510             :                                           coulmat(iy, ix) = coulomb(idum, kpts%nkpt) &
     511             :                                                             + EXP(CMPLX(0.0, 1.0)*2*pi_const* &
     512             :                                                                   dot_PRODUCT(kpts%bk(:, ikpt), &
     513             :                                                                               atoms%taual(:, ic2) - atoms%taual(:, ic1))) &
     514           0 :                                                             *rdum*structconst(lm, ic1, ic2, ikpt)
     515           0 :                                           coulmat(ix, iy) = CONJG(coulmat(iy, ix))
     516             :                                        END DO
     517             :                                     END DO
     518             :                                  END DO
     519             :                               END DO
     520             :                            END DO lp2
     521             : 
     522             :                         END DO
     523             :                      END DO
     524             :                   END DO
     525             :                END DO
     526             :             END DO
     527             : 
     528           0 :             IF (sym%invs) THEN
     529             :                !symmetrize makes the Coulomb matrix real symmetric
     530             :                CALL symmetrize(coulmat, hybrid%nbasp, hybrid%nbasp, 3, .FALSE., &
     531             :                                atoms, hybrid%lcutm1, hybrid%maxlcutm1, &
     532           0 :                                hybrid%nindxm1, sym)
     533             :             ENDIF
     534             : 
     535           0 :             coulomb(:hybrid%nbasp*(hybrid%nbasp + 1)/2, ikpt) = packmat(coulmat)
     536             : 
     537             :          END IF
     538           0 :          call timestop("MT-MT part")
     539             : 
     540             :       END DO
     541           0 :       IF (ANY(calc_mt)) DEALLOCATE (coulmat)
     542             : 
     543           0 :       IF (mpi%irank == 0) THEN
     544           0 :          WRITE (6, '(2X,A)', advance='no') 'done'
     545           0 :          CALL cpu_TIME(time2)
     546           0 :          WRITE (6, '(2X,A,F8.2,A)', advance='no') '( Timing:', time2 - time1, ' )'
     547           0 :          WRITE (6, *)
     548             :       END IF
     549             : 
     550           0 :       IF (hybrid%maxgptm == 0) GOTO 1 ! skip calculation of plane-wave contribution if mixed basis does not contain plane waves
     551             : 
     552             :       !
     553             :       !     (2) Case < MT | v | PW >
     554             :       !
     555             : 
     556           0 :       IF (mpi%irank == 0) WRITE (6, '(A)', advance='no') '< MT | v | PW > contribution...'
     557             : 
     558           0 :       CALL cpu_TIME(time1)
     559             : 
     560             :       !     (2a) r in MT, r' everywhere
     561             :       !     (2b) r,r' in same MT
     562             :       !     (2c) r,r' in different MT
     563             : 
     564           0 :       ALLOCATE (coulmat(hybrid%nbasp, hybrid%maxgptm), stat=ok)
     565           0 :       IF (ok /= 0) STOP 'coulombmatrix: failure allocation coulmat'
     566           0 :       coulmat = 0
     567             : 
     568           0 :       call timestart("loop over interst.")
     569           0 :       DO ikpt = ikptmin, ikptmax !1,kpts%nkpt
     570             : 
     571           0 :          coulmat = 0
     572             :          ! start to loop over interstitial plane waves
     573           0 :          DO igpt0 = igptmin(ikpt), igptmax(ikpt) !1,hybrid%ngptm1(ikpt)
     574           0 :             igpt = hybrid%pgptm1(igpt0, ikpt)
     575           0 :             igptp = hybrid%pgptm(igpt, ikpt)
     576           0 :             ix = hybrid%nbasp + igpt
     577           0 :             q = MATMUL(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp), cell%bmat)
     578           0 :             qnorm = SQRT(SUM(q**2))
     579           0 :             iqnrm = pqnrm(igpt, ikpt)
     580           0 :             IF (ABS(qnrm(iqnrm) - qnorm) > 1e-12) then
     581           0 :                STOP 'coulombmatrix: qnorm does not equal corresponding & element in qnrm (bug?)' ! We shouldn't stop here!
     582             :             endif
     583             : 
     584           0 :             call timestart("harmonics")
     585           0 :             call ylm4(2, MATMUL(kpts%bk(:, kpts%nkpt), cell%bmat), y1)
     586           0 :             call ylm4(2, MATMUL(hybrid%gptm(:, igptp), cell%bmat), y2)
     587           0 :             call ylm4(hybrid%lexp, q, y)
     588           0 :             call timestop("harmonics")
     589           0 :             y1 = CONJG(y1); y2 = CONJG(y2); y = CONJG(y)
     590             : 
     591           0 :             iy = 0
     592           0 :             ic = 0
     593           0 :             DO itype = 1, atoms%ntype
     594           0 :                DO ineq = 1, atoms%neq(itype)
     595           0 :                   ic = ic + 1
     596           0 :                   lm = 0
     597           0 :                   DO l = 0, hybrid%lcutm1(itype)
     598           0 :                      DO M = -l, l
     599           0 :                         lm = lm + 1
     600             : 
     601             :                         ! calculate sum over lm and centers for (2c) -> csum, csumf
     602           0 :                         csum = 0
     603             :                         csumf = 0
     604           0 :                         ic1 = 0
     605           0 :                         DO itype1 = 1, atoms%ntype
     606           0 :                            DO ineq1 = 1, atoms%neq(itype1)
     607           0 :                               ic1 = ic1 + 1
     608             :                               cexp = 4*pi_const*EXP(CMPLX(0.0, 1.0)*2*pi_const &
     609             :                                                     *(dot_PRODUCT(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp), atoms%taual(:, ic1)) &
     610           0 :                                                       - dot_PRODUCT(kpts%bk(:, ikpt), atoms%taual(:, ic))))
     611             : 
     612           0 :                               lm1 = 0
     613           0 :                               DO l1 = 0, hybrid%lexp
     614           0 :                                  l2 = l + l1 ! for structconst
     615           0 :                                  idum = 1
     616           0 :                                  cdum = sphbesmoment(l1, itype1, iqnrm)*CMPLX(0.0, 1.0)**(l1)*cexp
     617           0 :                                  DO m1 = -l1, l1
     618           0 :                                     lm1 = lm1 + 1
     619           0 :                                     m2 = M - m1              ! for structconst
     620           0 :                                     lm2 = l2**2 + l2 + m2 + 1 !
     621           0 :                                     csum = csum - idum*gmat(lm1, lm)*y(lm1)*cdum*structconst(lm2, ic, ic1, ikpt)
     622           0 :                                     idum = -idum ! factor (-1)*(l1+m1)
     623             :                                  END DO
     624             :                               END DO
     625             : 
     626             :                               ! add contribution of (2c) to csum and csumf coming from linear and quadratic orders of Y_lm*(G) / G * j_(l+1)(GS)
     627           0 :                               IF (ikpt == 1 .AND. l <= 2) THEN
     628             :                                  cexp = EXP(CMPLX(0.0, 1.0)*2*pi_const*dot_PRODUCT(hybrid%gptm(:, igptp), atoms%taual(:, ic1))) &
     629           0 :                                         *gmat(lm, 1)*4*pi_const/cell%vol
     630             :                                  csumf(lm) = csumf(lm) - cexp*SQRT(4*pi_const)* &
     631           0 :                                              CMPLX(0.0, 1.0)**l*sphbesmoment(0, itype1, iqnrm)/facC(l - 1)
     632           0 :                                  IF (l == 0) THEN
     633           0 :                                     IF (igpt /= 1) THEN
     634             :                                        csum = csum - cexp*(sphbesmoment(0, itype1, iqnrm)*atoms%rmt(itype1)**2 - &
     635           0 :                                                            sphbesmoment(2, itype1, iqnrm)*2.0/3)/10
     636             :                                     ELSE
     637           0 :                                        csum = csum - cexp*atoms%rmt(itype1)**5/30
     638             :                                     END IF
     639           0 :                                  ELSE IF (l == 1) THEN
     640             :                                     csum = csum + cexp*CMPLX(0.0, 1.0)*SQRT(4*pi_const) &
     641           0 :                                            *sphbesmoment(1, itype1, iqnrm)*y(lm)/3
     642             :                                  END IF
     643             :                               END IF
     644             : 
     645             :                            END DO
     646             :                         END DO
     647             : 
     648             :                         ! add contribution of (2a) to csumf
     649           0 :                         IF (ikpt == 1 .AND. igpt == 1 .AND. l <= 2) THEN
     650             :                            csumf(lm) = csumf(lm) + (4*pi_const)**2*CMPLX(0.0, 1.0)**l/facC(l)
     651             :                         END IF
     652             : 
     653             :                         ! finally define coulomb
     654           0 :                         idum = ix*(ix - 1)/2
     655             :                         cdum = (4*pi_const)**2*CMPLX(0.0, 1.0)**(l)*y(lm) &
     656             :                                *EXP(CMPLX(0.0, 1.0)*2*pi_const &
     657           0 :                                     *dot_PRODUCT(hybrid%gptm(:, igptp), atoms%taual(:, ic)))
     658           0 :                         DO n = 1, hybrid%nindxm1(l, itype)
     659           0 :                            iy = iy + 1
     660             : 
     661           0 :                            IF (ikpt == 1 .AND. igpt == 1) THEN
     662           0 :                               IF (l == 0) coulmat(iy, ix - hybrid%nbasp) = &
     663           0 :                                  -cdum*moment2(n, itype)/6/svol         ! (2a)
     664             :                               coulmat(iy, ix - hybrid%nbasp) = coulmat(iy, ix - hybrid%nbasp) &
     665             :                                                                + (-cdum/(2*l + 1)*integral(n, l, itype, iqnrm) & ! (2b)&
     666           0 :                                                                   + csum*moment(n, l, itype))/svol          ! (2c)
     667             :                            ELSE
     668             :                               coulmat(iy, ix - hybrid%nbasp) = &
     669             :                                  (cdum*olap(n, l, itype, iqnrm)/qnorm**2 &  ! (2a)&
     670             :                                   - cdum/(2*l + 1)*integral(n, l, itype, iqnrm) & ! (2b)&
     671           0 :                                   + csum*moment(n, l, itype))/svol          ! (2c)
     672             : 
     673             :                            END IF
     674             : 
     675             :                         END DO
     676             : 
     677             :                      END DO
     678             :                   END DO
     679             :                END DO
     680             : 
     681             :             END DO
     682             :          END DO
     683             : 
     684           0 :          IF (sym%invs) THEN
     685             :             CALL symmetrize(coulmat, hybrid%nbasp, hybrid%ngptm(ikpt), 1, .FALSE., &
     686           0 :                             atoms, hybrid%lcutm1, hybrid%maxlcutm1, hybrid%nindxm1, sym)
     687             :          ENDIF
     688             : 
     689           0 :          M = hybrid%nbasp*(hybrid%nbasp + 1)/2
     690           0 :          DO i = 1, hybrid%ngptm(ikpt)
     691           0 :             DO j = 1, hybrid%nbasp + i
     692           0 :                M = M + 1
     693           0 :                IF (j <= hybrid%nbasp) coulomb(M, ikpt) = coulmat(j, i)
     694             :             END DO
     695             :          END DO
     696             :       END DO
     697           0 :       call timestop("loop over interst.")
     698             : 
     699           0 :       DEALLOCATE (coulmat, olap, integral)
     700             : 
     701           0 :       IF (mpi%irank == 0) THEN
     702           0 :          WRITE (6, '(2X,A)', advance='no') 'done'
     703           0 :          CALL cpu_TIME(time2)
     704           0 :          WRITE (6, '(2X,A,F8.2,A)') '( Timing:', time2 - time1, ' )'
     705             :       END IF
     706             : 
     707             :       !
     708             :       !     (3) Case < PW | v | PW >
     709             :       !
     710             : 
     711           0 :       IF (mpi%irank == 0) WRITE (6, '(A)', advance='no') '< PW | v | PW > contribution...'
     712             : 
     713           0 :       CALL cpu_TIME(time1)
     714             : 
     715             :       !     (3a) r,r' everywhere; r everywhere, r' in MT; r in MT, r' everywhere
     716             : 
     717           0 :       CALL cpu_TIME(time1)
     718             :       ! Calculate the hermitian matrix smat(i,j) = sum(a) integral(MT(a)) exp[i(Gj-Gi)r] dr
     719           0 :       call timestart("calc smat")
     720           0 :       ALLOCATE (smat(hybrid%gptmd, hybrid%gptmd))
     721           0 :       smat = 0
     722           0 :       DO igpt2 = 1, hybrid%gptmd
     723           0 :          DO igpt1 = 1, igpt2
     724           0 :             g = hybrid%gptm(:, igpt2) - hybrid%gptm(:, igpt1)
     725           0 :             gnorm = gptnorm(g, cell%bmat)
     726           0 :             IF (gnorm == 0) THEN
     727           0 :                DO itype = 1, atoms%ntype
     728           0 :                   smat(igpt1, igpt2) = smat(igpt1, igpt2) + atoms%neq(itype)*4*pi_const*atoms%rmt(itype)**3/3
     729             :                END DO
     730             :             ELSE
     731           0 :                ic = 0
     732           0 :                DO itype = 1, atoms%ntype
     733           0 :                   rdum = atoms%rmt(itype)*gnorm
     734           0 :                   rdum = 4*pi_const*(SIN(rdum) - rdum*COS(rdum))/gnorm**3
     735           0 :                   DO ineq = 1, atoms%neq(itype)
     736           0 :                      ic = ic + 1
     737             :                      smat(igpt1, igpt2) = smat(igpt1, igpt2) &
     738           0 :                                           + rdum*EXP(CMPLX(0.0, 1.0)*2*pi_const*dot_PRODUCT(atoms%taual(:, ic), g))
     739             :                   END DO
     740             :                END DO
     741             :             END IF
     742           0 :             smat(igpt2, igpt1) = CONJG(smat(igpt1, igpt2))
     743             :          END DO
     744             :       END DO
     745           0 :       call timestop("calc smat")
     746             : 
     747             :       ! Coulomb matrix, contribution (3a)
     748           0 :       call timestart("coulomb matrix")
     749           0 :       DO ikpt = ikptmin, ikptmax
     750             : 
     751           0 :          DO igpt0 = igptmin(ikpt), igptmax(ikpt)
     752           0 :             igpt2 = hybrid%pgptm1(igpt0, ikpt)
     753           0 :             igptp2 = hybrid%pgptm(igpt2, ikpt)
     754           0 :             ix = hybrid%nbasp + igpt2
     755           0 :             iy = hybrid%nbasp
     756           0 :             q2 = MATMUL(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp2), cell%bmat)
     757           0 :             rdum2 = SUM(q2**2)
     758           0 :             IF (rdum2 /= 0) rdum2 = 4*pi_const/rdum2
     759             : 
     760           0 :             DO igpt1 = 1, igpt2
     761           0 :                igptp1 = hybrid%pgptm(igpt1, ikpt)
     762           0 :                iy = iy + 1
     763           0 :                q1 = MATMUL(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp1), cell%bmat)
     764           0 :                idum = ix*(ix - 1)/2 + iy
     765           0 :                rdum1 = SUM(q1**2)
     766           0 :                IF (rdum1 /= 0) rdum1 = 4*pi_const/rdum1
     767             : 
     768           0 :                IF (ikpt == 1) THEN
     769           0 :                   IF (igpt1 /= 1) THEN
     770           0 :                      coulomb(idum, 1) = -smat(igptp1, igptp2)*rdum1/cell%vol
     771             :                   END IF
     772           0 :                   IF (igpt2 /= 1) THEN
     773           0 :                      coulomb(idum, 1) = coulomb(idum, 1) - smat(igptp1, igptp2)*rdum2/cell%vol
     774             :                   END IF
     775             :                ELSE
     776           0 :                   coulomb(idum, ikpt) = -smat(igptp1, igptp2)*(rdum1 + rdum2)/cell%vol
     777             :                END IF
     778             :             END DO
     779           0 :             IF (ikpt /= 1 .OR. igpt2 /= 1) THEN                  !
     780           0 :                coulomb(idum, ikpt) = coulomb(idum, ikpt) + rdum2 ! diagonal term
     781             :             END IF                                            !
     782             :          END DO
     783             : 
     784             :       END DO
     785           0 :       call timestop("coulomb matrix")
     786             :       !     (3b) r,r' in different MT
     787             : 
     788           0 :       call timestart("loop 4:")
     789           0 :       DO ikpt = ikptmin, ikptmax!1,kpts%nkpt
     790             : 
     791             :          ! group together quantities which depend only on l,m and igpt -> carr2a
     792           0 :          ALLOCATE (carr2a((hybrid%lexp + 1)**2, hybrid%maxgptm), carr2b(atoms%nat, hybrid%maxgptm))
     793           0 :          carr2a = 0; carr2b = 0
     794           0 :          DO igpt = 1, hybrid%ngptm(ikpt)
     795           0 :             igptp = hybrid%pgptm(igpt, ikpt)
     796           0 :             iqnrm = pqnrm(igpt, ikpt)
     797           0 :             q = MATMUL(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp), cell%bmat)
     798             : 
     799           0 :             call timestart("harmonics")
     800           0 :             call ylm4(hybrid%lexp, q, y)
     801           0 :             call timestop("harmonics")
     802             : 
     803           0 :             y = CONJG(y)
     804           0 :             lm = 0
     805           0 :             DO l = 0, hybrid%lexp
     806           0 :                DO M = -l, l
     807           0 :                   lm = lm + 1
     808           0 :                   carr2a(lm, igpt) = 4*pi_const*CMPLX(0.0, 1.0)**(l)*y(lm)
     809             :                END DO
     810             :             END DO
     811           0 :             DO ic = 1, atoms%nat
     812             :                carr2b(ic, igpt) = EXP(-CMPLX(0.0, 1.0)*2*pi_const* &
     813           0 :                                       dot_PRODUCT(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp), atoms%taual(:, ic)))
     814             :             END DO
     815             :          END DO
     816             : 
     817             :          !finally we can loop over the plane waves (G: igpt1,igpt2)
     818           0 :          call timestart("loop over plane waves")
     819             :          ALLOCATE (carr2(atoms%nat, (hybrid%lexp + 1)**2), &
     820           0 :                    structconst1(atoms%nat, (2*hybrid%lexp + 1)**2))
     821           0 :          carr2 = 0; structconst1 = 0
     822           0 :          DO igpt0 = igptmin(ikpt), igptmax(ikpt)!1,hybrid%ngptm1(ikpt)
     823           0 :             igpt2 = hybrid%pgptm1(igpt0, ikpt)
     824           0 :             ix = hybrid%nbasp + igpt2
     825           0 :             igptp2 = hybrid%pgptm(igpt2, ikpt)
     826           0 :             iqnrm2 = pqnrm(igpt2, ikpt)
     827           0 :             ic2 = 0
     828           0 :             carr2 = 0
     829           0 :             DO itype2 = 1, atoms%ntype
     830           0 :                DO ineq2 = 1, atoms%neq(itype2)
     831           0 :                   ic2 = ic2 + 1
     832           0 :                   cexp = CONJG(carr2b(ic2, igpt2))
     833           0 :                   lm2 = 0
     834           0 :                   DO ic1 = 1, atoms%nat
     835           0 :                      structconst1(ic1, :) = structconst(:, ic1, ic2, ikpt)
     836             :                   END DO
     837           0 :                   DO l2 = 0, hybrid%lexp
     838           0 :                      idum = 1
     839           0 :                      DO m2 = -l2, l2
     840           0 :                         lm2 = lm2 + 1
     841           0 :                         cdum = idum*sphbesmoment(l2, itype2, iqnrm2)*cexp*carr2a(lm2, igpt2)
     842           0 :                         IF (cdum /= 0) THEN
     843             :                            lm1 = 0
     844           0 :                            DO l1 = 0, hybrid%lexp
     845           0 :                               l = l1 + l2
     846           0 :                               M = -l1 - m2 !first loop of m1
     847           0 :                               lm = l**2 + l + M
     848           0 :                               DO m1 = -l1, l1
     849           0 :                                  lm1 = lm1 + 1
     850           0 :                                  lm = lm + 1
     851           0 :                                  cdum1 = cdum*gmat(lm1, lm2)
     852           0 :                                  DO ic1 = 1, atoms%nat
     853           0 :                                     carr2(ic1, lm1) = carr2(ic1, lm1) + cdum1*structconst1(ic1, lm)
     854             :                                  END DO
     855             :                               END DO
     856             :                            END DO
     857             :                         END IF
     858           0 :                         idum = -idum !factor (-1)**(l+M)
     859             :                      END DO
     860             :                   END DO
     861             :                END DO
     862             :             END DO
     863             : 
     864           0 :             iy = hybrid%nbasp
     865           0 :             DO igpt1 = 1, igpt2
     866           0 :                iy = iy + 1
     867           0 :                igptp1 = hybrid%pgptm(igpt1, ikpt)
     868           0 :                iqnrm1 = pqnrm(igpt1, ikpt)
     869           0 :                csum = 0
     870           0 :                ic = 0
     871           0 :                DO itype = 1, atoms%ntype
     872           0 :                   DO ineq = 1, atoms%neq(itype)
     873           0 :                      ic = ic + 1
     874           0 :                      cexp = carr2b(ic, igpt1)
     875           0 :                      lm = 0
     876           0 :                      DO l = 0, hybrid%lexp
     877           0 :                         cdum = cexp*sphbesmoment(l, itype, iqnrm1)
     878           0 :                         DO M = -l, l
     879           0 :                            lm = lm + 1
     880           0 :                            csum = csum + cdum*carr2(ic, lm)*CONJG(carr2a(lm, igpt1)) ! for coulomb
     881             :                         END DO
     882             :                      END DO
     883             :                   END DO
     884             :                END DO
     885           0 :                idum = ix*(ix - 1)/2 + iy
     886           0 :                coulomb(idum, ikpt) = coulomb(idum, ikpt) + csum/cell%vol
     887             :             END DO
     888             :          END DO
     889           0 :          DEALLOCATE (carr2, carr2a, carr2b, structconst1)
     890           0 :          call timestop("loop over plane waves")
     891             :       END DO !ikpt
     892           0 :       call timestop("loop 4:")
     893             :       !     Add corrections from higher orders in (3b) to coulomb(:,1)
     894             :       ! (1) igpt1 > 1 , igpt2 > 1  (finite G vectors)
     895           0 :       call timestart("add corrections from higher orders")
     896           0 :       rdum = (4*pi_const)**(1.5)/cell%vol**2*gmat(1, 1)
     897           0 :       DO igpt0 = 1, hybrid%ngptm1(1)
     898           0 :          igpt2 = hybrid%pgptm1(igpt0, 1); IF (igpt2 == 1) CYCLE
     899           0 :          ix = hybrid%nbasp + igpt2
     900           0 :          iqnrm2 = pqnrm(igpt2, 1)
     901           0 :          igptp2 = hybrid%pgptm(igpt2, 1)
     902           0 :          q2 = MATMUL(hybrid%gptm(:, igptp2), cell%bmat)
     903           0 :          qnorm2 = SQRT(SUM(q2**2))
     904           0 :          iy = hybrid%nbasp + 1
     905           0 :          DO igpt1 = 2, igpt2
     906           0 :             iy = iy + 1
     907           0 :             idum = ix*(ix - 1)/2 + iy
     908           0 :             iqnrm1 = pqnrm(igpt1, 1)
     909           0 :             igptp1 = hybrid%pgptm(igpt1, 1)
     910           0 :             q1 = MATMUL(hybrid%gptm(:, igptp1), cell%bmat)
     911           0 :             qnorm1 = SQRT(SUM(q1**2))
     912           0 :             rdum1 = dot_PRODUCT(q1, q2)/(qnorm1*qnorm2)
     913           0 :             ic1 = 0
     914           0 :             DO itype1 = 1, atoms%ntype
     915           0 :                DO ineq1 = 1, atoms%neq(itype1)
     916           0 :                   ic1 = ic1 + 1
     917           0 :                   ic2 = 0
     918           0 :                   DO itype2 = 1, atoms%ntype
     919           0 :                      DO ineq2 = 1, atoms%neq(itype2)
     920           0 :                         ic2 = ic2 + 1
     921             :                         cdum = EXP(CMPLX(0.0, 1.0)*2*pi_const* &
     922             :                                    (-dot_PRODUCT(hybrid%gptm(:, igptp1), atoms%taual(:, ic1)) &
     923           0 :                                     + dot_PRODUCT(hybrid%gptm(:, igptp2), atoms%taual(:, ic2))))
     924             :                         coulomb(idum, 1) = coulomb(idum, 1) + rdum*cdum*( &
     925             :                                            -sphbesmoment(1, itype1, iqnrm1) &
     926             :                                            *sphbesmoment(1, itype2, iqnrm2)*rdum1/3 &
     927             :                                            - sphbesmoment(0, itype1, iqnrm1) &
     928             :                                            *sphbesmoment(2, itype2, iqnrm2)/6 &
     929             :                                            - sphbesmoment(2, itype1, iqnrm1) &
     930             :                                            *sphbesmoment(0, itype2, iqnrm2)/6 &
     931             :                                            + sphbesmoment(0, itype1, iqnrm1) &
     932             :                                            *sphbesmoment(1, itype2, iqnrm2)/qnorm2/2 &
     933             :                                            + sphbesmoment(1, itype1, iqnrm1) &
     934           0 :                                            *sphbesmoment(0, itype2, iqnrm2)/qnorm1/2)
     935             :                      END DO
     936             :                   END DO
     937             :                END DO
     938             :             END DO
     939             :          END DO
     940             :       END DO
     941             :       ! (2) igpt1 = 1 , igpt2 > 1  (first G vector vanishes, second finite)
     942           0 :       iy = hybrid%nbasp + 1
     943           0 :       DO igpt0 = 1, hybrid%ngptm1(1)
     944           0 :          igpt2 = hybrid%pgptm1(igpt0, 1); IF (igpt2 == 1) CYCLE
     945           0 :          ix = hybrid%nbasp + igpt2
     946           0 :          iqnrm2 = pqnrm(igpt2, 1)
     947           0 :          igptp2 = hybrid%pgptm(igpt2, 1)
     948           0 :          qnorm2 = qnrm(iqnrm2)
     949           0 :          idum = ix*(ix - 1)/2 + iy
     950           0 :          DO itype1 = 1, atoms%ntype
     951           0 :             DO ineq1 = 1, atoms%neq(itype1)
     952             :                ic2 = 0
     953           0 :                DO itype2 = 1, atoms%ntype
     954           0 :                   DO ineq2 = 1, atoms%neq(itype2)
     955           0 :                      ic2 = ic2 + 1
     956           0 :                      cdum = EXP(CMPLX(0.0, 1.0)*2*pi_const*dot_PRODUCT(hybrid%gptm(:, igptp2), atoms%taual(:, ic2)))
     957             :                      coulomb(idum, 1) = coulomb(idum, 1) &
     958             :                                         + rdum*cdum*atoms%rmt(itype1)**3*( &
     959             :                                         +sphbesmoment(0, itype2, iqnrm2)/30*atoms%rmt(itype1)**2 &
     960             :                                         - sphbesmoment(2, itype2, iqnrm2)/18 &
     961           0 :                                         + sphbesmoment(1, itype2, iqnrm2)/6/qnorm2)
     962             :                   END DO
     963             :                END DO
     964             :             END DO
     965             :          END DO
     966             :       END DO
     967             :       ! (2) igpt1 = 1 , igpt2 = 1  (vanishing G vectors)
     968           0 :       iy = hybrid%nbasp + 1
     969           0 :       ix = hybrid%nbasp + 1
     970           0 :       idum = ix*(ix - 1)/2 + iy
     971           0 :       DO itype1 = 1, atoms%ntype
     972           0 :          DO ineq1 = 1, atoms%neq(itype1)
     973           0 :             DO itype2 = 1, atoms%ntype
     974           0 :                DO ineq2 = 1, atoms%neq(itype2)
     975             :                   coulomb(idum, 1) = coulomb(idum, 1) &
     976             :                                      + rdum*atoms%rmt(itype1)**3*atoms%rmt(itype2)**3* &
     977           0 :                                      (atoms%rmt(itype1)**2 + atoms%rmt(itype2)**2)/90
     978             :                END DO
     979             :             END DO
     980             :          END DO
     981             :       END DO
     982           0 :       call timestop("add corrections from higher orders")
     983             : 
     984             :       !     (3c) r,r' in same MT
     985             : 
     986             :       ! Calculate sphbesintegral
     987           0 :       call timestart("sphbesintegral")
     988             :       ALLOCATE (sphbes0(-1:hybrid%lexp + 2, atoms%ntype, nqnrm),&
     989           0 :            &           carr2((hybrid%lexp + 1)**2, hybrid%maxgptm))
     990           0 :       sphbes0 = 0; carr2 = 0
     991           0 :       DO iqnrm = 1, nqnrm
     992           0 :          DO itype = 1, atoms%ntype
     993           0 :             rdum = qnrm(iqnrm)*atoms%rmt(itype)
     994           0 :             call sphbes(hybrid%lexp + 2, rdum, sphbes0(0, itype, iqnrm))
     995           0 :             IF (rdum /= 0) sphbes0(-1, itype, iqnrm) = COS(rdum)/rdum
     996             :          END DO
     997             :       END DO
     998           0 :       call timestop("sphbesintegral")
     999             : 
    1000           0 :       l_warn = (mpi%irank == 0)
    1001           0 :       call timestart("loop 2")
    1002           0 :       DO ikpt = ikptmin, ikptmax!1,nkpt
    1003             : 
    1004           0 :          DO igpt = 1, hybrid%ngptm(ikpt)
    1005           0 :             igptp = hybrid%pgptm(igpt, ikpt)
    1006           0 :             q = MATMUL(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp), cell%bmat)
    1007           0 :             call timestart("harmonics")
    1008           0 :             call ylm4(hybrid%lexp, q, carr2(:, igpt))
    1009           0 :             call timestop("harmonics")
    1010             :          END DO
    1011             : 
    1012           0 :          DO igpt0 = igptmin(ikpt), igptmax(ikpt)!1,hybrid%ngptm1(ikpt)
    1013           0 :             igpt2 = hybrid%pgptm1(igpt0, ikpt)
    1014           0 :             ix = hybrid%nbasp + igpt2
    1015           0 :             igptp2 = hybrid%pgptm(igpt2, ikpt)
    1016           0 :             iqnrm2 = pqnrm(igpt2, ikpt)
    1017           0 :             q2 = MATMUL(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp2), cell%bmat)
    1018           0 :             y2 = CONJG(carr2(:, igpt2))
    1019           0 :             iy = hybrid%nbasp
    1020           0 :             DO igpt1 = 1, igpt2
    1021           0 :                iy = iy + 1
    1022           0 :                igptp1 = hybrid%pgptm(igpt1, ikpt)
    1023           0 :                iqnrm1 = pqnrm(igpt1, ikpt)
    1024           0 :                q1 = MATMUL(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp1), cell%bmat)
    1025           0 :                y1 = carr2(:, igpt1)
    1026           0 :                cexp1 = 0
    1027           0 :                ic = 0
    1028           0 :                DO itype = 1, atoms%ntype
    1029           0 :                   DO ineq = 1, atoms%neq(itype)
    1030           0 :                      ic = ic + 1
    1031             :                      cexp1(itype) = cexp1(itype) + &
    1032             :                                     EXP(CMPLX(0.0, 1.0)*2*pi_const*dot_PRODUCT( &
    1033           0 :                                         (hybrid%gptm(:, igptp2) - hybrid%gptm(:, igptp1)), atoms%taual(:, ic)))
    1034             :                   ENDDO
    1035             :                ENDDO
    1036           0 :                lm = 0
    1037           0 :                cdum = 0
    1038           0 :                DO l = 0, hybrid%lexp
    1039           0 :                   cdum1 = 0
    1040           0 :                   DO itype = 1, atoms%ntype
    1041             :                      cdum1 = cdum1 + cexp1(itype)*sphbessel_integral( &
    1042             :                              atoms, itype, qnrm, nqnrm, &
    1043             :                              iqnrm1, iqnrm2, l, hybrid, &
    1044             :                              sphbes0, l_warn, l_warned) &
    1045           0 :                              /(2*l + 1)
    1046           0 :                      l_warn = l_warn .AND. .NOT. l_warned ! only warn once
    1047             :                   END DO
    1048           0 :                   DO M = -l, l
    1049           0 :                      lm = lm + 1
    1050           0 :                      cdum = cdum + cdum1*y1(lm)*y2(lm)
    1051             :                   ENDDO
    1052             :                ENDDO
    1053           0 :                idum = ix*(ix - 1)/2 + iy
    1054           0 :                coulomb(idum, ikpt) = coulomb(idum, ikpt) + (4*pi_const)**3*cdum/cell%vol
    1055             :             END DO
    1056             :          END DO
    1057             : 
    1058             :       END DO
    1059           0 :       call timestop("loop 2")
    1060           0 :       DEALLOCATE (carr2)
    1061             : 
    1062           0 :       IF (mpi%irank == 0) THEN
    1063           0 :          WRITE (6, '(2X,A)', advance='no') 'done'
    1064           0 :          CALL cpu_TIME(time2)
    1065           0 :          WRITE (6, '(2X,A,F8.2,A)') '( Timing:', time2 - time1, ' )'
    1066             :       END IF
    1067             : 
    1068             :       !
    1069             :       !     Symmetry-equivalent G vectors
    1070             :       !
    1071             : #ifndef CPP_NOCOULSYM
    1072             : 
    1073           0 :       IF (mpi%irank == 0) WRITE (6, '(A)', advance='no') 'Symm.-equiv. matrix elements...'
    1074           0 :       CALL cpu_TIME(time1)
    1075             :       ! All elements are needed so send all data to all processes treating the
    1076             :       ! respective k-points
    1077             : 
    1078           0 :       ALLOCATE (carr2(hybrid%maxbasm1, 2), iarr(hybrid%maxgptm))
    1079             :       ALLOCATE (nsym_gpt(hybrid%gptmd, kpts%nkpt), &
    1080           0 :                 sym_gpt(MAXVAL(nsym1), hybrid%gptmd, kpts%nkpt))
    1081           0 :       nsym_gpt = 0; sym_gpt = 0
    1082           0 :       call timestart("loop 3")
    1083           0 :       DO ikpt = ikptmin, ikptmax
    1084           0 :          carr2 = 0; iarr = 0
    1085           0 :          iarr(hybrid%pgptm1(:hybrid%ngptm1(ikpt), ikpt)) = 1
    1086           0 :          DO igpt0 = 1, hybrid%ngptm1(ikpt) !igptmin(ikpt),igptmax(ikpt)
    1087             :             lsym = ((igptmin(ikpt) <= igpt0) .AND. &
    1088           0 :                     (igptmax(ikpt) >= igpt0))
    1089           0 :             igpt2 = hybrid%pgptm1(igpt0, ikpt)
    1090           0 :             j = (hybrid%nbasp + igpt2 - 1)*(hybrid%nbasp + igpt2)/2
    1091           0 :             i = hybrid%nbasp + igpt2
    1092           0 :             carr2(1:i, 2) = coulomb(j + 1:j + i, ikpt)
    1093           0 :             j = j + i
    1094           0 :             DO i = hybrid%nbasp + igpt2 + 1, nbasm1(ikpt)
    1095           0 :                j = j + i - 1
    1096           0 :                IF (sym%invs) THEN
    1097           0 :                   carr2(i, 2) = coulomb(j, ikpt)
    1098             :                ELSE
    1099           0 :                   carr2(i, 2) = CONJG(coulomb(j, ikpt))
    1100             :                ENDIF
    1101             :             END DO
    1102           0 :             IF (lsym) THEN
    1103           0 :                ic = 1
    1104           0 :                sym_gpt(ic, igpt0, ikpt) = igpt2
    1105             :             END IF
    1106           0 :             DO isym1 = 2, nsym1(ikpt)
    1107           0 :                isym = sym1(isym1, ikpt)
    1108             :                CALL bramat_trafo( &
    1109             :                   carr2(:, 1), igpt1, &
    1110             :                   carr2(:, 2), igpt2, ikpt, isym, .FALSE., POINTER(ikpt, :, :, :), &
    1111             :                   sym, rrot(:, :, isym), invrrot(:, :, isym), hybrid, &
    1112             :                   kpts, hybrid%maxlcutm1, atoms, hybrid%lcutm1, &
    1113             :                   hybrid%nindxm1, hybrid%maxindxm1, dwgn(:, :, :, isym), &
    1114           0 :                   hybrid%nbasp, nbasm1)
    1115           0 :                IF (iarr(igpt1) == 0) THEN
    1116             :                   CALL bramat_trafo( &
    1117             :                      carr2(:, 1), igpt1, &
    1118             :                      carr2(:, 2), igpt2, ikpt, isym, .TRUE., POINTER(ikpt, :, :, :), &
    1119             :                      sym, rrot(:, :, isym), invrrot(:, :, isym), hybrid, &
    1120             :                      kpts, hybrid%maxlcutm1, atoms, hybrid%lcutm1, &
    1121             :                      hybrid%nindxm1, hybrid%maxindxm1, &
    1122           0 :                      dwgn(:, :, :, isym), hybrid%nbasp, nbasm1)
    1123           0 :                   l = (hybrid%nbasp + igpt1 - 1)*(hybrid%nbasp + igpt1)/2
    1124           0 :                   coulomb(l + 1:l + hybrid%nbasp + igpt1, ikpt) = carr2(:hybrid%nbasp + igpt1, 1)
    1125           0 :                   iarr(igpt1) = 1
    1126           0 :                   IF (lsym) THEN
    1127           0 :                      ic = ic + 1
    1128           0 :                      sym_gpt(ic, igpt0, ikpt) = igpt1
    1129             :                   END IF
    1130             :                END IF
    1131             :             END DO
    1132           0 :             nsym_gpt(igpt0, ikpt) = ic
    1133             :          END DO ! igpt0
    1134             :       END DO ! ikpt
    1135           0 :       call timestop("loop 3")
    1136           0 :       call timestart("gap 1:")
    1137           0 :       DEALLOCATE (carr2, iarr, hybrid%pgptm1)
    1138           0 :       IF (mpi%irank == 0) THEN
    1139           0 :          WRITE (6, '(2X,A)', advance='no') 'done'
    1140           0 :          CALL cpu_TIME(time2)
    1141           0 :          WRITE (6, '(2X,A,F8.2,A)') '( Timing:', time2 - time1, ' )'
    1142             :       END IF
    1143             : 
    1144             :       ! no symmetry used
    1145             : #else
    1146             : 
    1147             :       ALLOCATE (nsym_gpt(hybrid%gptmd, kpts%nkpt), sym_gpt(1, hybrid%gptmd, kpts%nkpt))
    1148             :       nsym_gpt = 1
    1149             :       DO ikpt = 1, kpts%nkpt
    1150             :          sym_gpt(1, :, ikpt) = (/(igpt0, igpt0=1, hybrid%gptmd)/)
    1151             :       END DO
    1152             : 
    1153             : #endif
    1154             : 
    1155           0 : 1     DEALLOCATE (qnrm, pqnrm)
    1156             : 
    1157           0 :       CALL cpu_TIME(time1)
    1158           0 :       IF (xcpot%is_name("hse") .OR. xcpot%is_name("vhse")) THEN
    1159             :          !
    1160             :          ! The HSE functional is realized subtracting erf/r from
    1161             :          ! the normal Coulomb matrix
    1162             :          !
    1163             : #ifdef CPP_NOSPMVEC
    1164             :          CALL change_coulombmatrix( &
    1165             :             atoms, kpts, kpts, kpts%nkpt, &
    1166             :             cell, cell, hybrid%lcutm1, hybrid%maxlcutm1, &
    1167             :             hybrid%nindxm1, hybrid%maxindxm1, hybrid, &
    1168             :             hybrid%basm1, hybrid%maxbasm1, nbasm1, sym, mpi, &
    1169             :             coulomb)
    1170             : #endif
    1171             :       ELSE
    1172           0 :          IF (ikptmin == 1) CALL subtract_sphaverage(sym, cell, atoms, kpts, hybrid, nbasm1, gridf, coulomb)
    1173             :       END IF
    1174             : 
    1175             :       ! transform Coulomb matrix to the biorthogonal set
    1176           0 :       IF (mpi%irank == 0) WRITE (6, '(A)', advance='no') 'Transform to biorthogonal set...'
    1177           0 :       CALL cpu_TIME(time1)
    1178           0 :       call timestop("gap 1:")
    1179           0 :       call timestart("calc eigenvalues olap_pw")
    1180           0 :       DO ikpt = ikptmin, ikptmax
    1181             : 
    1182             :          !calculate IR overlap-matrix
    1183           0 :          CALL olapm%alloc(sym%invs, hybrid%ngptm(ikpt), hybrid%ngptm(ikpt), 0.0)
    1184             : 
    1185           0 :          CALL olap_pw(olapm, hybrid%gptm(:, hybrid%pgptm(:hybrid%ngptm(ikpt), ikpt)), hybrid%ngptm(ikpt), atoms, cell)
    1186             : 
    1187             :          !         !calculate eigenvalues of olapm
    1188             :          !         ALLOCATE( eval(ngptm(ikpt)),evec(ngptm(ikpt),ngptm(ikpt)) )
    1189             :          !         CALL diagonalize(evec,eval,olapm)
    1190             :          !
    1191             :          !         !
    1192             :          !         ! small eigenvalues lead to inaccuries in the inversion
    1193             :          !         ! however it seems that these do not play an important role
    1194             :          !         ! for the HF exchange
    1195             :          !         ! thus we do not do a SingularValueDecomposition
    1196             :          !         !
    1197             :          !
    1198             :          !         IF( any(eval .le. 1E-06 ) ) THEN
    1199             :          ! !           WRITE(*,*) count( eval .le. 1E-06 )
    1200             :          ! !           WRITE(*,*) 'eval .le. 1E-06'
    1201             :          ! !           ALLOCATE( involapm(ngptm(ikpt),ngptm(ikpt)) )
    1202             :          !           olapm = 0
    1203             :          !           DO i = 1,ngptm(ikpt)
    1204             :          !             IF( eval(i) .le. 1E-06) CYCLE
    1205             :          !             olapm(i,i) = 1/eval(i)
    1206             :          !           END DO
    1207             :          !
    1208             :          !           ALLOCATE( invevec(ngptm(ikpt),ngptm(ikpt)) )
    1209             :          ! if (sym%invs) then
    1210             :          !           invevec =         transpose(evec)
    1211             :          ! else
    1212             :          !           invevec = conjg( transpose(evec) )
    1213             :          ! endif
    1214             :          !
    1215             :          !           olapm   = matmul(evec,matmul(olapm,invevec) )
    1216             :          !
    1217             :          !           DEALLOCATE(invevec)!,involapm)
    1218             :          !         ELSE
    1219             :          !calculate inverse overlap-matrix
    1220           0 :          CALL olapm%inverse()
    1221             :          !         END IF
    1222             : 
    1223             :          !unpack matrix coulomb
    1224           0 :          CALL coulhlp%from_packed(sym%invs, nbasm1(ikpt), REAL(coulomb(:, ikpt)), coulomb(:, ikpt))
    1225             : 
    1226           0 :          if (olapm%l_real) THEN
    1227             :             !multiply with inverse olap from right hand side
    1228           0 :             coulhlp%data_r(:, hybrid%nbasp + 1:) = MATMUL(coulhlp%data_r(:, hybrid%nbasp + 1:), olapm%data_r)
    1229             :             !multiply with inverse olap from left side
    1230           0 :             coulhlp%data_r(hybrid%nbasp + 1:, :) = MATMUL(olapm%data_r, coulhlp%data_r(hybrid%nbasp + 1:, :))
    1231             :          else
    1232             :             !multiply with inverse olap from right hand side
    1233           0 :             coulhlp%data_c(:, hybrid%nbasp + 1:) = MATMUL(coulhlp%data_c(:, hybrid%nbasp + 1:), olapm%data_c)
    1234             :             !multiply with inverse olap from left side
    1235           0 :             coulhlp%data_c(hybrid%nbasp + 1:, :) = MATMUL(olapm%data_c, coulhlp%data_c(hybrid%nbasp + 1:, :))
    1236             :          end if
    1237           0 :          coulomb(:(nbasm1(ikpt)*(nbasm1(ikpt) + 1))/2, ikpt) = coulhlp%to_packed()
    1238             : 
    1239             :       END DO
    1240           0 :       call timestop("calc eigenvalues olap_pw")
    1241             : 
    1242           0 :       IF (mpi%irank == 0) THEN
    1243           0 :          WRITE (6, '(1X,A)', advance='no') 'done'
    1244           0 :          CALL cpu_TIME(time2)
    1245           0 :          WRITE (6, '(2X,A,F8.2,A)') '( Timing:', time2 - time1, ' )'
    1246             :       END IF
    1247             : 
    1248             :       !call plot_coulombmatrix() -> code was shifted to plot_coulombmatrix.F90
    1249             : 
    1250           0 :       IF (mpi%irank == 0) WRITE (6, '(A)', advance='no') 'Writing of data to file...'
    1251           0 :       CALL cpu_TIME(time1)
    1252             : #if( !defined CPP_NOSPMVEC && !defined CPP_IRAPPROX )
    1253             :       !
    1254             :       ! rearrange coulomb matrix
    1255             :       !
    1256             : 
    1257           0 :       ALLOCATE (coulomb_mt1(hybrid%maxindxm1 - 1, hybrid%maxindxm1 - 1, 0:hybrid%maxlcutm1, atoms%ntype, 1))
    1258           0 :       ic = (hybrid%maxlcutm1 + 1)**2*atoms%nat
    1259           0 :       idum = ic + hybrid%maxgptm
    1260           0 :       idum = (idum*(idum + 1))/2
    1261           0 :       if (sym%invs) THEN
    1262           0 :          ALLOCATE (coulomb_mt2_r(hybrid%maxindxm1 - 1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1 + 1, atoms%nat, 1))
    1263           0 :          ALLOCATE (coulomb_mt3_r(hybrid%maxindxm1 - 1, atoms%nat, atoms%nat, 1))
    1264             : #ifdef CPP_IRCOULOMBAPPROX
    1265             :          ALLOCATE (coulomb_mtir_r(ic, ic + hybrid%maxgptm, 1))
    1266             :          coulomb_mtir_r = 0
    1267             :          ALLOCATE (coulombp_mtir_r(0, 0))
    1268             : #else
    1269           0 :          ALLOCATE (coulomb_mtir_r(ic + hybrid%maxgptm, ic + hybrid%maxgptm, 1))
    1270           0 :          ALLOCATE (coulombp_mtir_r(idum, 1))
    1271             : #endif
    1272             :       else
    1273           0 :          ALLOCATE (coulomb_mt2_c(hybrid%maxindxm1 - 1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1 + 1, atoms%nat, 1))
    1274           0 :          ALLOCATE (coulomb_mt3_c(hybrid%maxindxm1 - 1, atoms%nat, atoms%nat, 1))
    1275             : #ifdef CPP_IRCOULOMBAPPROX
    1276             :          ALLOCATE (coulomb_mtir_c(ic, ic + hybrid%maxgptm, 1))
    1277             :          ALLOCATE (coulombp_mtir_c(0, 0))
    1278             : #else
    1279           0 :          ALLOCATE (coulomb_mtir_c(ic + hybrid%maxgptm, ic + hybrid%maxgptm, 1))
    1280           0 :          ALLOCATE (coulombp_mtir_c(idum, 1))
    1281             : #endif
    1282             :       endif
    1283           0 :       call timestart("loop bla")
    1284           0 :       DO ikpt = ikptmin, ikptmax
    1285           0 :          ikpt0 = 1
    1286           0 :          ikpt1 = 1
    1287             :          ! initialize arrays
    1288           0 :          if (sym%invs) THEN
    1289           0 :             coulomb_mt1 = 0; coulomb_mt2_r = 0
    1290           0 :             coulomb_mt3_r = 0; coulombp_mtir_r = 0
    1291             :          else
    1292           0 :             coulomb_mt1 = 0; coulomb_mt2_c = 0
    1293           0 :             coulomb_mt3_c = 0; coulombp_mtir_c = 0
    1294             :          endif
    1295             :          ! unpack coulomb into coulhlp
    1296             : 
    1297           0 :          call coulhlp%from_packed(sym%invs, nbasm1(ikpt), real(coulomb(:, ikpt)), coulomb(:, ikpt))
    1298             : 
    1299             :          ! only one processor per k-point calculates MT convolution
    1300           0 :          IF (calc_mt(ikpt)) THEN
    1301             : 
    1302             :             !
    1303             :             ! store m-independent part of Coulomb matrix in MT spheres
    1304             :             ! in coulomb_mt1(:hybrid%nindxm1(l,itype)-1,:hybrid%nindxm1(l,itype)-1,l,itype)
    1305             :             !
    1306           0 :             call timestart("m-indep. part of coulomb mtx")
    1307           0 :             indx1 = 0
    1308           0 :             DO itype = 1, atoms%ntype
    1309           0 :                DO ineq = 1, atoms%neq(itype)
    1310           0 :                   DO l = 0, hybrid%lcutm1(itype)
    1311             : 
    1312           0 :                      IF (ineq == 1) THEN
    1313           0 :                         DO n = 1, hybrid%nindxm1(l, itype) - 1
    1314           0 :                            if (coulhlp%l_real) THEN
    1315             :                               coulomb_mt1(n, 1:hybrid%nindxm1(l, itype) - 1, l, itype, ikpt0) &
    1316           0 :                                  = coulhlp%data_r(indx1 + n, indx1 + 1:indx1 + hybrid%nindxm1(l, itype) - 1)
    1317             :                            else
    1318             :                               coulomb_mt1(n, 1:hybrid%nindxm1(l, itype) - 1, l, itype, ikpt0) &
    1319           0 :                                  = coulhlp%data_c(indx1 + n, indx1 + 1:indx1 + hybrid%nindxm1(l, itype) - 1)
    1320             :                            end if
    1321             :                         END DO
    1322             :                      END IF
    1323             : 
    1324           0 :                      indx1 = indx1 + (2*l + 1)*hybrid%nindxm1(l, itype)
    1325             :                   END DO
    1326             :                END DO
    1327             :             END DO
    1328           0 :             call timestop("m-indep. part of coulomb mtx")
    1329             : 
    1330             :             !
    1331             :             ! store m-dependent and atom-dependent part of Coulomb matrix in MT spheres
    1332             :             ! in coulomb_mt2(:hybrid%nindxm1(l,itype)-1,-l:l,l,iatom)
    1333             :             !
    1334           0 :             call timestart("m-dep. part of coulomb mtx")
    1335           0 :             indx1 = 0
    1336           0 :             iatom = 0
    1337           0 :             DO itype = 1, atoms%ntype
    1338           0 :                DO ineq = 1, atoms%neq(itype)
    1339           0 :                   iatom = iatom + 1
    1340           0 :                   DO l = 0, hybrid%lcutm1(itype)
    1341           0 :                      DO M = -l, l
    1342           0 :                         if (coulhlp%l_real) THEN
    1343             :                            coulomb_mt2_r(:hybrid%nindxm1(l, itype) - 1, M, l, iatom, ikpt0) &
    1344           0 :                               = coulhlp%data_r(indx1 + 1:indx1 + hybrid%nindxm1(l, itype) - 1, indx1 + hybrid%nindxm1(l, itype))
    1345             :                         else
    1346             :                            coulomb_mt2_c(:hybrid%nindxm1(l, itype) - 1, M, l, iatom, ikpt0) &
    1347           0 :                               = coulhlp%data_c(indx1 + 1:indx1 + hybrid%nindxm1(l, itype) - 1, indx1 + hybrid%nindxm1(l, itype))
    1348             :                         endif
    1349             : 
    1350           0 :                         indx1 = indx1 + hybrid%nindxm1(l, itype)
    1351             : 
    1352             :                      END DO
    1353             :                   END DO
    1354             :                END DO
    1355             :             END DO
    1356           0 :             call timestop("m-dep. part of coulomb mtx")
    1357             : 
    1358             :             !
    1359             :             ! due to the subtraction of the divergent part at the Gamma point
    1360             :             ! additional contributions occur
    1361             :             !
    1362           0 :             call timestart("gamma point treatment")
    1363           0 :             IF (ikpt == 1) THEN
    1364             :                !
    1365             :                ! store the contribution of the G=0 plane wave with the MT l=0 functions in
    1366             :                ! coulomb_mt2(:hybrid%nindxm1(l=0,itype),0,hybrid%maxlcutm1+1,iatom)
    1367             :                !
    1368           0 :                ic = 0
    1369           0 :                iatom = 0
    1370           0 :                DO itype = 1, atoms%ntype
    1371           0 :                   DO ineq = 1, atoms%neq(itype)
    1372           0 :                      iatom = iatom + 1
    1373           0 :                      DO n = 1, hybrid%nindxm1(0, itype) - 1
    1374           0 :                         if (coulhlp%l_real) THEN
    1375           0 :                            coulomb_mt2_r(n, 0, hybrid%maxlcutm1 + 1, iatom, ikpt0) = coulhlp%data_r(ic + n, hybrid%nbasp + 1)
    1376             :                         else
    1377           0 :                            coulomb_mt2_c(n, 0, hybrid%maxlcutm1 + 1, iatom, ikpt0) = coulhlp%data_c(ic + n, hybrid%nbasp + 1)
    1378             :                         endif
    1379             :                      END DO
    1380           0 :                      ic = ic + SUM((/((2*l + 1)*hybrid%nindxm1(l, itype), l=0, hybrid%lcutm1(itype))/))
    1381             :                   END DO
    1382             :                END DO
    1383             : 
    1384             :                !
    1385             :                ! store the contributions between the MT s-like functions at atom1 and
    1386             :                ! and the constant function at a different atom2
    1387             :                !
    1388           0 :                iatom = 0
    1389           0 :                ic = 0
    1390           0 :                DO itype = 1, atoms%ntype
    1391           0 :                   ishift = SUM((/((2*l + 1)*hybrid%nindxm1(l, itype), l=0, hybrid%lcutm1(itype))/))
    1392           0 :                   DO ineq = 1, atoms%neq(itype)
    1393           0 :                      iatom = iatom + 1
    1394           0 :                      ic1 = ic + hybrid%nindxm1(0, itype)
    1395             : 
    1396           0 :                      iatom1 = 0
    1397           0 :                      ic2 = 0
    1398           0 :                      DO itype1 = 1, atoms%ntype
    1399           0 :                         ishift1 = SUM((/((2*l1 + 1)*hybrid%nindxm1(l1, itype1), l1=0, hybrid%lcutm1(itype1))/))
    1400           0 :                         DO ineq1 = 1, atoms%neq(itype1)
    1401           0 :                            iatom1 = iatom1 + 1
    1402           0 :                            ic3 = ic2 + 1
    1403           0 :                            ic4 = ic3 + hybrid%nindxm1(0, itype1) - 2
    1404             : 
    1405           0 :                            IF (sym%invs) THEN
    1406           0 :                               coulomb_mt3_r(:hybrid%nindxm1(0, itype1) - 1, iatom, iatom1, ikpt0) = coulhlp%data_r(ic1, ic3:ic4)
    1407             :                            ELSE
    1408             :                               coulomb_mt3_c(:hybrid%nindxm1(0, itype1) - 1, iatom, iatom1, ikpt0) &
    1409           0 :                                  = CONJG(coulhlp%data_c(ic1, ic3:ic4))
    1410             :                            ENDIF
    1411           0 :                            ic2 = ic2 + ishift1
    1412             :                         END DO
    1413             :                      END DO
    1414             : 
    1415           0 :                      ic = ic + ishift
    1416             :                   END DO
    1417             :                END DO
    1418             : 
    1419             :                !test
    1420           0 :                iatom = 0
    1421           0 :                DO itype = 1, atoms%ntype
    1422           0 :                   DO ineq = 1, atoms%neq(itype)
    1423           0 :                      iatom = iatom + 1
    1424           0 :                      if (sym%invs) THEN
    1425           0 :                         IF (MAXVAL(ABS(coulomb_mt2_r(:hybrid%nindxm1(0, itype) - 1, 0, 0, &
    1426             :                                                      iatom, ikpt0) &
    1427             :                                        - coulomb_mt3_r(:hybrid%nindxm1(0, itype) - 1, iatom, &
    1428             :                                                        iatom, ikpt0))) &
    1429             :                             > 1E-08) &
    1430           0 :                            call judft_error('coulombmatrix: coulomb_mt2 and coulomb_mt3 are inconsistent')
    1431             : 
    1432             :                      else
    1433           0 :                         IF (MAXVAL(ABS(coulomb_mt2_c(:hybrid%nindxm1(0, itype) - 1, 0, 0, &
    1434             :                                                      iatom, ikpt0) &
    1435             :                                        - coulomb_mt3_c(:hybrid%nindxm1(0, itype) - 1, iatom, &
    1436             :                                                        iatom, ikpt0))) &
    1437             :                             > 1E-08) &
    1438           0 :                            call judft_error('coulombmatrix: coulomb_mt2 and coulomb_mt3 are inconsistent')
    1439             :                      endif
    1440             :                   END DO
    1441             :                END DO
    1442             :             END IF
    1443           0 :             call timestop("gamma point treatment")
    1444             : 
    1445             :          END IF ! calc_mt
    1446             : 
    1447             :          !
    1448             :          ! add the residual MT contributions, i.e. those functions with an moment,
    1449             :          ! to the matrix coulomb_mtir, which is fully occupied
    1450             :          !
    1451           0 :          call timestart("residual MT contributions")
    1452           0 :          ic = 0
    1453           0 :          DO itype = 1, atoms%ntype
    1454           0 :             DO ineq = 1, atoms%neq(itype)
    1455           0 :                DO l = 0, hybrid%lcutm1(itype)
    1456           0 :                   DO M = -l, l
    1457           0 :                      ic = ic + 1
    1458             :                   END DO
    1459             :                END DO
    1460             :             END DO
    1461             :          END DO
    1462             : 
    1463           0 :          indx1 = 0; indx2 = 0; indx3 = 0; indx4 = 0
    1464           0 :          DO itype = 1, atoms%ntype
    1465           0 :             DO ineq = 1, atoms%neq(itype)
    1466           0 :                DO l = 0, hybrid%lcutm1(itype)
    1467           0 :                   DO M = -l, l
    1468           0 :                      indx1 = indx1 + 1
    1469           0 :                      indx3 = indx3 + hybrid%nindxm1(l, itype)
    1470             : 
    1471           0 :                      indx2 = 0
    1472           0 :                      indx4 = 0
    1473           0 :                      DO itype1 = 1, atoms%ntype
    1474           0 :                         DO ineq1 = 1, atoms%neq(itype1)
    1475           0 :                            DO l1 = 0, hybrid%lcutm1(itype1)
    1476           0 :                               DO m1 = -l1, l1
    1477           0 :                                  indx2 = indx2 + 1
    1478           0 :                                  indx4 = indx4 + hybrid%nindxm1(l1, itype1)
    1479           0 :                                  IF (indx4 < indx3) CYCLE
    1480           0 :                                  IF (calc_mt(ikpt)) THEN
    1481           0 :                                     IF (sym%invs) THEN
    1482           0 :                                        coulomb_mtir_r(indx1, indx2, ikpt1) = coulhlp%data_r(indx3, indx4)
    1483           0 :                                        coulomb_mtir_r(indx2, indx1, ikpt1) = coulomb_mtir_r(indx1, indx2, ikpt1)
    1484             :                                     ELSE
    1485           0 :                                        coulomb_mtir_c(indx1, indx2, ikpt1) = coulhlp%data_c(indx3, indx4)
    1486           0 :                                        coulomb_mtir_c(indx2, indx1, ikpt1) = CONJG(coulomb_mtir_c(indx1, indx2, ikpt1))
    1487             :                                     ENDIF
    1488             :                                  END IF
    1489             :                               END DO
    1490             :                            END DO
    1491             :                         END DO
    1492             :                      END DO
    1493             : 
    1494           0 :                      DO igpt = 1, hybrid%ngptm(ikpt)
    1495           0 :                         indx2 = indx2 + 1
    1496           0 :                         IF (sym%invs) THEN
    1497           0 :                            coulomb_mtir_r(indx1, indx2, ikpt1) = coulhlp%data_r(indx3, hybrid%nbasp + igpt)
    1498             : #if !defined CPP_IRCOULOMBAPPROX
    1499           0 :                            coulomb_mtir_r(indx2, indx1, ikpt1) = coulomb_mtir_r(indx1, indx2, ikpt1)
    1500             : #endif
    1501             :                         ELSE
    1502           0 :                            coulomb_mtir_c(indx1, indx2, ikpt1) = coulhlp%data_c(indx3, hybrid%nbasp + igpt)
    1503             : #if !defined CPP_IRCOULOMBAPPROX
    1504           0 :                            coulomb_mtir_c(indx2, indx1, ikpt1) = CONJG(coulomb_mtir_c(indx1, indx2, ikpt1))
    1505             : #endif
    1506             :                         ENDIF
    1507             : 
    1508             :                      END DO
    1509             : 
    1510             :                   END DO
    1511             :                END DO
    1512             :             END DO
    1513             :          END DO
    1514           0 :          call timestop("residual MT contributions")
    1515             : 
    1516           0 :          IF (indx1 /= ic) STOP 'coulombmatrix: error index counting'
    1517             : 
    1518             :          !
    1519             :          ! add ir part to the matrix coulomb_mtir
    1520             :          !
    1521             : #ifndef CPP_IRCOULOMBAPPROX
    1522           0 :          if (sym%invs) THEN
    1523             :             coulomb_mtir_r(ic + 1:ic + hybrid%ngptm(ikpt), ic + 1:ic + hybrid%ngptm(ikpt), ikpt1) &
    1524           0 :                = coulhlp%data_r(hybrid%nbasp + 1:nbasm1(ikpt), hybrid%nbasp + 1:nbasm1(ikpt))
    1525           0 :             ic2 = indx1 + hybrid%ngptm(ikpt)
    1526           0 :             coulombp_mtir_r(:ic2*(ic2 + 1)/2, ikpt0) = packmat(coulomb_mtir_r(:ic2, :ic2, ikpt1))
    1527             :          else
    1528             :             coulomb_mtir_c(ic + 1:ic + hybrid%ngptm(ikpt), ic + 1:ic + hybrid%ngptm(ikpt), ikpt1) &
    1529           0 :                = coulhlp%data_c(hybrid%nbasp + 1:nbasm1(ikpt), hybrid%nbasp + 1:nbasm1(ikpt))
    1530           0 :             ic2 = indx1 + hybrid%ngptm(ikpt)
    1531           0 :             coulombp_mtir_c(:ic2*(ic2 + 1)/2, ikpt0) = packmat(coulomb_mtir_c(:ic2, :ic2, ikpt1))
    1532             :          end if
    1533             : #endif
    1534           0 :          call timestart("write coulomb_spm")
    1535           0 :          if (sym%invs) THEN
    1536             : #ifdef CPP_IRCOULOMBAPPROX
    1537             :             call write_coulomb_spm_r(ikpt, coulomb_mt1(:, :, :, :, 1), coulomb_mt2_r(:, :, :, :, 1), &
    1538             :                                      coulomb_mt3_r(:, :, :, 1), coulomb_mtir_r(:, 1))
    1539             : #else
    1540             :             CALL write_coulomb_spm_r(ikpt, coulomb_mt1(:, :, :, :, 1), coulomb_mt2_r(:, :, :, :, 1), &
    1541           0 :                                      coulomb_mt3_r(:, :, :, 1), coulombp_mtir_r(:, 1))
    1542             : !!$       print *,"DEBUG"
    1543             : !!$       DO n1=1,SIZE(coulomb_mt1,1)
    1544             : !!$          DO n2=1,SIZE(coulomb_mt1,2)
    1545             : !!$             DO i=1,SIZE(coulomb_mt1,3)
    1546             : !!$                DO j=1,SIZE(coulomb_mt1,4)
    1547             : !!$                   WRITE(732,*) n1,n2,i-1,j,coulomb_mt2_r(n1,n2,i-1,j,1)
    1548             : !!$                ENDDO
    1549             : !!$             ENDDO
    1550             : !!$          ENDDO
    1551             : !!$       ENDDO
    1552             : #endif
    1553             :          else
    1554             : #ifdef CPP_IRCOULOMBAPPROX
    1555             :             call write_coulomb_spm_c(ikpt, coulomb_mt1(:, :, :, :, 1), coulomb_mt2_c(:, :, :, :, 1), &
    1556             :                                      coulomb_mt3_c(:, :, :, 1), coulomb_mtir_c(:, 1))
    1557             : #else
    1558             :             call write_coulomb_spm_c(ikpt, coulomb_mt1(:, :, :, :, 1), coulomb_mt2_c(:, :, :, :, 1), &
    1559           0 :                                      coulomb_mt3_c(:, :, :, 1), coulombp_mtir_c(:, 1))
    1560             : #endif
    1561             :          endif
    1562           0 :          call timestop("write coulomb_spm")
    1563             : 
    1564             :       END DO ! ikpt
    1565           0 :       call timestop("loop bla")
    1566             : 
    1567           0 :       if (sym%invs) THEN
    1568           0 :          DEALLOCATE (coulomb_mt1, coulomb_mt2_r, coulomb_mt3_r, coulomb_mtir_r, coulombp_mtir_r)
    1569             :       else
    1570           0 :          DEALLOCATE (coulomb_mt1, coulomb_mt2_c, coulomb_mt3_c, coulomb_mtir_c, coulombp_mtir_c)
    1571             :       end if
    1572             : #else
    1573             :       !write coulomb matrix to direct access file coulomb
    1574             :       DO i = 1, kpts%nkpt
    1575             :          call write_coulomb(i, sym%invs, coulomb(:, i))
    1576             :       END DO
    1577             : #endif
    1578             : 
    1579           0 :       IF (mpi%irank == 0) THEN
    1580           0 :          WRITE (6, '(7X,A)', advance='no') 'done'
    1581           0 :          CALL cpu_TIME(time2)
    1582           0 :          WRITE (6, '(2X,A,F8.2,A)') '( Timing:', time2 - time1, ' )'
    1583             :       END IF
    1584             : 
    1585           0 :       CALL timestop("Coulomb matrix setup")
    1586             : 
    1587           0 :    END SUBROUTINE coulombmatrix
    1588             : 
    1589             :    !     Calculate body of Coulomb matrix at Gamma point: v_IJ = SUM(G) c^*_IG c_JG 4*pi/G**2 .
    1590             :    !     For this we must subtract from coulomb(:,1) the spherical average of a term that comes
    1591             :    !     from the fact that MT functions have k-dependent Fourier coefficients (see script).
    1592           0 :    SUBROUTINE subtract_sphaverage(sym, cell, atoms, kpts, hybrid, nbasm1, gridf, coulomb)
    1593             : 
    1594             :       USE m_types
    1595             :       USE m_constants
    1596             :       USE m_wrapper
    1597             :       USE m_trafo
    1598             :       USE m_util
    1599             :       USE m_olap
    1600             :       IMPLICIT NONE
    1601             : 
    1602             :       TYPE(t_sym), INTENT(IN)    :: sym
    1603             :       TYPE(t_cell), INTENT(IN)    :: cell
    1604             :       TYPE(t_atoms), INTENT(IN)    :: atoms
    1605             :       TYPE(t_kpts), INTENT(IN)    :: kpts
    1606             :       TYPE(t_hybrid), INTENT(IN)    :: hybrid
    1607             : 
    1608             :       INTEGER, INTENT(IN)    :: nbasm1(kpts%nkptf)
    1609             :       REAL, INTENT(IN)    :: gridf(:, :)
    1610             :       COMPLEX, INTENT(INOUT) :: coulomb(hybrid%maxbasm1*(hybrid%maxbasm1 + 1)/2, kpts%nkpt)
    1611             : 
    1612             :       ! - local scalars -
    1613             :       INTEGER               :: l, i, j, n, nn, itype, ieq, M
    1614             : 
    1615             :       ! - local arrays -
    1616           0 :       TYPE(t_mat) :: olap
    1617             :       !COMPLEX , ALLOCATABLE :: constfunc(:)  !can also be real in inversion case
    1618           0 :       COMPLEX      :: coeff(nbasm1(1)), cderiv(nbasm1(1), -1:1), claplace(nbasm1(1))
    1619             : 
    1620           0 :       CALL olap%alloc(sym%invs, hybrid%ngptm(1), hybrid%ngptm(1), 0.)
    1621             : 
    1622           0 :       n = nbasm1(1)
    1623           0 :       nn = n*(n + 1)/2
    1624           0 :       CALL olap_pw(olap, hybrid%gptm(:, hybrid%pgptm(:hybrid%ngptm(1), 1)), hybrid%ngptm(1), atoms, cell)
    1625             : 
    1626             :       ! Define coefficients (coeff) and their derivatives (cderiv,claplace)
    1627           0 :       coeff = 0
    1628           0 :       cderiv = 0
    1629           0 :       claplace = 0
    1630           0 :       j = 0
    1631           0 :       DO itype = 1, atoms%ntype
    1632           0 :          DO ieq = 1, atoms%neq(itype)
    1633           0 :             DO l = 0, hybrid%lcutm1(itype)
    1634           0 :                DO M = -l, l
    1635           0 :                   DO i = 1, hybrid%nindxm1(l, itype)
    1636           0 :                      j = j + 1
    1637           0 :                      IF (l == 0) THEN
    1638             :                         coeff(j) = SQRT(4*pi_const) &
    1639             :                                    *intgrf(atoms%rmsh(:, itype)*hybrid%basm1(:, i, 0, itype), &
    1640             :                                            atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf) &
    1641           0 :                                    /SQRT(cell%vol)
    1642             : 
    1643             :                         claplace(j) = -SQRT(4*pi_const) &
    1644             :                                       *intgrf(atoms%rmsh(:, itype)**3*hybrid%basm1(:, i, 0, itype), &
    1645             :                                               atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf) &
    1646           0 :                                       /SQRT(cell%vol)
    1647             : 
    1648           0 :                      ELSE IF (l == 1) THEN
    1649             :                         cderiv(j, M) = -SQRT(4*pi_const/3)*CMPLX(0.0, 1.0) &
    1650             :                                        *intgrf(atoms%rmsh(:, itype)**2*hybrid%basm1(:, i, 1, itype), &
    1651             :                                                atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf) &
    1652           0 :                                        /SQRT(cell%vol)
    1653             :                      END IF
    1654             :                   END DO
    1655             :                END DO
    1656             :             END DO
    1657             :          END DO
    1658             :       END DO
    1659           0 :       IF (olap%l_real) THEN
    1660           0 :          coeff(hybrid%nbasp + 1:n) = olap%data_r(1, 1:n - hybrid%nbasp)
    1661             :       else
    1662           0 :          coeff(hybrid%nbasp + 1:n) = olap%data_c(1, 1:n - hybrid%nbasp)
    1663             :       END IF
    1664           0 :       IF (sym%invs) THEN
    1665             :          CALL symmetrize(coeff, 1, nbasm1(1), 2, .FALSE., &
    1666             :                          atoms, hybrid%lcutm1, hybrid%maxlcutm1, &
    1667           0 :                          hybrid%nindxm1, sym)
    1668             :          CALL symmetrize(claplace, 1, nbasm1(1), 2, .FALSE., &
    1669             :                          atoms, hybrid%lcutm1, hybrid%maxlcutm1, &
    1670           0 :                          hybrid%nindxm1, sym)
    1671             :          CALL symmetrize(cderiv(:, -1), 1, nbasm1(1), 2, .FALSE., &
    1672             :                          atoms, hybrid%lcutm1, hybrid%maxlcutm1, &
    1673           0 :                          hybrid%nindxm1, sym)
    1674             :          CALL symmetrize(cderiv(:, 0), 1, nbasm1(1), 2, .FALSE., &
    1675             :                          atoms, hybrid%lcutm1, hybrid%maxlcutm1, &
    1676           0 :                          hybrid%nindxm1, sym)
    1677             :          CALL symmetrize(cderiv(:, 1), 1, nbasm1(1), 2, .FALSE., &
    1678             :                          atoms, hybrid%lcutm1, hybrid%maxlcutm1, &
    1679           0 :                          hybrid%nindxm1, sym)
    1680             :       ENDIF
    1681             :       ! Subtract head contributions from coulomb(:nn,1) to obtain the body
    1682             :       l = 0
    1683           0 :       DO j = 1, n
    1684           0 :          DO i = 1, j
    1685           0 :             l = l + 1
    1686             :             coulomb(l, 1) = coulomb(l, 1) - 4*pi_const/3 &
    1687             :                             *(dot_PRODUCT(cderiv(i, :), cderiv(j, :)) &
    1688             :                               + (CONJG(coeff(i))*claplace(j) &
    1689           0 :                                  + CONJG(claplace(i))*coeff(j))/2)
    1690             :          END DO
    1691             :       END DO
    1692           0 :       coeff(hybrid%nbasp + 1) = 1.0
    1693           0 :       coeff(hybrid%nbasp + 2:) = 0.0
    1694           0 :       IF (sym%invs) THEN
    1695             : 
    1696             :          CALL desymmetrize(coeff, 1, nbasm1(1), 2, &
    1697             :                            atoms, hybrid%lcutm1, hybrid%maxlcutm1, &
    1698           0 :                            hybrid%nindxm1, sym)
    1699             :          CALL symmetrize(coeff, nbasm1(1), 1, 1, .FALSE., &
    1700             :                          atoms, hybrid%lcutm1, hybrid%maxlcutm1, &
    1701           0 :                          hybrid%nindxm1, sym)
    1702             :       ENDIF
    1703             :       ! Explicit normalization here in order to prevent failure of the diagonalization in diagonalize_coulomb
    1704             :       ! due to inaccuracies in the overlap matrix (which can make it singular).
    1705             :       !constfunc = coeff / SQRT ( ( SUM(ABS(coeff(:hybrid%nbasp))**2) + dotprod ( coeff(hybrid%nbasp+1:), MATMUL(olap,coeff(hybrid%nbasp+1:)) ) ) )
    1706             : 
    1707           0 :    END SUBROUTINE subtract_sphaverage
    1708             : 
    1709             :    !     -----------------------------------------------------------------------------------------------
    1710             : 
    1711             :    !     Calculates the structure constant
    1712             :    !                                                        1               *      ^
    1713             :    !     structconst(lm,ic1,ic2,k) = SUM exp(ikT) -----------------------  Y  ( T + R(ic) )
    1714             :    !                                  T           | T + R(ic1) - R(ic2) |   lm
    1715             :    !
    1716             :    !     with T = lattice vectors
    1717             :    !
    1718             :    !     An Ewald summation method devised by O.K. Andersen is used for l<5
    1719             :    !     (see e.g. H.L. Skriver, "The LMTO method", Springer 1984).
    1720             :    !     (The real-space function G can be calculated with gfunction.f)
    1721             :    !
    1722             : 
    1723             :    ! Convergence parameter
    1724             : #define CONVPARAM 1e-18
    1725             :    ! Do some additional shells ( real-space and Fourier-space sum )
    1726             : #define ADDSHELL1 40
    1727             : #define ADDSHELL2 0
    1728             : 
    1729           0 :    SUBROUTINE structureconstant(structconst, cell, hybrid, atoms, kpts, mpi)
    1730             : 
    1731             :       USE m_constants, ONLY: pi_const
    1732             :       USE m_util, ONLY: rorderp, rorderpf
    1733             :       USE m_types
    1734             :       USE m_juDFT
    1735             :       use m_ylm
    1736             :       IMPLICIT NONE
    1737             : 
    1738             :       TYPE(t_mpi), INTENT(IN)   :: mpi
    1739             : 
    1740             :       TYPE(t_hybrid), INTENT(INOUT)   :: hybrid
    1741             : 
    1742             :       TYPE(t_cell), INTENT(IN)   :: cell
    1743             : 
    1744             :       TYPE(t_atoms), INTENT(IN)   :: atoms
    1745             :       TYPE(t_kpts), INTENT(IN)    :: kpts
    1746             :       ! - scalars -
    1747             : 
    1748             :       ! - arrays -
    1749             :       COMPLEX, INTENT(INOUT)   ::  structconst((2*hybrid%lexp + 1)**2, atoms%nat, atoms%nat, kpts%nkpt)
    1750             : 
    1751             :       ! - local scalars -
    1752             :       INTEGER                   ::  i, ic1, ic2, lm, ikpt, l, ishell, nshell
    1753             :       INTEGER                   ::  m
    1754             :       INTEGER                   ::  nptsh, maxl
    1755             : 
    1756             :       REAL                      ::  rad, rrad, rdum
    1757             :       REAL                      ::  a, a1, aa
    1758             :       REAL                      ::  pref, rexp
    1759             :       REAL                      ::  time1, time2
    1760             :       REAL                      ::  scale
    1761             : 
    1762             :       COMPLEX                   ::  cdum, cexp
    1763             : 
    1764             :       LOGICAL, SAVE          ::  first = .TRUE.
    1765             :       ! - local arrays -
    1766           0 :       INTEGER                   ::  conv(0:2*hybrid%lexp)
    1767           0 :       INTEGER, ALLOCATABLE     ::  pnt(:), ptsh(:, :)
    1768             : 
    1769             :       REAL                      ::  rc(3), ra(3), k(3), ki(3), ka(3)
    1770           0 :       REAL                      ::  convpar(0:2*hybrid%lexp), g(0:2*hybrid%lexp)
    1771           0 :       REAL, ALLOCATABLE     ::  radsh(:)
    1772             : 
    1773           0 :       COMPLEX                   ::  y((2*hybrid%lexp + 1)**2)
    1774           0 :       COMPLEX                   ::  shlp((2*hybrid%lexp + 1)**2, kpts%nkpt)
    1775             : 
    1776           0 :       call timestart("calc struc_const.")
    1777             : 
    1778           0 :       IF (mpi%irank /= 0) first = .FALSE.
    1779             : 
    1780           0 :       rdum = cell%vol**(1.0/3) ! define "average lattice parameter"
    1781             : 
    1782             :       ! ewaldlambda = ewaldscale
    1783           0 :       scale = hybrid%ewaldlambda/rdum
    1784             : 
    1785             :       !       lambda = ewaldlambda / rdum
    1786             : 
    1787           0 :       pref = 4*pi_const/(scale**3*cell%vol)
    1788             : 
    1789           0 :       DO l = 0, 2*hybrid%lexp
    1790           0 :          convpar(l) = CONVPARAM/scale**(l + 1)
    1791             :       END DO
    1792             : 
    1793           0 :       IF (first) THEN
    1794           0 :          WRITE (6, '(//A)') '### subroutine: structureconstant ###'
    1795           0 :          WRITE (6, '(/A)') 'Real-space sum:'
    1796             :       END IF
    1797             : 
    1798             :       !
    1799             :       !     Determine cutoff radii for real-space and Fourier-space summation
    1800             :       ! (1) real space
    1801           0 :       call timestart("determine cutoff radii")
    1802           0 :       a = 1
    1803           0 : 1     rexp = EXP(-a)
    1804           0 :       g(0) = rexp/a*(1 + a*11/16*(1 + a*3/11*(1 + a/9)))
    1805           0 :       g(1) = rexp/a**2*(1 + a*(1 + a/2*(1 + a*7/24*(1 + a/7))))
    1806             :       g(2) = rexp/a**3*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a*3/16 &
    1807           0 :                                                           *(1 + a/9))))))
    1808             :       g(3) = rexp/a**4*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6 &
    1809           0 :                                                                    *(1 + a/8)))))))
    1810             :       g(4) = rexp/a**5*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6 &
    1811           0 :                                                                    *(1 + a/7*(1 + a/8*(1 + a/10)))))))))
    1812             :       g(5) = rexp/a**6*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6 &
    1813           0 :                                                                    *(1 + a/7*(1 + a/8*(1 + a/9*(1 + a/10))))))))))
    1814             :       g(6) = rexp/a**7*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6 &
    1815           0 :                                                                    *(1 + a/7*(1 + a/8*(1 + a/9*(1 + a/10*(1 + a/11*(1 + a/12))))))))))))
    1816             :       g(7) = rexp/a**8*(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6 &
    1817           0 :                                                                    *(1 + a/7*(1 + a/8*(1 + a/9*(1 + a/10*(1 + a/11*(1 + a/12*(1 + a/13)))))))))))))
    1818           0 :       DO l = 8, 2*hybrid%lexp
    1819           0 :          g(l) = a**(-l - 1)
    1820             :       END DO
    1821           0 :       IF (ANY(g > convpar/10)) THEN ! one digit more accuracy for real-space sum
    1822           0 :          a = a + 1
    1823           0 :          GOTO 1
    1824             :       END IF
    1825           0 :       rad = a/scale
    1826           0 :       call timestop("determine cutoff radii")
    1827             : 
    1828             :       ! (2) Fourier space
    1829           0 :       call timestart("fourier space")
    1830           0 :       a = 1
    1831           0 : 2     aa = (1 + a**2)**(-1)
    1832           0 :       g(0) = pref*aa**4/a**2
    1833           0 :       g(1) = pref*aa**4/a
    1834           0 :       g(2) = pref*aa**5/3
    1835           0 :       g(3) = pref*aa**5*a/15
    1836           0 :       g(4) = pref*aa**6*a**2/105
    1837           0 :       g(5) = pref*aa**6*a**3/945
    1838           0 :       g(6) = pref*aa**7*a**4/10395
    1839           0 :       g(7) = pref*aa**7*a**5/135135
    1840           0 :       IF (ANY(g > convpar)) THEN
    1841           0 :          a = a + 1
    1842           0 :          GOTO 2
    1843             :       END IF
    1844           0 :       rrad = a*scale
    1845           0 :       call timestop("fourier space")
    1846             : 
    1847           0 :       IF (first) THEN
    1848           0 :          WRITE (6, '(/A,2F10.5)') 'Cutoff radii: ', rad, rrad
    1849           0 :          WRITE (6, '(/A)') 'Real-space sum'
    1850             :       END IF
    1851             : 
    1852             :       !
    1853             :       !     Determine atomic shells
    1854           0 :       call timestart("determine atomic shell")
    1855           0 :       CALL getshells(ptsh, nptsh, radsh, nshell, rad, cell%amat, first)
    1856           0 :       call timestop("determine atomic shell")
    1857             : 
    1858           0 :       ALLOCATE (pnt(nptsh))
    1859           0 :       structconst = 0
    1860             : 
    1861             :       !
    1862             :       !     Real-space sum
    1863             :       !
    1864           0 :       CALL cpu_TIME(time1)
    1865             : 
    1866           0 :       call timestart("realspace sum")
    1867           0 :       DO ic2 = 1, atoms%nat
    1868           0 :          DO ic1 = 1, atoms%nat
    1869           0 :             IF (ic2 /= 1 .AND. ic1 == ic2) CYCLE
    1870           0 :             rc = MATMUL(cell%amat, (atoms%taual(:, ic2) - atoms%taual(:, ic1)))
    1871           0 :             DO i = 1, nptsh
    1872           0 :                ra = MATMUL(cell%amat, ptsh(:, i)) + rc
    1873           0 :                a = SQRT(SUM(ra**2))
    1874           0 :                radsh(i) = a
    1875             :             END DO
    1876           0 :             call timestart("rorderpf")
    1877           0 :             CALL rorderpf(pnt, radsh, nptsh, MAX(0, INT(LOG(nptsh*0.001)/LOG(2.0))))
    1878           0 :             call timestop("rorderpf")
    1879           0 :             ptsh = ptsh(:, pnt)
    1880           0 :             radsh = radsh(pnt)
    1881           0 :             maxl = 2*hybrid%lexp
    1882           0 :             a1 = HUGE(a1)  ! stupid initial value
    1883           0 :             ishell = 1
    1884           0 :             conv = HUGE(i)
    1885           0 :             shlp = 0
    1886           0 :             DO i = 1, nptsh
    1887           0 :                IF (ALL(conv /= HUGE(i))) EXIT
    1888           0 :                IF (i /= 1) THEN
    1889           0 :                   IF (ABS(radsh(i) - radsh(i - 1)) > 1e-10) ishell = ishell + 1
    1890             :                ENDIF
    1891           0 :                ra = MATMUL(cell%amat, ptsh(:, i)) + rc
    1892           0 :                a = scale*SQRT(SUM(ra**2))
    1893           0 :                IF (a == 0) THEN
    1894             :                   CYCLE
    1895           0 :                ELSE IF (ABS(a - a1) > 1e-10) THEN
    1896           0 :                   a1 = a
    1897           0 :                   rexp = EXP(-a)
    1898           0 :                   IF (ishell <= conv(0)) g(0) = rexp/a &
    1899           0 :                                                 *(1 + a*11/16*(1 + a*3/11*(1 + a/9)))
    1900           0 :                   IF (ishell <= conv(1)) g(1) = rexp/a**2 &
    1901           0 :                                                 *(1 + a*(1 + a/2*(1 + a*7/24*(1 + a/7))))
    1902           0 :                   IF (ishell <= conv(2)) g(2) = rexp/a**3 &
    1903           0 :                                                 *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a*3/16*(1 + a/9))))))
    1904           0 :                   IF (ishell <= conv(3)) g(3) = rexp/a**4 &
    1905           0 :                                                 *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6*(1 + a/8)))))))
    1906           0 :                   IF (ishell <= conv(4)) g(4) = rexp/a**5 &
    1907             :                                                 *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6*(1 + a/7*(1 + a/8 &
    1908           0 :                                                                                                                *(1 + a/10)))))))))
    1909           0 :                   IF (ishell <= conv(5)) g(5) = rexp/a**6 &
    1910             :                                                 *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6*(1 + a/7*(1 + a/8*(1 + a/9 &
    1911           0 :                                                                                                                         *(1 + a/10))))))))))
    1912           0 :                   IF (ishell <= conv(6)) g(6) = rexp/a**7 &
    1913             :                                                 *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6*(1 + a/7*(1 + a/8*(1 + a/9 &
    1914           0 :                                                                                                                         *(1 + a/10*(1 + a/11*(1 + a/12))))))))))))
    1915           0 :                   IF (ishell <= conv(7)) g(7) = rexp/a**8 &
    1916             :                                                 *(1 + a*(1 + a/2*(1 + a/3*(1 + a/4*(1 + a/5*(1 + a/6*(1 + a/7*(1 + a/8*(1 + a/9 &
    1917           0 :                                                                                                                         *(1 + a/10*(1 + a/11*(1 + a/12*(1 + a/13)))))))))))))
    1918           0 :                   DO l = 8, maxl
    1919           0 :                      IF (ishell <= conv(l)) g(l) = a**(-l - 1)
    1920             :                   END DO
    1921           0 :                   DO l = 0, maxl
    1922           0 :                      IF (conv(l) == HUGE(i) .AND. g(l) < convpar(l)/10) conv(l) = ishell + ADDSHELL1
    1923             :                   END DO
    1924             :                END IF
    1925           0 :                IF (ishell > conv(maxl) .AND. maxl /= 0) maxl = maxl - 1
    1926           0 :                call timestart("harmonics")
    1927           0 :                call ylm4(maxl, ra, y)
    1928           0 :                call timestop("harmonics")
    1929           0 :                y = CONJG(y)
    1930           0 :                call timestart("kloop")
    1931           0 :                DO ikpt = 1, kpts%nkpt
    1932           0 :                   rdum = kpts%bk(1, ikpt)*ptsh(1, i) + kpts%bk(2, ikpt)*ptsh(2, i) + kpts%bk(3, ikpt)*ptsh(3, i)
    1933           0 :                   cexp = EXP(CMPLX(0.0, 1.0)*2*pi_const*rdum)
    1934           0 :                   lm = 0
    1935           0 :                   DO l = 0, maxl
    1936           0 :                      IF (ishell <= conv(l)) THEN
    1937           0 :                         cdum = cexp*g(l)
    1938           0 :                         DO M = -l, l
    1939           0 :                            lm = lm + 1
    1940           0 :                            shlp(lm, ikpt) = shlp(lm, ikpt) + cdum*y(lm)
    1941             :                         END DO
    1942             :                      ELSE
    1943           0 :                         lm = lm + 2*l + 1
    1944             :                      END IF
    1945             :                   END DO
    1946             :                END DO
    1947           0 :                call timestop("kloop")
    1948             :             END DO
    1949           0 :             structconst(:, ic1, ic2, :) = shlp
    1950             :          END DO
    1951             :       END DO
    1952           0 :       call timestop("realspace sum")
    1953             : 
    1954           0 :       DEALLOCATE (ptsh, radsh)
    1955             : 
    1956           0 :       CALL cpu_TIME(time2)
    1957           0 :       IF (first) WRITE (6, '(A,F7.2)') '  Timing: ', time2 - time1
    1958           0 :       CALL cpu_TIME(time1)
    1959             : 
    1960           0 :       IF (first) WRITE (6, '(/A)') 'Fourier-space sum'
    1961             : 
    1962             :       !
    1963             :       !     Determine reciprocal shells
    1964             :       !
    1965           0 :       call timestart("determince reciproc. shell")
    1966           0 :       CALL getshells(ptsh, nptsh, radsh, nshell, rrad, cell%bmat, first)
    1967           0 :       call timestop("determince reciproc. shell")
    1968             :       ! minimum nonzero reciprocal-shell radius (needed in routines concerning the non-local hartree-fock exchange)
    1969             :       !hybrid%radshmin = radsh(2)
    1970             :       !
    1971             :       !     Fourier-space sum
    1972             :       !
    1973           0 :       call timestart("fourierspace sum")
    1974           0 :       DO ikpt = 1, kpts%nkpt
    1975           0 :          k = kpts%bk(:, ikpt)
    1976           0 :          maxl = MIN(7, hybrid%lexp*2)
    1977           0 :          ishell = 1
    1978           0 :          conv = HUGE(i)
    1979           0 :          DO i = 1, nptsh
    1980           0 :             IF (i > 1) THEN
    1981           0 :                IF (ABS(radsh(i) - radsh(i - 1)) > 1e-10) ishell = ishell + 1
    1982             :             ENDIF
    1983           0 :             ki = ptsh(:, i) + k - NINT(k) ! -nint(...) transforms to Wigner-Seitz cell ( i.e. -0.5 <= x,y,z < 0.5 )
    1984           0 :             ka = MATMUL(ki, cell%bmat)
    1985           0 :             a = SQRT(SUM(ka**2))/scale
    1986           0 :             aa = (1 + a**2)**(-1)
    1987           0 :             IF (ABS(a - a1) > 1e-10) THEN
    1988           0 :                a1 = a
    1989           0 :                IF (a == 0) THEN
    1990           0 :                   g(0) = pref*(-4)
    1991           0 :                   g(1) = 0
    1992             :                ELSE
    1993           0 :                   IF (ishell <= conv(0)) g(0) = pref*aa**4/a**2
    1994           0 :                   IF (ishell <= conv(1)) g(1) = pref*aa**4/a
    1995             :                END IF
    1996           0 :                IF (ishell <= conv(2)) g(2) = pref*aa**5/3
    1997           0 :                IF (ishell <= conv(3)) g(3) = pref*aa**5*a/15
    1998           0 :                IF (ishell <= conv(4)) g(4) = pref*aa**6*a**2/105
    1999           0 :                IF (ishell <= conv(5)) g(5) = pref*aa**6*a**3/945
    2000           0 :                IF (ishell <= conv(6)) g(6) = pref*aa**7*a**4/10395
    2001           0 :                IF (ishell <= conv(7)) g(7) = pref*aa**7*a**5/135135
    2002           0 :                IF (ishell > 1) THEN
    2003           0 :                   DO l = 0, 7
    2004           0 :                      IF (conv(l) == HUGE(i) .AND. g(l) < convpar(l)) conv(l) = ishell + ADDSHELL2
    2005             :                   END DO
    2006             :                END IF
    2007             :             END IF
    2008             : 
    2009           0 :             IF (ishell > conv(maxl) .AND. maxl /= 0) maxl = maxl - 1
    2010           0 :             call timestart("harmonics")
    2011           0 :             call ylm4(maxl, ka, y)
    2012           0 :             call timestop("harmonics")
    2013           0 :             cdum = 1.0
    2014           0 :             lm = 0
    2015           0 :             DO l = 0, maxl
    2016           0 :                IF (ishell <= conv(l)) THEN
    2017           0 :                   DO M = -l, l
    2018           0 :                      lm = lm + 1
    2019           0 :                      y(lm) = CONJG(y(lm))*cdum*g(l)
    2020             :                   END DO
    2021             :                ELSE
    2022           0 :                   y(lm + 1:lm + 2*l + 1) = 0
    2023             :                   lm = lm + 2*l + 1
    2024             :                END IF
    2025           0 :                cdum = cdum*CMPLX(0.0, 1.0)
    2026             :             END DO
    2027           0 :             DO ic2 = 1, atoms%nat
    2028           0 :                DO ic1 = 1, atoms%nat
    2029           0 :                   IF (ic2 /= 1 .AND. ic1 == ic2) CYCLE
    2030           0 :                   cexp = EXP(CMPLX(0.0, 1.0)*2*pi_const*dot_PRODUCT(ki, atoms%taual(:, ic1) - atoms%taual(:, ic2)))
    2031           0 :                   DO lm = 1, (maxl + 1)**2
    2032           0 :                      structconst(lm, ic1, ic2, ikpt) = structconst(lm, ic1, ic2, ikpt) + cexp*y(lm)
    2033             :                   END DO
    2034             :                END DO
    2035             :             END DO
    2036             :          END DO
    2037             :       END DO
    2038           0 :       call timestop("fourierspace sum")
    2039             : 
    2040           0 :       CALL cpu_TIME(time2)
    2041           0 :       IF (first) WRITE (6, '(A,F7.2)') '  Timing: ', time2 - time1
    2042             : 
    2043             :       !
    2044             :       !     Add contribution for l=0 to diagonal elements and rescale structure constants
    2045             :       !
    2046           0 :       structconst(1, 1, 1, :) = structconst(1, 1, 1, :) - 5.0/16/SQRT(4*pi_const)
    2047           0 :       DO i = 2, atoms%nat
    2048           0 :          structconst(:, i, i, :) = structconst(:, 1, 1, :)
    2049             :       END DO
    2050           0 :       DO l = 0, 2*hybrid%lexp
    2051           0 :          structconst(l**2 + 1:(l + 1)**2, :, :, :) = structconst(l**2 + 1:(l + 1)**2, :, :, :)*scale**(l + 1)
    2052             :       END DO
    2053             : 
    2054           0 :       rad = (cell%vol*3/4/pi_const)**(1.0/3) ! Wigner-Seitz radius (rad is recycled)
    2055             : 
    2056             :       !     Calculate accuracy of Gamma-decomposition
    2057           0 :       IF (ALL(kpts%bk == 0)) GOTO 4
    2058           0 :       a = 1e30 ! ikpt = index of shortest non-zero k-point
    2059           0 :       DO i = 2, kpts%nkpt
    2060           0 :          rdum = SUM(MATMUL(kpts%bk(:, i), cell%bmat)**2)
    2061           0 :          IF (rdum < a) THEN
    2062           0 :             ikpt = i
    2063           0 :             a = rdum
    2064             :          END IF
    2065             :       END DO
    2066           0 :       rdum = SQRT(SUM(MATMUL(kpts%bk(:, ikpt), cell%bmat)**2))
    2067           0 :       a = 0
    2068           0 :       DO ic2 = 1, atoms%nat
    2069           0 :          DO ic1 = 1, MAX(1, ic2 - 1)
    2070             :             a = a + ABS(structconst(1, ic1, ic2, ikpt) - &
    2071             :                         (structconst(1, ic1, ic2, 1) + SQRT(4*pi_const)/cell%vol/rdum**2* &
    2072             :                          EXP(-CMPLX(0.0, 1.0)*2*pi_const*dot_PRODUCT( &
    2073           0 :                              kpts%bk(:, ikpt), atoms%taual(:, ic2) - atoms%taual(:, ic1)))))**2
    2074             :          END DO
    2075             :       END DO
    2076           0 :       a = SQRT(a/atoms%nat**2)
    2077           0 :       aa = SQRT(SUM(ABS(structconst(1, :, :, ikpt))**2)/atoms%nat**2)
    2078           0 :       IF (first) WRITE (6, '(/A,F8.5,A,F8.5,A)') 'Accuracy of Gamma-decomposition (structureconstant):', a, ' (abs)', a/aa, ' (rel)'
    2079             : 
    2080           0 : 4     DEALLOCATE (ptsh, radsh)
    2081             : 
    2082           0 :       first = .FALSE.
    2083             : 
    2084           0 :       call timestop("calc struc_const.")
    2085           0 :    END SUBROUTINE structureconstant
    2086             : 
    2087             :    !     -----------------
    2088             : 
    2089             :    !     Determines all shells of the crystal defined by lat and vol with radii smaller than rad.
    2090             :    !     The lattice points (number = nptsh) are stored in ptsh, their corresponding lengths (shell radii) in radsh.
    2091             : 
    2092           0 :    SUBROUTINE getshells(ptsh, nptsh, radsh, nshell, rad, lat, lwrite)
    2093             :       USE m_util, ONLY: rorderpf
    2094             :       IMPLICIT NONE
    2095             :       LOGICAL, INTENT(IN)    :: lwrite
    2096             :       INTEGER, INTENT(OUT)   :: nptsh, nshell
    2097             :       INTEGER, ALLOCATABLE   :: ptsh(:, :)
    2098             :       REAL, ALLOCATABLE   :: radsh(:)
    2099             :       REAL, INTENT(IN)    :: rad, lat(3, 3)
    2100             :       REAL                    :: r(3), rdum
    2101           0 :       INTEGER, ALLOCATABLE   :: pnt(:)
    2102             :       INTEGER                 :: n, i, ix, iy, iz, ok
    2103             :       LOGICAL                 :: found
    2104           0 :       INTEGER, ALLOCATABLE   :: ihelp(:, :)
    2105           0 :       REAL, ALLOCATABLE   :: rhelp(:)
    2106             : 
    2107           0 :       ALLOCATE (ptsh(3, 100000), radsh(100000), stat=ok)
    2108           0 :       IF (ok /= 0) STOP 'getshells: failure allocation ptsh/radsh'
    2109             : 
    2110           0 :       ptsh = 0
    2111           0 :       radsh = 0
    2112             : 
    2113             :       n = 0
    2114             :       i = 0
    2115           0 :       DO
    2116           0 :          found = .FALSE.
    2117           0 :          DO ix = -n, n
    2118           0 :             DO iy = -(n - ABS(ix)), n - ABS(ix)
    2119           0 :                iz = n - ABS(ix) - ABS(iy)
    2120           0 : 1              r = ix*lat(:, 1) + iy*lat(:, 2) + iz*lat(:, 3)
    2121           0 :                rdum = SUM(r**2)
    2122           0 :                IF (rdum < rad**2) THEN
    2123           0 :                   found = .TRUE.
    2124           0 :                   i = i + 1
    2125           0 :                   IF (i > SIZE(radsh)) THEN
    2126           0 :                      ALLOCATE (rhelp(SIZE(radsh)), ihelp(3, SIZE(ptsh, 2)), stat=ok)
    2127           0 :                      IF (ok /= 0) STOP 'getshells: failure allocation rhelp/ihelp'
    2128           0 :                      rhelp = radsh
    2129           0 :                      ihelp = ptsh
    2130           0 :                      DEALLOCATE (radsh, ptsh)
    2131           0 :                      ALLOCATE (radsh(SIZE(rhelp) + 100000), ptsh(3, SIZE(ihelp, 2) + 100000), stat=ok)
    2132           0 :                      IF (ok /= 0) STOP 'getshells: failure re-allocation ptsh/radsh'
    2133           0 :                      radsh(1:SIZE(rhelp)) = rhelp
    2134           0 :                      ptsh(:, 1:SIZE(ihelp, 2)) = ihelp
    2135           0 :                      DEALLOCATE (rhelp, ihelp)
    2136             :                   END IF
    2137           0 :                   ptsh(:, i) = (/ix, iy, iz/)
    2138           0 :                   radsh(i) = SQRT(rdum)
    2139             :                END IF
    2140           0 :                IF (iz > 0) THEN
    2141           0 :                   iz = -iz
    2142           0 :                   GOTO 1
    2143             :                END IF
    2144             :             END DO
    2145             :          END DO
    2146           0 :          IF (.NOT. found) EXIT
    2147           0 :          n = n + 1
    2148             :       END DO
    2149           0 :       nptsh = i
    2150             : 
    2151           0 :       ALLOCATE (pnt(nptsh))
    2152             : 
    2153             :       !reallocate radsh ptsh
    2154           0 :       ALLOCATE (rhelp(nptsh), ihelp(3, nptsh))
    2155           0 :       rhelp = radsh(1:nptsh)
    2156           0 :       ihelp = ptsh(:, 1:nptsh)
    2157           0 :       DEALLOCATE (radsh, ptsh)
    2158           0 :       ALLOCATE (radsh(nptsh), ptsh(3, nptsh))
    2159           0 :       radsh = rhelp
    2160           0 :       ptsh = ihelp
    2161           0 :       DEALLOCATE (rhelp, ihelp)
    2162             : 
    2163           0 :       CALL rorderpf(pnt, radsh, nptsh, MAX(0, INT(LOG(nptsh*0.001)/LOG(2.0))))
    2164           0 :       radsh = radsh(pnt)
    2165           0 :       ptsh = ptsh(:, pnt)
    2166           0 :       nshell = 1
    2167           0 :       DO i = 2, nptsh
    2168           0 :          IF (radsh(i) - radsh(i - 1) > 1e-10) nshell = nshell + 1
    2169             :       END DO
    2170             : 
    2171           0 :       IF (lwrite) &
    2172             :          WRITE (6, '(A,F10.5,A,I7,A,I5,A)') &
    2173           0 :          '  Sphere of radius', rad, ' contains', &
    2174           0 :          nptsh, ' lattice points and', nshell, ' shells.'
    2175             : 
    2176           0 :    END SUBROUTINE getshells
    2177             : 
    2178             :    !     ---------
    2179             : 
    2180             :    !     Returns a list of (k+G) vector lengths in qnrm(1:nqnrm) and the corresponding pointer pqnrm(1:ngpt(ikpt),ikpt)
    2181           0 :    SUBROUTINE getnorm(kpts, gpt, ngpt, pgpt, qnrm, nqnrm, pqnrm, cell)
    2182             :       USE m_types
    2183             :       IMPLICIT NONE
    2184             :       TYPE(t_cell), INTENT(IN)   :: cell
    2185             :       TYPE(t_kpts), INTENT(IN)   :: kpts
    2186             : 
    2187             :       INTEGER, INTENT(IN)  :: ngpt(kpts%nkpt), gpt(3, *), pgpt(:, :)!(dim,kpts%nkpt)
    2188           0 :       REAL, ALLOCATABLE :: qnrm(:), help(:)
    2189             :       INTEGER, INTENT(OUT) :: nqnrm
    2190             :       INTEGER, ALLOCATABLE :: pqnrm(:, :)
    2191             :       INTEGER               :: i, j, ikpt, igpt, igptp
    2192             :       REAL                  :: q(3), qnorm
    2193             : 
    2194           0 :       ALLOCATE (qnrm(MAXVAL(ngpt)*kpts%nkpt), pqnrm(MAXVAL(ngpt), kpts%nkpt))
    2195           0 :       i = 0
    2196           0 :       DO ikpt = 1, kpts%nkpt
    2197           0 :          igptloop: DO igpt = 1, ngpt(ikpt)
    2198           0 :             igptp = pgpt(igpt, ikpt)
    2199           0 :             IF (igptp == 0) STOP 'getnorm: zero pointer (bug?)'
    2200           0 :             q = MATMUL(kpts%bk(:, ikpt) + gpt(:, igptp), cell%bmat)
    2201           0 :             qnorm = SQRT(SUM(q**2))
    2202           0 :             DO j = 1, i
    2203           0 :                IF (ABS(qnrm(j) - qnorm) < 1e-12) THEN
    2204           0 :                   pqnrm(igpt, ikpt) = j
    2205           0 :                   CYCLE igptloop
    2206             :                END IF
    2207             :             END DO
    2208           0 :             i = i + 1
    2209           0 :             qnrm(i) = qnorm
    2210           0 :             pqnrm(igpt, ikpt) = i
    2211             :          END DO igptloop
    2212             :       END DO
    2213           0 :       nqnrm = i
    2214             : 
    2215           0 :       ALLOCATE (help(nqnrm))
    2216           0 :       help(1:nqnrm) = qnrm(1:nqnrm)
    2217           0 :       DEALLOCATE (qnrm)
    2218           0 :       ALLOCATE (qnrm(1:nqnrm))
    2219           0 :       qnrm = help
    2220             : 
    2221           0 :    END SUBROUTINE getnorm
    2222             : 
    2223           0 :    FUNCTION sphbessel_integral(atoms, itype, qnrm, nqnrm, iqnrm1, iqnrm2, l, hybrid, &
    2224           0 :                                sphbes0, l_warnin, l_warnout)
    2225             : 
    2226             :       USE m_types
    2227             :       IMPLICIT NONE
    2228             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
    2229             :       TYPE(t_atoms), INTENT(IN)   :: atoms
    2230             : 
    2231             :       INTEGER, INTENT(IN)  :: itype, nqnrm, iqnrm1, iqnrm2, l
    2232             :       REAL, INTENT(IN)  :: qnrm(nqnrm), sphbes0(-1:hybrid%lexp + 2, atoms%ntype, nqnrm)
    2233             :       LOGICAL, INTENT(IN), OPTIONAL  ::  l_warnin
    2234             :       LOGICAL, INTENT(OUT), OPTIONAL  ::  l_warnout
    2235             :       REAL                  :: sphbessel_integral
    2236             :       REAL                  :: q1, q2, dq, s, sb01, sb11, sb21, sb31, sb02, sb12
    2237             :       REAL                  :: sb22, sb32, a1, a2, da, b1, b2, db, c1, c2, dc, r1, r2
    2238             :       LOGICAL               :: l_warn, l_warned
    2239             : 
    2240           0 :       IF (PRESENT(l_warnin)) THEN
    2241           0 :          l_warn = l_warnin
    2242             :       ELSE
    2243             :          l_warn = .TRUE.
    2244             :       END IF
    2245           0 :       l_warned = .FALSE.
    2246             : 
    2247           0 :       q1 = qnrm(iqnrm1)
    2248           0 :       q2 = qnrm(iqnrm2)
    2249           0 :       s = atoms%rmt(itype)
    2250           0 :       IF (q1 == 0 .AND. q2 == 0) THEN
    2251           0 :          IF (l > 0) THEN
    2252             :             sphbessel_integral = 0
    2253             :          ELSE
    2254           0 :             sphbessel_integral = 2*s**5/15
    2255             :          ENDIF
    2256           0 :       ELSE IF (q1 == 0 .OR. q2 == 0) THEN
    2257           0 :          IF (l > 0) THEN
    2258             :             sphbessel_integral = 0
    2259           0 :          ELSE IF (q1 == 0) THEN
    2260             :             sphbessel_integral = s**3/(3*q2**2)*(q2*s*sphbes0(1, itype, iqnrm2) &
    2261           0 :                                                  + sphbes0(2, itype, iqnrm2))
    2262             :          ELSE
    2263             :             sphbessel_integral = s**3/(3*q1**2)*(q1*s*sphbes0(1, itype, iqnrm1) &
    2264           0 :                                                  + sphbes0(2, itype, iqnrm1))
    2265             :          ENDIF
    2266           0 :       ELSE IF (q1 == q2) THEN
    2267             :          sphbessel_integral = s**3/(2*q1**2)*((2*l + 3)*sphbes0(l + 1, itype, iqnrm1)**2 - &
    2268           0 :                                               (2*l + 1)*sphbes0(l, itype, iqnrm1)*sphbes0(l + 2, itype, iqnrm1))
    2269             :       ELSE ! We use either if two fromulas that are stable for high and small q1/q2 respectively
    2270           0 :          sb01 = sphbes0(l - 1, itype, iqnrm1)
    2271           0 :          sb11 = sphbes0(l, itype, iqnrm1)
    2272           0 :          sb21 = sphbes0(l + 1, itype, iqnrm1)
    2273           0 :          sb31 = sphbes0(l + 2, itype, iqnrm1)
    2274           0 :          sb02 = sphbes0(l - 1, itype, iqnrm2)
    2275           0 :          sb12 = sphbes0(l, itype, iqnrm2)
    2276           0 :          sb22 = sphbes0(l + 1, itype, iqnrm2)
    2277           0 :          sb32 = sphbes0(l + 2, itype, iqnrm2)
    2278           0 :          dq = q1**2 - q2**2
    2279           0 :          a1 = q2/q1*sb21*sb02
    2280           0 :          a2 = q1/q2*sb22*sb01
    2281           0 :          da = a1 - a2
    2282           0 :          b1 = sb31*sb12
    2283           0 :          b2 = sb32*sb11
    2284           0 :          db = b1 - b2
    2285           0 :          c1 = sb21*sb22/(q1*q2)
    2286           0 :          c2 = db/dq*(2*l + 1)/(2*l + 3)
    2287           0 :          dc = c1 + c2
    2288           0 :          r1 = ABS(da/a1)
    2289           0 :          r2 = MIN(ABS(db/b1), ABS(dc/c1))
    2290             :          ! Ensure numerical stability. If both formulas are not sufficiently stable, the program stops.
    2291           0 :          IF (r1 > r2) THEN
    2292           0 :             IF (r1 < 1e-6 .AND. l_warn) THEN
    2293           0 :                WRITE (6, '(A,E12.5,A,E12.5,A)') 'sphbessel_integral: Warning! Formula One possibly unstable. Ratios:', &
    2294           0 :                   r1, '(', r2, ')'
    2295           0 :                WRITE (6, '(A,2F15.10,I4)') '                    Current qnorms and atom type:', q1, q2, itype
    2296           0 :                l_warned = .TRUE.
    2297             :             END IF
    2298           0 :             sphbessel_integral = s**3/dq*da
    2299             :          ELSE
    2300           0 :             IF (r2 < 1e-6 .AND. l_warn) THEN
    2301           0 :                WRITE (6, '(A,E13.5,A,E13.5,A)') 'sphbessel_integral: Warning! Formula Two possibly unstable. Ratios:', &
    2302           0 :                   r2, '(', r1, ')'
    2303           0 :                WRITE (6, '(A,2F15.10,I4)') '                    Current qnorms and atom type:', &
    2304           0 :                   q1, q2, itype
    2305           0 :                l_warned = .TRUE.
    2306             :             END IF
    2307           0 :             sphbessel_integral = s**3*dc
    2308             :          END IF
    2309             :       END IF
    2310             : 
    2311           0 :       IF (PRESENT(l_warnout)) l_warnout = l_warned
    2312             : 
    2313           0 :    END FUNCTION sphbessel_integral
    2314             : 
    2315             :    !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    2316             : 
    2317             : END MODULE m_coulombmatrix

Generated by: LCOV version 1.13