subvxc.F90 22.8 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

Daniel Wortmann's avatar
Daniel Wortmann committed
7
MODULE m_subvxc
8

Daniel Wortmann's avatar
Daniel Wortmann committed
9 10
CONTAINS

Matthias Redies's avatar
Matthias Redies committed
11 12
   SUBROUTINE subvxc(lapw, bk, DIMENSION, input, jsp, vr0, atoms, usdus, hybrid, el, ello, sym, &
                     cell, sphhar, stars, xcpot, mpi, oneD, hmat, vx)
13

14 15
      USE m_types
      USE m_judft
Matthias Redies's avatar
Matthias Redies committed
16
      USE m_intgr, ONLY: intgr3
17
      USE m_constants
Matthias Redies's avatar
Matthias Redies committed
18
      USE m_gaunt, ONLY: gaunt1
19 20 21 22 23 24 25 26
      USE m_wrapper
      USE m_loddop
      USE m_radflo
      USE m_radfun
      USE m_abcof3

      IMPLICIT NONE

Matthias Redies's avatar
Matthias Redies committed
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
      CLASS(t_xcpot), INTENT(IN)    :: xcpot
      TYPE(t_mpi), INTENT(IN)    :: mpi
      TYPE(t_dimension), INTENT(IN)    :: dimension
      TYPE(t_oneD), INTENT(IN)    :: oneD
      TYPE(t_hybrid), INTENT(INOUT) :: hybrid
      TYPE(t_input), INTENT(IN)    :: input
      TYPE(t_sym), INTENT(IN)    :: sym
      TYPE(t_stars), INTENT(IN)    :: stars
      TYPE(t_cell), INTENT(IN)    :: cell
      TYPE(t_sphhar), INTENT(IN)    :: sphhar
      TYPE(t_atoms), INTENT(IN)    :: atoms
      TYPE(t_lapw), INTENT(IN)    :: lapw
      TYPE(t_usdus), INTENT(INOUT) :: usdus
      TYPE(t_potden), INTENT(IN)    :: vx
      TYPE(t_mat), INTENT(INOUT) :: hmat
42 43

      ! Scalar Arguments
Matthias Redies's avatar
Matthias Redies committed
44
      INTEGER, INTENT(IN) :: jsp
45 46

      ! Array Arguments
Matthias Redies's avatar
Matthias Redies committed
47 48 49 50
      REAL, INTENT(IN) :: vr0(atoms%jmtd, atoms%ntype, input%jspins)               ! just for radial functions
      REAL, INTENT(IN) :: el(0:atoms%lmaxd, atoms%ntype, input%jspins)
      REAL, INTENT(IN) :: ello(atoms%nlod, atoms%ntype, input%jspins)
      REAL, INTENT(IN) :: bk(3)
51 52

      ! Local Scalars
Matthias Redies's avatar
Matthias Redies committed
53 54 55
      INTEGER               ::  ic, indx, m, ig1, ig2, n, nn
      INTEGER               ::  nlharm, nnbas, typsym, lm
      INTEGER               ::  noded, nodeu
56
      INTEGER               ::  nbasf0
Matthias Redies's avatar
Matthias Redies committed
57 58 59 60 61 62
      INTEGER               ::  i, j, l, ll, l1, l2, m1, m2, j1, j2
      INTEGER               ::  ok, p1, p2, lh, mh, pp1, pp2
      INTEGER               ::  igrid, itype, ilharm, istar
      INTEGER               ::  ineq, iatom, ilo, ilop, ieq, icentry
      INTEGER               ::  ikvecat, ikvecprevat, invsfct, ikvec, ikvecp
      INTEGER               ::  lp, mp, pp
63 64
      REAL                  ::  a_ex
      REAL                  ::  wronk
Matthias Redies's avatar
Matthias Redies committed
65
      COMPLEX               ::  rc, rr
66 67 68

      ! Local Arrays
      INTEGER               ::  gg(3)
Matthias Redies's avatar
Matthias Redies committed
69
      INTEGER               ::  pointer_lo(atoms%nlod, atoms%ntype)
