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

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_subvxc
       8             : 
       9             : CONTAINS
      10             : 
      11           0 :    SUBROUTINE subvxc(lapw, bk, DIMENSION, input, jsp, vr0, atoms, usdus, hybrid, el, ello, sym, &
      12             :                      cell, sphhar, stars, xcpot, mpi, oneD, hmat, vx)
      13             : 
      14             :       USE m_types
      15             :       USE m_judft
      16             :       USE m_intgr, ONLY: intgr3
      17             :       USE m_constants
      18             :       USE m_gaunt, ONLY: gaunt1
      19             :       USE m_wrapper
      20             :       USE m_loddop
      21             :       USE m_radflo
      22             :       USE m_radfun
      23             :       USE m_abcof3
      24             : 
      25             :       IMPLICIT NONE
      26             : 
      27             :       CLASS(t_xcpot), INTENT(IN)    :: xcpot
      28             :       TYPE(t_mpi), INTENT(IN)    :: mpi
      29             :       TYPE(t_dimension), INTENT(IN)    :: dimension
      30             :       TYPE(t_oneD), INTENT(IN)    :: oneD
      31             :       TYPE(t_hybrid), INTENT(INOUT) :: hybrid
      32             :       TYPE(t_input), INTENT(IN)    :: input
      33             :       TYPE(t_sym), INTENT(IN)    :: sym
      34             :       TYPE(t_stars), INTENT(IN)    :: stars
      35             :       TYPE(t_cell), INTENT(IN)    :: cell
      36             :       TYPE(t_sphhar), INTENT(IN)    :: sphhar
      37             :       TYPE(t_atoms), INTENT(IN)    :: atoms
      38             :       TYPE(t_lapw), INTENT(IN)    :: lapw
      39             :       TYPE(t_usdus), INTENT(INOUT) :: usdus
      40             :       TYPE(t_potden), INTENT(IN)    :: vx
      41             :       TYPE(t_mat), INTENT(INOUT) :: hmat
      42             : 
      43             :       ! Scalar Arguments
      44             :       INTEGER, INTENT(IN) :: jsp
      45             : 
      46             :       ! Array Arguments
      47             :       REAL, INTENT(IN) :: vr0(atoms%jmtd, atoms%ntype, input%jspins)               ! just for radial functions
      48             :       REAL, INTENT(IN) :: el(0:atoms%lmaxd, atoms%ntype, input%jspins)
      49             :       REAL, INTENT(IN) :: ello(atoms%nlod, atoms%ntype, input%jspins)
      50             :       REAL, INTENT(IN) :: bk(3)
      51             : 
      52             :       ! Local Scalars
      53             :       INTEGER               ::  ic, indx, m, ig1, ig2, n, nn
      54             :       INTEGER               ::  nlharm, nnbas, typsym, lm
      55             :       INTEGER               ::  noded, nodeu
      56             :       INTEGER               ::  nbasf0
      57             :       INTEGER               ::  i, j, l, ll, l1, l2, m1, m2, j1, j2
      58             :       INTEGER               ::  ok, p1, p2, lh, mh, pp1, pp2
      59             :       INTEGER               ::  igrid, itype, ilharm, istar
      60             :       INTEGER               ::  ineq, iatom, ilo, ilop, ieq, icentry
      61             :       INTEGER               ::  ikvecat, ikvecprevat, invsfct, ikvec, ikvecp
      62             :       INTEGER               ::  lp, mp, pp
      63             :       REAL                  ::  a_ex
      64             :       REAL                  ::  wronk
      65             :       COMPLEX               ::  rc, rr
      66             : 
      67             :       ! Local Arrays
      68             :       INTEGER               ::  gg(3)
      69           0 :       INTEGER               ::  pointer_lo(atoms%nlod, atoms%ntype)
      70             : 
      71           0 :       REAL                  ::  integ(0:sphhar%nlhd, hybrid%maxindx, 0:atoms%lmaxd, hybrid%maxindx, 0:atoms%lmaxd)
      72           0 :       REAL                  ::  grid(atoms%jmtd)
      73           0 :       REAL                  ::  vr(atoms%jmtd, 0:sphhar%nlhd)
      74           0 :       REAL                  ::  f(atoms%jmtd, 2, 0:atoms%lmaxd), g(atoms%jmtd, 2, 0:atoms%lmaxd)
      75           0 :       REAL                  ::  flo(atoms%jmtd, 2, atoms%nlod)
      76           0 :       REAL                  ::  uuilon(atoms%nlod, atoms%ntype), duilon(atoms%nlod, atoms%ntype)
      77           0 :       REAL                  ::  ulouilopn(atoms%nlod, atoms%nlod, atoms%ntype)
      78             : 
      79           0 :       REAL                  ::  bas1(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype)
      80           0 :       REAL                  ::  bas2(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype)
      81             : 
      82           0 :       COMPLEX               ::  vpw(stars%ng3)
      83           0 :       COMPLEX               ::  vxc(hmat%matsize1*(hmat%matsize1 + 1)/2)
      84           0 :       COMPLEX               ::  vrmat(hybrid%maxlmindx, hybrid%maxlmindx)
      85           0 :       COMPLEX               ::  carr(hybrid%maxlmindx, DIMENSION%nvd), carr1(DIMENSION%nvd, DIMENSION%nvd)
      86           0 :       COMPLEX, ALLOCATABLE  ::  ahlp(:, :, :), bhlp(:, :, :)
      87           0 :       COMPLEX, ALLOCATABLE  ::  bascof(:, :, :)
      88             : #ifndef CPP_OLDINTEL
      89           0 :       COMPLEX               ::  bascof_lo(3, -atoms%llod:atoms%llod, 4*atoms%llod + 2, atoms%nlod, atoms%nat)
      90             : #endif
      91             : 
      92           0 :       CALL timestart("subvxc")
      93           0 :       vxc = 0
      94             : 
      95             :       ! Calculate radial functions
      96           0 :       hybrid%nindx = 2
      97           0 :       DO itype = 1, atoms%ntype
      98             : 
      99             :          ! Generate the radial basis-functions for each l
     100           0 :          WRITE (6, '(a,i3,a)') new_LINE('n')//new_LINE('n')//' wavefunction parameters for atom type', itype, ':'
     101           0 :          WRITE (6, '(31x,a,32x,a)') 'radial function', 'energy derivative'
     102             :          WRITE (6, '(a)') '  l    energy            value        '// &
     103           0 :             'derivative    nodes          value        derivative    nodes       norm        wronskian'
     104           0 :          DO l = 0, atoms%lmax(itype)
     105           0 :             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)
     106           0 :             WRITE (6, FMT=8010) l, el(l, itype, jsp), usdus%us(l, itype, jsp), &
     107           0 :                usdus%dus(l, itype, jsp), nodeu, usdus%uds(l, itype, jsp), usdus%duds(l, itype, jsp), noded, &
     108           0 :                usdus%ddn(l, itype, jsp), wronk
     109             :          END DO
     110             : 8010     FORMAT(i3, f10.5, 2(5x, 1p, 2e16.7, i5), 1p, 2e16.7)
     111             : 
     112           0 :          bas1(:, 1, :, itype) = f(:, 1, :)
     113           0 :          bas1(:, 2, :, itype) = g(:, 1, :)
     114           0 :          bas2(:, 1, :, itype) = f(:, 2, :)
     115           0 :          bas2(:, 2, :, itype) = g(:, 2, :)
     116             : 
     117             :          ! Generate the extra radial basis-functions for the local orbitals, if there are any.
     118           0 :          IF (atoms%nlo(itype) >= 1) THEN
     119             :             CALL radflo(atoms, itype, jsp, ello(1, 1, jsp), vr0(1, itype, jsp), f, g, mpi, &
     120           0 :                         usdus, uuilon, duilon, ulouilopn, flo, .TRUE.)
     121             : 
     122           0 :             DO i = 1, atoms%nlo(itype)
     123           0 :                hybrid%nindx(atoms%llo(i, itype), itype) = hybrid%nindx(atoms%llo(i, itype), itype) + 1
     124           0 :                pointer_lo(i, itype) = hybrid%nindx(atoms%llo(i, itype), itype)
     125           0 :                bas1(:, hybrid%nindx(atoms%llo(i, itype), itype), atoms%llo(i, itype), itype) = flo(:, 1, i)
     126           0 :                bas2(:, hybrid%nindx(atoms%llo(i, itype), itype), atoms%llo(i, itype), itype) = flo(:, 2, i)
     127             :             END DO
     128             :          END IF
     129             :       END DO
     130             : 
     131             :       ! Compute APW coefficients
     132             : 
     133             :       ! Calculate bascof
     134           0 :       ALLOCATE (ahlp(DIMENSION%nvd, 0:DIMENSION%lmd, atoms%nat), bhlp(DIMENSION%nvd, 0:DIMENSION%lmd, atoms%nat), stat=ok)
     135           0 :       IF (ok /= 0) STOP 'subvxc: error in allocation of ahlp/bhlp'
     136             : #ifndef CPP_OLDINTEL
     137           0 :       CALL abcof3(input, atoms, sym, jsp, cell, bk, lapw, usdus, oneD, ahlp, bhlp, bascof_lo)
     138             : #endif
     139           0 :       ALLOCATE (bascof(DIMENSION%nvd, 2*(DIMENSION%lmd + 1), atoms%nat), stat=ok)
     140           0 :       IF (ok /= 0) STOP 'subvxc: error in allocation of bascof'
     141             : 
     142           0 :       bascof = 0
     143           0 :       ic = 0
     144           0 :       DO itype = 1, atoms%ntype
     145           0 :          DO ieq = 1, atoms%neq(itype)
     146           0 :             ic = ic + 1
     147           0 :             indx = 0
     148           0 :             DO l = 0, atoms%lmax(itype)
     149           0 :                ll = l*(l + 1)
     150           0 :                DO M = -l, l
     151           0 :                   lm = ll + M
     152           0 :                   DO i = 1, 2
     153           0 :                      indx = indx + 1
     154           0 :                      IF (i == 1) THEN
     155           0 :                         bascof(:, indx, ic) = ahlp(:, lm, ic)
     156             :                      ELSE IF (i == 2) THEN
     157           0 :                         bascof(:, indx, ic) = bhlp(:, lm, ic)
     158             :                      END IF
     159             :                   END DO
     160             :                END DO
     161             :             END DO
     162             :          END DO
     163             :       END DO
     164             : 
     165           0 :       DEALLOCATE (ahlp, bhlp)
     166             : 
     167             :       ! Loop over atom types
     168           0 :       iatom = 0
     169           0 :       DO itype = 1, atoms%ntype
     170             : 
     171           0 :          typsym = atoms%ntypsy(SUM(atoms%neq(:itype - 1)) + 1)
     172           0 :          nlharm = sphhar%nlh(typsym)
     173             : 
     174             :          ! Calculate vxc = vtot - vcoul
     175           0 :          DO l = 0, nlharm
     176           0 :             DO i = 1, atoms%jri(itype)
     177           0 :                IF (l == 0) THEN
     178             :                   ! vr(i,0)= vrtot(i,0,itype)*sfp/rmsh(i,itype) - vrcou(i,0,itype,jsp)
     179           0 :                   vr(i, 0) = vx%mt(i, 0, itype, jsp)! * sfp_const / atoms%rmsh(i,itype)
     180             :                ELSE ! vxc = vtot - vcoul
     181             :                   ! vr(i,l)= vrtot(i,l,itype)-vrcou(i,l,itype,jsp)
     182           0 :                   vr(i, l) = vx%mt(i, l, itype, jsp)
     183             :                END IF
     184             :             END DO
     185             :          END DO
     186             : 
     187             :          ! Calculate MT contribution to vxc matrix elements
     188             :          ! Precompute auxiliary radial integrals
     189           0 :          DO ilharm = 0, nlharm
     190           0 :             i = 0
     191           0 :             DO l1 = 0, atoms%lmax(itype)
     192           0 :                DO p1 = 1, 2
     193           0 :                   i = i + 1
     194           0 :                   j = 0
     195           0 :                   DO l2 = 0, atoms%lmax(itype)
     196           0 :                      DO p2 = 1, 2
     197           0 :                         j = j + 1
     198           0 :                         IF (j <= i) THEN
     199           0 :                            DO igrid = 1, atoms%jri(itype)
     200             :                               grid(igrid) = vr(igrid, ilharm)*(bas1(igrid, p1, l1, itype)*bas1(igrid, p2, l2, itype) + &
     201           0 :                                                                bas2(igrid, p1, l1, itype)*bas2(igrid, p2, l2, itype))
     202             :                            END DO
     203           0 :                            CALL intgr3(grid, atoms%rmsh(:, itype), atoms%dx(itype), atoms%jri(itype), integ(ilharm, p1, l1, p2, l2))
     204           0 :                            integ(ilharm, p2, l2, p1, l1) = integ(ilharm, p1, l1, p2, l2)
     205             :                         END IF
     206             :                      END DO
     207             :                   END DO
     208             :                END DO
     209             :             END DO
     210             :          END DO
     211             : 
     212             :          ! Calculate muffin tin contribution to vxc matrix
     213           0 :          vrmat = 0
     214             : 
     215           0 :          j1 = 0
     216           0 :          DO l1 = 0, atoms%lmax(itype) ! loop: left basis function
     217           0 :             DO m1 = -l1, l1
     218           0 :                DO p1 = 1, 2
     219           0 :                   j1 = j1 + 1
     220           0 :                   j2 = 0
     221           0 :                   DO l2 = 0, atoms%lmax(itype) ! loop: right basis function
     222           0 :                      DO m2 = -l2, l2
     223           0 :                         DO p2 = 1, 2
     224           0 :                            j2 = j2 + 1
     225           0 :                            rr = 0
     226           0 :                            DO ilharm = 0, nlharm ! loop: lattice harmonics of vxc
     227           0 :                               l = sphhar%llh(ilharm, typsym)
     228           0 :                               DO i = 1, sphhar%nmem(ilharm, typsym)
     229           0 :                                  M = sphhar%mlh(i, ilharm, typsym)
     230           0 :                                  rc = sphhar%clnu(i, ilharm, typsym)*gaunt1(l1, l, l2, m1, M, m2, atoms%lmaxd)
     231           0 :                                  rr = rr + integ(ilharm, p1, l1, p2, l2)*rc
     232             :                               END DO
     233             :                            END DO
     234           0 :                            rc = CMPLX(0, 1)**(l2 - l1) ! adjusts to a/b/ccof-scaling
     235           0 :                            vrmat(j1, j2) = rr*rc
     236             :                         END DO
     237             :                      END DO
     238             :                   END DO
     239             :                END DO
     240             :             END DO
     241             :          END DO
     242           0 :          nnbas = j1
     243             : 
     244             :          ! Project on bascof
     245           0 :          DO ineq = 1, atoms%neq(itype)
     246           0 :             iatom = iatom + 1
     247             :             carr(:nnbas, :lapw%nv(jsp)) = CONJG(MATMUL(vrmat(:nnbas, :nnbas), &
     248           0 :                                                        TRANSPOSE(bascof(:lapw%nv(jsp), :nnbas, iatom))))
     249           0 :             carr1(:lapw%nv(jsp), :lapw%nv(jsp)) = MATMUL(bascof(:lapw%nv(jsp), :nnbas, iatom), carr(:nnbas, :lapw%nv(jsp)))
     250           0 :             ic = 0
     251           0 :             DO j = 1, lapw%nv(jsp)
     252             :                ! carr(:nnbas) =  matmul(vrmat(:nnbas,:nnbas),bascof(j,:nnbas,iatom) )
     253           0 :                DO i = 1, j
     254           0 :                   ic = ic + 1
     255           0 :                   vxc(ic) = vxc(ic) + carr1(i, j)
     256             :                   ! vxc(ic) = vxc(ic) + conjg(dotprod ( bascof(i,:nnbas,iatom),carr(:nnbas) ))
     257             :                END DO
     258             :             END DO
     259             :          END DO
     260             :       END DO ! End loop over atom types
     261             : 
     262             :       ! Calculate plane wave contribution
     263           0 :       DO i = 1, stars%ng3
     264           0 :          vpw(i) = vx%pw_w(i, jsp)
     265             :          ! vpw(i)=vpwtot(i)-vpwcou(i,jsp)
     266             :       END DO
     267             : 
     268             :       ! Calculate vxc-matrix,  left basis function (ig1)
     269             :       !                        right basis function (ig2)
     270           0 :       ic = 0
     271           0 :       DO ig1 = 1, lapw%nv(jsp)
     272           0 :          DO ig2 = 1, ig1
     273           0 :             ic = ic + 1
     274           0 :             gg(1) = lapw%k1(ig1, jsp) - lapw%k1(ig2, jsp)
     275           0 :             gg(2) = lapw%k2(ig1, jsp) - lapw%k2(ig2, jsp)
     276           0 :             gg(3) = lapw%k3(ig1, jsp) - lapw%k3(ig2, jsp)
     277           0 :             istar = stars%ig(gg(1), gg(2), gg(3))
     278           0 :             IF (istar /= 0) THEN
     279           0 :                vxc(ic) = vxc(ic) + stars%rgphs(gg(1), gg(2), gg(3))*vpw(istar)
     280             :             ELSE
     281           0 :                IF (mpi%irank == 0) THEN
     282           0 :                   WRITE (6, '(A,/6I5)') 'Warning: Gi-Gj not in any star:', &
     283           0 :                      lapw%k1(ig1, jsp), lapw%k2(ig1, jsp), lapw%k3(ig1, jsp), &
     284           0 :                      lapw%k1(ig2, jsp), lapw%k2(ig2, jsp), lapw%k3(ig2, jsp)
     285             :                END IF
     286             :             END IF
     287             :          END DO
     288             :       END DO
     289             : 
     290             :       ! Calculate local orbital contribution
     291           0 :       IF (ANY(atoms%nlo /= 0)) THEN
     292             : 
     293           0 :          nbasf0 = lapw%nv(jsp)*(lapw%nv(jsp) + 1)/2  ! number of pure APW contributions
     294           0 :          icentry = nbasf0                          ! icentry counts the entry in the matrix vxc
     295           0 :          iatom = 0
     296           0 :          ikvecat = 0
     297           0 :          ikvecprevat = 0
     298             : 
     299           0 :          DO itype = 1, atoms%ntype
     300             : 
     301           0 :             typsym = atoms%ntypsy(SUM(atoms%neq(:itype - 1)) + 1)
     302           0 :             nlharm = sphhar%nlh(typsym)
     303             : 
     304             :             ! Calculate vxc = vtot - vcoul
     305           0 :             DO l = 0, nlharm
     306           0 :                DO i = 1, atoms%jri(itype)
     307           0 :                   IF (l == 0) THEN
     308             :                      ! vr(i,0)= vrtot(i,0,itype)*sfp/rmsh(i,itype) -  vrcou(i,0,itype,jsp)
     309           0 :                      vr(i, 0) = vx%mt(i, 0, itype, jsp)*sfp_const/atoms%rmsh(i, itype)
     310             :                   ELSE ! vxc = vtot - vcoul
     311           0 :                      vr(i, l) = vx%mt(i, l, itype, jsp)                    !
     312             :                      ! vr(i,l)=  vrtot(i,l,itype)-vrcou(i,l,itype,jsp)
     313             :                   END IF
     314             :                END DO
     315             :             END DO
     316             : 
     317             :             ! Precompute auxiliary radial integrals
     318           0 :             DO ilharm = 0, nlharm
     319           0 :                i = 0
     320           0 :                DO l1 = 0, atoms%lmax(itype)
     321           0 :                   DO p1 = 1, hybrid%nindx(l1, itype)
     322           0 :                      i = i + 1
     323           0 :                      j = 0
     324           0 :                      DO l2 = 0, atoms%lmax(itype)
     325           0 :                         DO p2 = 1, hybrid%nindx(l2, itype)
     326           0 :                            j = j + 1
     327           0 :                            IF (j <= i) THEN
     328           0 :                               DO igrid = 1, atoms%jri(itype)
     329             :                                  grid(igrid) = vr(igrid, ilharm)*(bas1(igrid, p1, l1, itype)*bas1(igrid, p2, l2, itype) + &
     330           0 :                                                                   bas2(igrid, p1, l1, itype)*bas2(igrid, p2, l2, itype))
     331             :                               END DO
     332           0 :                               CALL intgr3(grid, atoms%rmsh(:, itype), atoms%dx(itype), atoms%jri(itype), integ(ilharm, p1, l1, p2, l2))
     333           0 :                               integ(ilharm, p2, l2, p1, l1) = integ(ilharm, p1, l1, p2, l2)
     334             :                            END IF
     335             :                         END DO
     336             :                      END DO
     337             :                   END DO
     338             :                END DO
     339             :             END DO
     340             : 
     341           0 :             DO ieq = 1, atoms%neq(itype)
     342           0 :                iatom = iatom + 1
     343           0 :                IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN
     344             : 
     345           0 :                   IF (atoms%invsat(iatom) == 0) invsfct = 1
     346           0 :                   IF (atoms%invsat(iatom) == 1) invsfct = 2
     347             : 
     348           0 :                   DO ilo = 1, atoms%nlo(itype)
     349             : #ifdef CPP_OLDINTEL
     350             :                      CALL judft_error("no LOs & hybrid with old intel compiler!", calledby="subvxc.F90")
     351             : #else
     352           0 :                      l1 = atoms%llo(ilo, itype)
     353           0 :                      DO ikvec = 1, invsfct*(2*l1 + 1)
     354           0 :                         DO m1 = -l1, l1
     355           0 :                            DO p1 = 1, 3
     356           0 :                               IF (p1 == 3) THEN
     357           0 :                                  pp1 = pointer_lo(ilo, itype)
     358             :                               ELSE
     359             :                                  pp1 = p1
     360             :                               END IF
     361             : 
     362           0 :                               IF (hybrid%nindx(l1, itype) <= 2) STOP 'subvxc: error hybrid%nindx'
     363             : 
     364           0 :                               lm = 0
     365             : 
     366             :                               !loop over APW
     367           0 :                               DO l2 = 0, atoms%lmax(itype)
     368           0 :                                  DO m2 = -l2, l2
     369           0 :                                     DO p2 = 1, 2
     370           0 :                                        lm = lm + 1
     371             : 
     372           0 :                                        rr = 0
     373           0 :                                        DO ilharm = 0, nlharm
     374           0 :                                           lh = sphhar%llh(ilharm, typsym)
     375           0 :                                           DO i = 1, sphhar%nmem(ilharm, typsym)
     376           0 :                                              mh = sphhar%mlh(i, ilharm, typsym)
     377           0 :                                              rc = sphhar%clnu(i, ilharm, typsym)*gaunt1(l1, lh, l2, m1, mh, m2, atoms%lmaxd)
     378           0 :                                              rr = rr + integ(ilharm, p2, l2, pp1, l1)*rc
     379             :                                           END DO
     380             :                                        END DO
     381             : 
     382           0 :                                        rc = CMPLX(0.0, 1.0)**(l2 - l1) ! adjusts to a/b/ccof-scaling
     383             : 
     384             :                                        ! ic counts the entry in vxc
     385           0 :                                        ic = icentry
     386           0 :                                        DO i = 1, lapw%nv(jsp)
     387           0 :                                           ic = ic + 1
     388           0 :                                           IF (hmat%l_real) THEN
     389             :                                              vxc(ic) = vxc(ic) + invsfct*REAL(rr*rc*bascof(i, lm, iatom)* &
     390           0 :                                                                               CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom)))
     391             :                                           ELSE
     392             :                                              vxc(ic) = vxc(ic) + rr*rc*bascof(i, lm, iatom)* &
     393           0 :                                                        CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom))
     394             :                                           END IF
     395             :                                        END DO
     396             :                                     END DO  !p2
     397             :                                  END DO  ! m2
     398             :                               END DO ! l2 ->  loop over APW
     399             : 
     400             :                               ! calcualte matrix-elements with local orbitals at the same atom
     401           0 :                               IF (ic /= icentry + lapw%nv(jsp)) STOP 'subvxc: error counting ic'
     402             : 
     403           0 :                               ic = ic + ikvecprevat
     404             : 
     405           0 :                               DO ilop = 1, ilo - 1
     406           0 :                                  lp = atoms%llo(ilop, itype)
     407           0 :                                  DO ikvecp = 1, invsfct*(2*lp + 1)
     408           0 :                                     ic = ic + 1
     409           0 :                                     DO mp = -lp, lp
     410           0 :                                        DO pp = 1, 3
     411           0 :                                           IF (pp == 3) THEN
     412           0 :                                              pp2 = pointer_lo(ilop, itype)
     413             :                                           ELSE
     414             :                                              pp2 = pp
     415             :                                           END IF
     416             : 
     417           0 :                                           rr = 0
     418           0 :                                           DO ilharm = 0, nlharm
     419           0 :                                              lh = sphhar%llh(ilharm, typsym)
     420           0 :                                              DO i = 1, sphhar%nmem(ilharm, typsym)
     421           0 :                                                 mh = sphhar%mlh(i, ilharm, typsym)
     422           0 :                                                 rc = sphhar%clnu(i, ilharm, typsym)*gaunt1(l1, lh, lp, m1, mh, mp, atoms%lmaxd)
     423           0 :                                                 rr = rr + integ(ilharm, pp2, lp, pp1, l1)*rc
     424             :                                              END DO
     425             :                                           END DO
     426             : 
     427           0 :                                           rc = CMPLX(0.0, 1.0)**(lp - l1) ! adjusts to a/b/ccof-scaling
     428             : 
     429           0 :                                           IF (hmat%l_real) THEN
     430             :                                              vxc(ic) = vxc(ic) + invsfct*REAL(rr*rc*bascof_lo(pp, mp, ikvecp, ilop, iatom)* &
     431           0 :                                                                               CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom)))
     432             :                                           ELSE
     433             :                                              vxc(ic) = vxc(ic) + rr*rc*bascof_lo(pp, mp, ikvecp, ilop, iatom)* &
     434           0 :                                                        CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom))
     435             :                                           END IF
     436             :                                        END DO ! pp
     437             :                                     END DO ! mp
     438             : 
     439             :                                  END DO !ikvecp
     440             :                               END DO ! ilop
     441             : 
     442             :                               ! calculate matrix-elements of one local orbital with itself
     443           0 :                               DO ikvecp = 1, ikvec
     444           0 :                                  ic = ic + 1
     445             : 
     446           0 :                                  lp = l1
     447           0 :                                  ilop = ilo
     448           0 :                                  DO mp = -lp, lp
     449           0 :                                     DO pp = 1, 3
     450           0 :                                        IF (pp == 3) THEN
     451           0 :                                           pp2 = pointer_lo(ilop, itype)
     452             :                                        ELSE
     453             :                                           pp2 = pp
     454             :                                        END IF
     455             : 
     456           0 :                                        rr = 0
     457           0 :                                        DO ilharm = 0, nlharm
     458           0 :                                           lh = sphhar%llh(ilharm, typsym)
     459           0 :                                           DO i = 1, sphhar%nmem(ilharm, typsym)
     460           0 :                                              mh = sphhar%mlh(i, ilharm, typsym)
     461           0 :                                              rc = sphhar%clnu(i, ilharm, typsym)*gaunt1(l1, lh, lp, m1, mh, mp, atoms%lmaxd)
     462           0 :                                              rr = rr + integ(ilharm, pp2, lp, pp1, l1)*rc
     463             :                                           END DO
     464             :                                        END DO
     465             : 
     466           0 :                                        rc = CMPLX(0.0, 1.0)**(lp - l1) ! adjusts to a/b/ccof-scaling
     467             : 
     468           0 :                                        IF (hmat%l_real) THEN
     469             :                                           vxc(ic) = vxc(ic) + invsfct*REAL(rr*rc*bascof_lo(pp, mp, ikvecp, ilop, iatom)* &
     470           0 :                                                                            CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom)))
     471             :                                        ELSE
     472             :                                           vxc(ic) = vxc(ic) + rr*rc*bascof_lo(pp, mp, ikvecp, ilop, iatom)* &
     473           0 :                                                     CONJG(bascof_lo(p1, m1, ikvec, ilo, iatom))
     474             :                                        END IF
     475             :                                     END DO ! pp
     476             :                                  END DO ! mp
     477             :                               END DO ! ikvecp
     478             :                            END DO  ! p1
     479             :                         END DO  ! m1
     480           0 :                         icentry = ic
     481             :                      END DO !ikvec
     482           0 :                      ikvecat = ikvecat + invsfct*(2*l1 + 1)
     483             : #endif
     484             :                   END DO  ! ilo
     485           0 :                   ikvecprevat = ikvecprevat + ikvecat
     486           0 :                   ikvecat = 0
     487             :                END IF ! atoms%invsat(iatom)
     488             :             END DO ! ieq
     489             :          END DO !itype
     490             :       END IF ! if any atoms%llo
     491             : 
     492             :       !initialize weighting factor
     493           0 :       a_ex = xcpot%get_exchange_weight()
     494             : 
     495           0 :       i = 0
     496           0 :       DO n = 1, hmat%matsize1
     497           0 :          DO nn = 1, n
     498           0 :             i = i + 1
     499           0 :             IF (hmat%l_real) THEN
     500           0 :                hmat%data_r(nn, n) = hmat%data_r(nn, n) - a_ex*REAL(vxc(i))
     501           0 :                IF ((n <= 5) .AND. (nn <= 5)) THEN
     502           0 :                   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
     505           0 :                hmat%data_c(nn, n) = hmat%data_c(nn, n) - a_ex*vxc(i)
     506             :             ENDIF
     507             :          END DO
     508             :       END DO
     509             : 
     510           0 :       CALL timestop("subvxc")
     511             : 
     512           0 :       DEALLOCATE (bascof)
     513             : 
     514           0 :    END SUBROUTINE subvxc
     515             : END MODULE m_subvxc
     516             : 

Generated by: LCOV version 1.13