IffGit has a new shared runner for building Docker images in GitLab CI. Visit https://iffgit.fz-juelich.de/examples/ci-docker-in-docker for more details.

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