70

Matthias Redies's avatar
Matthias Redies committed
71
      REAL                  ::  integ(0:sphhar%nlhd, hybrid%maxindx, 0:atoms%lmaxd, hybrid%maxindx, 0:atoms%lmaxd)
72
      REAL                  ::  grid(atoms%jmtd)
Matthias Redies's avatar
Matthias Redies committed
73 74 75 76 77
      REAL                  ::  vr(atoms%jmtd, 0:sphhar%nlhd)
      REAL                  ::  f(atoms%jmtd, 2, 0:atoms%lmaxd), g(atoms%jmtd, 2, 0:atoms%lmaxd)
      REAL                  ::  flo(atoms%jmtd, 2, atoms%nlod)
      REAL                  ::  uuilon(atoms%nlod, atoms%ntype), duilon(atoms%nlod, atoms%ntype)
      REAL                  ::  ulouilopn(atoms%nlod, atoms%nlod, atoms%ntype)
78

Matthias Redies's avatar
Matthias Redies committed
79 80
      REAL                  ::  bas1(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype)
      REAL                  ::  bas2(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype)
81 82

      COMPLEX               ::  vpw(stars%ng3)
Matthias Redies's avatar
Matthias Redies committed
83 84 85 86 87
      COMPLEX               ::  vxc(hmat%matsize1*(hmat%matsize1 + 1)/2)
      COMPLEX               ::  vrmat(hybrid%maxlmindx, hybrid%maxlmindx)
      COMPLEX               ::  carr(hybrid%maxlmindx, DIMENSION%nvd), carr1(DIMENSION%nvd, DIMENSION%nvd)
      COMPLEX, ALLOCATABLE  ::  ahlp(:, :, :), bhlp(:, :, :)
      COMPLEX, ALLOCATABLE  ::  bascof(:, :, :)
88
#ifndef CPP_OLDINTEL
Matthias Redies's avatar
Matthias Redies committed
89
      COMPLEX               ::  bascof_lo(3, -atoms%llod:atoms%llod, 4*atoms%llod + 2, atoms%nlod, atoms%nat)
90
#endif
91 92

      CALL timestart("subvxc")
Matthias Redies's avatar
Matthias Redies committed
93
      vxc = 0
94 95

      ! Calculate radial functions
Matthias Redies's avatar
Matthias Redies committed
96 97
      hybrid%nindx = 2
      DO itype = 1, atoms%ntype
98 99

         ! Generate the radial basis-functions for each l
Matthias Redies's avatar
Matthias Redies committed
100 101 102 103 104 105 106 107 108
         WRITE (6, '(a,i3,a)') new_LINE('n')//new_LINE('n')//' wavefunction parameters for atom type', itype, ':'
         WRITE (6, '(31x,a,32x,a)') 'radial function', 'energy derivative'
         WRITE (6, '(a)') '  l    energy            value        '// &
            'derivative    nodes          value        derivative    nodes       norm        wronskian'
         DO l = 0, atoms%lmax(itype)
            CALL radfun(l, itype, jsp, el(l, itype, jsp), vr0(1, itype, jsp), atoms, f(1, 1, l), g(1, 1, l), usdus, nodeu, noded, wronk)
            WRITE (6, FMT=8010) l, el(l, itype, jsp), usdus%us(l, itype, jsp), &
               usdus%dus(l, itype, jsp), nodeu, usdus%uds(l, itype, jsp), usdus%duds(l, itype, jsp), noded, &
               usdus%ddn(l, itype, jsp), wronk
109
         END DO
Matthias Redies's avatar
Matthias Redies committed
110
8010     FORMAT(i3, f10.5, 2(5x, 1p, 2e16.7, i5), 1p, 2e16.7)
111

Matthias Redies's avatar
Matthias Redies committed
112 113 114 115
         bas1(:, 1, :, itype) = f(:, 1, :)
         bas1(:, 2, :, itype) = g(:, 1, :)
         bas2(:, 1, :, itype) = f(:, 2, :)
         bas2(:, 2, :, itype) = g(:, 2, :)
116 117

         ! Generate the extra radial basis-functions for the local orbitals, if there are any.
Matthias Redies's avatar
Matthias Redies committed
118 119 120
         IF (atoms%nlo(itype) >= 1) THEN
            CALL radflo(atoms, itype, jsp, ello(1, 1, jsp), vr0(1, itype, jsp), f, g, mpi, &
                        usdus, uuilon, duilon, ulouilopn, flo, .TRUE.)
121 122

            DO i = 1, atoms%nlo(itype)
Matthias Redies's avatar
Matthias Redies committed
123 124 125 126
               hybrid%nindx(atoms%llo(i, itype), itype) = hybrid%nindx(atoms%llo(i, itype), itype) + 1
               pointer_lo(i, itype) = hybrid%nindx(atoms%llo(i, itype), itype)
               bas1(:, hybrid%nindx(atoms%llo(i, itype), itype), atoms%llo(i, itype), itype) = flo(:, 1, i)
               bas2(:, hybrid%nindx(atoms%llo(i, itype), itype), atoms%llo(i, itype), itype) = flo(:, 2, i)
127 128 129 130 131 132 133
            END DO
         END IF
      END DO

      ! Compute APW coefficients

      ! Calculate bascof
Matthias Redies's avatar
Matthias Redies committed
134 135
      ALLOCATE (ahlp(DIMENSION%nvd, 0:DIMENSION%lmd, atoms%nat), bhlp(DIMENSION%nvd, 0:DIMENSION%lmd, atoms%nat), stat=ok)
      IF (ok /= 0) STOP 'subvxc: error in allocation of ahlp/bhlp'
136
#ifndef CPP_OLDINTEL
Matthias Redies's avatar
Matthias Redies committed
137
      CALL abcof3(input, atoms, sym, jsp, cell, bk, lapw, usdus, oneD, ahlp, bhlp, bascof_lo)
138
#endif
Matthias Redies's avatar
Matthias Redies committed
139 140
      ALLOCATE (bascof(DIMENSION%nvd, 2*(DIMENSION%lmd + 1), atoms%nat), stat=ok)
      IF (ok /= 0) STOP 'subvxc: error in allocation of bascof'
141 142

      bascof = 0
Matthias Redies's avatar
Matthias Redies committed
143
      ic = 0
144 145 146 147 148
      DO itype = 1, atoms%ntype
         DO ieq = 1, atoms%neq(itype)
            ic = ic + 1
            indx = 0
            DO l = 0, atoms%lmax(itype)
Matthias Redies's avatar
Matthias Redies committed
149
               ll = l*(l + 1)
150
               DO M = -l, l
Matthias Redies's avatar
Matthias Redies committed
151
                  lm = ll + M
152 153
                  DO i = 1, 2
                     indx = indx + 1
Matthias Redies's avatar
Matthias Redies committed
154 155 156 157
                     IF (i == 1) THEN
                        bascof(:, indx, ic) = ahlp(:, lm, ic)
                     ELSE IF (i == 2) THEN
                        bascof(:, indx, ic) = bhlp(:, lm, ic)
158 159 160 161 162 163 164
                     END IF
                  END DO
               END DO
            END DO
         END DO
      END DO

Matthias Redies's avatar
Matthias Redies committed
165
      DEALLOCATE (ahlp, bhlp)
166 167 168 169 170

      ! Loop over atom types
      iatom = 0
      DO itype = 1, atoms%ntype

Matthias Redies's avatar
Matthias Redies committed
171
         typsym = atoms%ntypsy(SUM(atoms%neq(:itype - 1)) + 1)
172 173 174
         nlharm = sphhar%nlh(typsym)

         ! Calculate vxc = vtot - vcoul
Matthias Redies's avatar
Matthias Redies committed
175 176 177 178 179
         DO l = 0, nlharm
            DO i = 1, atoms%jri(itype)
               IF (l == 0) THEN
                  ! vr(i,0)= vrtot(i,0,itype)*sfp/rmsh(i,itype) - vrcou(i,0,itype,jsp)
                  vr(i, 0) = vx%mt(i, 0, itype, jsp)! * sfp_const / atoms%rmsh(i,itype)
180 181
               ELSE ! vxc = vtot - vcoul
                  ! vr(i,l)= vrtot(i,l,itype)-vrcou(i,l,itype,jsp)
Matthias Redies's avatar
Matthias Redies committed
182
                  vr(i, l) = vx%mt(i, l, itype, jsp)
183 184 185 186 187 188 189 190 191 192 193 194
               END IF
            END DO
         END DO

         ! Calculate MT contribution to vxc matrix elements
         ! Precompute auxiliary radial integrals
         DO ilharm = 0, nlharm
            i = 0
            DO l1 = 0, atoms%lmax(itype)
               DO p1 = 1, 2
                  i = i + 1
                  j = 0
Matthias Redies's avatar
Matthias Redies committed
195
                  DO l2 = 0, atoms%lmax(itype)
196 197
                     DO p2 = 1, 2
                        j = j + 1
Matthias Redies's avatar
Matthias Redies committed
198
                        IF (j <= i) THEN
199
                           DO igrid = 1, atoms%jri(itype)
Matthias Redies's avatar
Matthias Redies committed
200 201
                              grid(igrid) = vr(igrid, ilharm)*(bas1(igrid, p1, l1, itype)*bas1(igrid, p2, l2, itype) + &
                                                               bas2(igrid, p1, l1, itype)*bas2(igrid, p2, l2, itype))
202
                           END DO
Matthias Redies's avatar
Matthias Redies committed
203 204
                           CALL intgr3(grid, atoms%rmsh(:, itype), atoms%dx(itype), atoms%jri(itype), integ(ilharm, p1, l1, p2, l2))
                           integ(ilharm, p2, l2, p1, l1) = integ(ilharm, p1, l1, p2, l2)
205 206 207 208 209 210 211 212
                        END IF
                     END DO
                  END DO
               END DO
            END DO
         END DO

         ! Calculate muffin tin contribution to vxc matrix
Matthias Redies's avatar
Matthias Redies committed
213
         vrmat = 0
214

Matthias Redies's avatar
Matthias Redies committed
215
         j1 = 0
216
         DO l1 = 0, atoms%lmax(itype) ! loop: left basis function
Matthias Redies's avatar
Matthias Redies committed
217
            DO m1 = -l1, l1
218 219 220 221 222 223 224 225 226
               DO p1 = 1, 2
                  j1 = j1 + 1
                  j2 = 0
                  DO l2 = 0, atoms%lmax(itype) ! loop: right basis function
                     DO m2 = -l2, l2
                        DO p2 = 1, 2
                           j2 = j2 + 1
                           rr = 0
                           DO ilharm = 0, nlharm ! loop: lattice harmonics of vxc
Matthias Redies's avatar
Matthias Redies committed
227 228 229 230 231
                              l = sphhar%llh(ilharm, typsym)
                              DO i = 1, sphhar%nmem(ilharm, typsym)
                                 M = sphhar%mlh(i, ilharm, typsym)
                                 rc = sphhar%clnu(i, ilharm, typsym)*gaunt1(l1, l, l2, m1, M, m2, atoms%lmaxd)
                                 rr = rr + integ(ilharm, p1, l1, p2, l2)*rc
232 233
                              END DO
                           END DO
Matthias Redies's avatar
Matthias Redies committed
234 235
                           rc = CMPLX(0, 1)**(l2 - l1) ! adjusts to a/b/ccof-scaling
                           vrmat(j1, j2) = rr*rc
236 237 238 239 240 241 242 243 244
                        END DO
                     END DO
                  END DO
               END DO
            END DO
         END DO
         nnbas = j1

         ! Project on bascof
Matthias Redies's avatar
Matthias Redies committed
245 246 247 248 249
         DO ineq = 1, atoms%neq(itype)
            iatom = iatom + 1
            carr(:nnbas, :lapw%nv(jsp)) = CONJG(MATMUL(vrmat(:nnbas, :nnbas), &
                                                       TRANSPOSE(bascof(:lapw%nv(jsp), :nnbas, iatom))))
            carr1(:lapw%nv(jsp), :lapw%nv(jsp)) = MATMUL(bascof(:lapw%nv(jsp), :nnbas, iatom), carr(:nnbas, :lapw%nv(jsp)))
250
            ic = 0
Matthias Redies's avatar
Matthias Redies committed
251
            DO j = 1, lapw%nv(jsp)
252 253 254
               ! carr(:nnbas) =  matmul(vrmat(:nnbas,:nnbas),bascof(j,:nnbas,iatom) )
               DO i = 1, j
                  ic = ic + 1
Matthias Redies's avatar
Matthias Redies committed
255
                  vxc(ic) = vxc(ic) + carr1(i, j)
256 257 258 259 260 261 262
                  ! vxc(ic) = vxc(ic) + conjg(dotprod ( bascof(i,:nnbas,iatom),carr(:nnbas) ))
               END DO
            END DO
         END DO
      END DO ! End loop over atom types

      ! Calculate plane wave contribution
Matthias Redies's avatar
Matthias Redies committed
263 264 265
      DO i = 1, stars%ng3
         vpw(i) = vx%pw_w(i, jsp)
         ! vpw(i)=vpwtot(i)-vpwcou(i,jsp)
266 267 268 269 270 271 272 273
      END DO

      ! Calculate vxc-matrix,  left basis function (ig1)
      !                        right basis function (ig2)
      ic = 0
      DO ig1 = 1, lapw%nv(jsp)
         DO ig2 = 1, ig1
            ic = ic + 1
Matthias Redies's avatar
Matthias Redies committed
274 275 276 277 278 279
            gg(1) = lapw%k1(ig1, jsp) - lapw%k1(ig2, jsp)
            gg(2) = lapw%k2(ig1, jsp) - lapw%k2(ig2, jsp)
            gg(3) = lapw%k3(ig1, jsp) - lapw%k3(ig2, jsp)
            istar = stars%ig(gg(1), gg(2), gg(3))
            IF (istar /= 0) THEN
               vxc(ic) = vxc(ic) + stars%rgphs(gg(1), gg(2), gg(3))*vpw(istar)
280
            ELSE
Matthias Redies's avatar
Matthias Redies committed
281 282 283 284
               IF (mpi%irank == 0) THEN
                  WRITE (6, '(A,/6I5)') 'Warning: Gi-Gj not in any star:', &
                     lapw%k1(ig1, jsp), lapw%k2(ig1, jsp), lapw%k3(ig1, jsp), &
                     lapw%k1(ig2, jsp), lapw%k2(ig2, jsp), lapw%k3(ig2, jsp)
285 286 287 288 289 290
               END IF
            END IF
         END DO
      END DO

      ! Calculate local orbital contribution
Matthias Redies's avatar
Matthias Redies committed
291
      IF (ANY(atoms%nlo /= 0)) THEN
292

Matthias Redies's avatar
Matthias Redies committed
293
         nbasf0 = lapw%nv(jsp)*(lapw%nv(jsp) + 1)/2  ! number of pure APW contributions
294 295 296 297 298 299 300
         icentry = nbasf0                          ! icentry counts the entry in the matrix vxc
         iatom = 0
         ikvecat = 0
         ikvecprevat = 0

         DO itype = 1, atoms%ntype

Matthias Redies's avatar
Matthias Redies committed
301
            typsym = atoms%ntypsy(SUM(atoms%neq(:itype - 1)) + 1)
302 303 304 305 306
            nlharm = sphhar%nlh(typsym)

            ! Calculate vxc = vtot - vcoul
            DO l = 0, nlharm
               DO i = 1, atoms%jri(itype)
Matthias Redies's avatar
Matthias Redies committed
307
                  IF (l == 0) THEN
308
                     ! vr(i,0)= vrtot(i,0,itype)*sfp/rmsh(i,itype) -  vrcou(i,0,itype,jsp)
Matthias Redies's avatar
Matthias Redies committed
309
                     vr(i, 0) = vx%mt(i, 0, itype, jsp)*sfp_const/atoms%rmsh(i, itype)
310
                  ELSE ! vxc = vtot - vcoul
Matthias Redies's avatar
Matthias Redies committed
311
                     vr(i, l) = vx%mt(i, l, itype, jsp)                    !
312 313 314 315 316 317 318 319 320
                     ! vr(i,l)=  vrtot(i,l,itype)-vrcou(i,l,itype,jsp)
                  END IF
               END DO
            END DO

            ! Precompute auxiliary radial integrals
            DO ilharm = 0, nlharm
               i = 0
               DO l1 = 0, atoms%lmax(itype)
Matthias Redies's avatar
Matthias Redies committed
321
                  DO p1 = 1, hybrid%nindx(l1, itype)
322 323
                     i = i + 1
                     j = 0
Matthias Redies's avatar
Matthias Redies committed
324 325
                     DO l2 = 0, atoms%lmax(itype)
                        DO p2 = 1, hybrid%nindx(l2, itype)
326
                           j = j + 1
Matthias Redies's avatar
Matthias Redies committed
327
                           IF (j <= i) THEN
328
                              DO igrid = 1, atoms%jri(itype)
Matthias Redies's avatar
Matthias Redies committed
329 330
                                 grid(igrid) = vr(igrid, ilharm)*(bas1(igrid, p1, l1, itype)*bas1(igrid, p2, l2, itype) + &
                                                                  bas2(igrid, p1, l1, itype)*bas2(igrid, p2, l2, itype))
331
                              END DO
Matthias Redies's avatar
Matthias Redies committed
332 333
                              CALL intgr3(grid, atoms%rmsh(:, itype), atoms%dx(itype), atoms%jri(itype), integ(ilharm, p1, l1, p2, l2))
                              integ(ilharm, p2, l2, p1, l1) = integ(ilharm, p1, l1, p2, l2)
334 335 336 337 338 339 340 341 342
                           END IF
                        END DO
                     END DO
                  END DO
               END DO
            END DO

            DO ieq = 1, atoms%neq(itype)
               iatom = iatom + 1
Matthias Redies's avatar
Matthias Redies committed
343
               IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN
344

Matthias Redies's avatar
Matthias Redies committed
345 346
                  IF (atoms%invsat(iatom) == 0) invsfct = 1
                  IF (atoms%invsat(iatom) == 1) invsfct = 2
347 348

                  DO ilo = 1, atoms%nlo(itype)
349
#ifdef CPP_OLDINTEL
Matthias Redies's avatar
Matthias Redies committed
350
                     CALL judft_error("no LOs & hybrid with old intel compiler!", calledby="subvxc.F90")
351
#else
Matthias Redies's avatar
Matthias Redies committed
352 353
                     l1 = atoms%llo(ilo, itype)
                     DO ikvec = 1, invsfct*(2*l1 + 1)
354 355
                        DO m1 = -l1, l1
                           DO p1 = 1, 3
Matthias Redies's avatar
Matthias Redies committed
356 357
                              IF (p1 == 3) THEN
                                 pp1 = pointer_lo(ilo, itype)
358 359 360 361
                              ELSE
                                 pp1 = p1
                              END IF

Matthias Redies's avatar
Matthias Redies committed
362
                              IF (hybrid%nindx(l1, itype) <= 2) STOP 'subvxc: error hybrid%nindx'
363 364 365 366 367 368 369 370 371 372 373

                              lm = 0

                              !loop over APW
                              DO l2 = 0, atoms%lmax(itype)
                                 DO m2 = -l2, l2
                                    DO p2 = 1, 2
                                       lm = lm + 1

                                       rr = 0
                                       DO ilharm = 0, nlharm
Matthias Redies's avatar
Matthias Redies committed
374 375 376 377 378
                                          lh = sphhar%llh(ilharm, typsym)
                                          DO i = 1, sphhar%nmem(ilharm, typsym)
                                             mh = sphhar%mlh(i, ilharm, typsym)
                                             rc = sphhar%clnu(i, ilharm, typsym)*gaunt1(l1, lh, l2, m1, mh, m2, atoms%lmaxd)
                                             rr = rr + integ(ilharm, p2, l2, pp1, l1)*rc
379 380 381
                                          END DO
                                       END DO

Matthias Redies's avatar
Matthias Redies committed
382
                                       rc = CMPLX(0.0, 1.0)**(l2 - l1) ! adjusts to a/b/ccof-scaling
383 384 385 386 387

                                       ! ic counts the entry in vxc
                                       ic = icentry
                                       DO i = 1, lapw%nv(jsp)
                                          ic = ic + 1
388
                                          IF (hmat%l_real) THEN
Matthias Redies's avatar
Matthias Redies committed
389 390
                                             vxc(ic) = vxc(ic) + invsfct*REAL(rr*rc*bascof(i, lm, iatom)* &
                                                                              CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom)))
391
                                          ELSE
Matthias Redies's avatar
Matthias Redies committed
392 393
                                             vxc(ic) = vxc(ic) + rr*rc*bascof(i, lm, iatom)* &
                                                       CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom))
394 395 396 397 398 399 400
                                          END IF
                                       END DO
                                    END DO  !p2
                                 END DO  ! m2
                              END DO ! l2 ->  loop over APW

                              ! calcualte matrix-elements with local orbitals at the same atom
Matthias Redies's avatar
Matthias Redies committed
401
                              IF (ic /= icentry + lapw%nv(jsp)) STOP 'subvxc: error counting ic'
402 403 404

                              ic = ic + ikvecprevat

Matthias Redies's avatar
Matthias Redies committed
405 406 407
                              DO ilop = 1, ilo - 1
                                 lp = atoms%llo(ilop, itype)
                                 DO ikvecp = 1, invsfct*(2*lp + 1)
408 409 410
                                    ic = ic + 1
                                    DO mp = -lp, lp
                                       DO pp = 1, 3
Matthias Redies's avatar
Matthias Redies committed
411 412
                                          IF (pp == 3) THEN
                                             pp2 = pointer_lo(ilop, itype)
413 414 415 416 417 418
                                          ELSE
                                             pp2 = pp
                                          END IF

                                          rr = 0
                                          DO ilharm = 0, nlharm
Matthias Redies's avatar
Matthias Redies committed
419 420 421 422 423
                                             lh = sphhar%llh(ilharm, typsym)
                                             DO i = 1, sphhar%nmem(ilharm, typsym)
                                                mh = sphhar%mlh(i, ilharm, typsym)
                                                rc = sphhar%clnu(i, ilharm, typsym)*gaunt1(l1, lh, lp, m1, mh, mp, atoms%lmaxd)
                                                rr = rr + integ(ilharm, pp2, lp, pp1, l1)*rc
424 425 426
                                             END DO
                                          END DO

Matthias Redies's avatar
Matthias Redies committed
427
                                          rc = CMPLX(0.0, 1.0)**(lp - l1) ! adjusts to a/b/ccof-scaling
428

429
                                          IF (hmat%l_real) THEN
Matthias Redies's avatar
Matthias Redies committed
430 431
                                             vxc(ic) = vxc(ic) + invsfct*REAL(rr*rc*bascof_lo(pp, mp, ikvecp, ilop, iatom)* &
                                                                              CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom)))
432
                                          ELSE
Matthias Redies's avatar
Matthias Redies committed
433 434
                                             vxc(ic) = vxc(ic) + rr*rc*bascof_lo(pp, mp, ikvecp, ilop, iatom)* &
                                                       CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom))
435 436 437 438 439 440 441 442 443 444 445
                                          END IF
                                       END DO ! pp
                                    END DO ! mp

                                 END DO !ikvecp
                              END DO ! ilop

                              ! calculate matrix-elements of one local orbital with itself
                              DO ikvecp = 1, ikvec
                                 ic = ic + 1

Matthias Redies's avatar
Matthias Redies committed
446
                                 lp = l1
447 448 449
                                 ilop = ilo
                                 DO mp = -lp, lp
                                    DO pp = 1, 3
Matthias Redies's avatar
Matthias Redies committed
450 451
                                       IF (pp == 3) THEN
                                          pp2 = pointer_lo(ilop, itype)
452 453 454 455 456 457
                                       ELSE
                                          pp2 = pp
                                       END IF

                                       rr = 0
                                       DO ilharm = 0, nlharm
Matthias Redies's avatar
Matthias Redies committed
458 459 460 461 462
                                          lh = sphhar%llh(ilharm, typsym)
                                          DO i = 1, sphhar%nmem(ilharm, typsym)
                                             mh = sphhar%mlh(i, ilharm, typsym)
                                             rc = sphhar%clnu(i, ilharm, typsym)*gaunt1(l1, lh, lp, m1, mh, mp, atoms%lmaxd)
                                             rr = rr + integ(ilharm, pp2, lp, pp1, l1)*rc
463 464 465
                                          END DO
                                       END DO

Matthias Redies's avatar
Matthias Redies committed
466
                                       rc = CMPLX(0.0, 1.0)**(lp - l1) ! adjusts to a/b/ccof-scaling
467

468
                                       IF (hmat%l_real) THEN
Matthias Redies's avatar
Matthias Redies committed
469 470
                                          vxc(ic) = vxc(ic) + invsfct*REAL(rr*rc*bascof_lo(pp, mp, ikvecp, ilop, iatom)* &
                                                                           CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom)))
471
                                       ELSE
Matthias Redies's avatar
Matthias Redies committed
472 473
                                          vxc(ic) = vxc(ic) + rr*rc*bascof_lo(pp, mp, ikvecp, ilop, iatom)* &
                                                    CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom))
474 475 476 477 478 479
                                       END IF
                                    END DO ! pp
                                 END DO ! mp
                              END DO ! ikvecp
                           END DO  ! p1
                        END DO  ! m1
Matthias Redies's avatar
Matthias Redies committed
480
                        icentry = ic
481
                     END DO !ikvec
Matthias Redies's avatar
Matthias Redies committed
482
                     ikvecat = ikvecat + invsfct*(2*l1 + 1)
483
#endif
484 485 486 487 488 489 490 491 492 493
                  END DO  ! ilo
                  ikvecprevat = ikvecprevat + ikvecat
                  ikvecat = 0
               END IF ! atoms%invsat(iatom)
            END DO ! ieq
         END DO !itype
      END IF ! if any atoms%llo

      !initialize weighting factor
      a_ex = xcpot%get_exchange_weight()
494 495

      i = 0
Matthias Redies's avatar
Matthias Redies committed
496
      DO n = 1, hmat%matsize1
497 498 499
         DO nn = 1, n
            i = i + 1
            IF (hmat%l_real) THEN
Matthias Redies's avatar
Matthias Redies committed
500 501 502
               hmat%data_r(nn, n) = hmat%data_r(nn, n) - a_ex*REAL(vxc(i))
               IF ((n <= 5) .AND. (nn <= 5)) THEN
                  WRITE (1235, '(2i7,3f15.8)') n, nn, hmat%data_r(n, nn), hmat%data_r(nn, n), REAL(vxc(i))
503
               END IF
504
            ELSE
Matthias Redies's avatar
Matthias Redies committed
505
               hmat%data_c(nn, n) = hmat%data_c(nn, n) - a_ex*vxc(i)
506
            ENDIF
507
         END DO
508
      END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
509

510
      CALL timestop("subvxc")
Daniel Wortmann's avatar
Daniel Wortmann committed
511

Matthias Redies's avatar
Matthias Redies committed
512
      DEALLOCATE (bascof)
Daniel Wortmann's avatar
Daniel Wortmann committed
513

514
   END SUBROUTINE subvxc
Daniel Wortmann's avatar
Daniel Wortmann committed
515
END MODULE m_subvxc
516