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

          Line data    Source code
       1             : !     Calculates the HF exchange term
       2             : !
       3             : !                                          s          s*          s            s*
       4             : !                                       phi    (r) phi     (r) phi     (r') phi    (r')
       5             : !                         occ.             n_1k       n'k+q       n'k+q        n_2k
       6             : !     exchange(n,q)  =  - SUM  INT INT  ------------------------------------------- dr dr'
       7             : !                         k,n'                           | r - r' |
       8             : !
       9             : !                         occ                  s          s    ~        ~       s         s
      10             : !                    =  - SUM  SUM  v     < phi      | phi     M    > < M    phi     | phi      >
      11             : !                         k,n' I,J   k,IJ      n'k+q      n_1k  q,I      q,J    n_2k      n'k+q
      12             : !
      13             : !     for the different combinations of n_1 and n_2 and  n' runs over the core states.
      14             : !     ( n_1,n_2:  valence-valence, core-core,core-valence )
      15             : !
      16             : !     It is done directly without employing the mixed basis set.
      17             : 
      18             : MODULE m_exchange_core
      19             : 
      20             : CONTAINS
      21             : 
      22           0 :    SUBROUTINE exchange_vccv(nk, atoms, hybrid, hybdat, DIMENSION, jsp, lapw, &
      23           0 :                             maxbands, mnobd, mpi, degenerat, symequivalent, results, &
      24           0 :                             ex_vv_r, ex_vv_c, l_real)
      25             : 
      26             :       USE m_constants
      27             :       USE m_util
      28             :       USE m_wrapper
      29             :       USE m_types
      30             :       USE m_io_hybrid
      31             :       IMPLICIT NONE
      32             : 
      33             :       TYPE(t_hybdat), INTENT(IN)   :: hybdat
      34             :       TYPE(t_results), INTENT(INOUT)   :: results
      35             :       TYPE(t_mpi), INTENT(IN)   :: mpi
      36             :       TYPE(t_dimension), INTENT(IN)   :: DIMENSION
      37             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
      38             :       TYPE(t_atoms), INTENT(IN)   :: atoms
      39             :       TYPE(t_lapw), INTENT(IN)   :: lapw
      40             : 
      41             :       !     -scalars -
      42             :       INTEGER, INTENT(IN)      :: jsp
      43             :       INTEGER, INTENT(IN)      ::nk, maxbands, mnobd
      44             :       !     - arays -
      45             :       INTEGER, INTENT(IN)      ::  degenerat(hybrid%ne_eig(nk))
      46             :       LOGICAL, INTENT(IN)      :: l_real
      47             :       REAL, INTENT(INOUT)  ::  ex_vv_r(:, :, :)!(maxbands,mnobd,nkpti)
      48             :       COMPLEX, INTENT(INOUT)  ::  ex_vv_c(:, :, :)!(maxbands,mnobd,nkpti)
      49             :       LOGICAL                 ::  symequivalent(COUNT(degenerat >= 1), COUNT(degenerat >= 1))
      50             : 
      51             :       !     - local scalars -
      52             :       INTEGER                 ::  iatom, ieq, itype, ic, l, l1, l2, ll, lm, m1, m2, p1, p2, n, n1, n2, i, j
      53             :       INTEGER                 ::  iband1, iband2, ndb1, ndb2, ic1, ic2
      54             :       INTEGER                 ::  m
      55             : 
      56             :       REAL                    ::  time1, time2
      57             :       REAL                    ::  rdum
      58             :       REAL                    ::  sum_offdia
      59             : 
      60             :       COMPLEX                 ::  cdum
      61             : 
      62             :       !     - local arrays -
      63           0 :       INTEGER, ALLOCATABLE     ::  larr(:), larr2(:)
      64           0 :       INTEGER, ALLOCATABLE     ::  parr(:), parr2(:)
      65             : 
      66           0 :       REAL                    ::  integrand(atoms%jmtd)
      67           0 :       REAL                    ::  primf1(atoms%jmtd), primf2(atoms%jmtd)
      68           0 :       REAL, ALLOCATABLE        ::  fprod(:, :), fprod2(:, :)
      69           0 :       REAL, ALLOCATABLE        ::  integral(:, :)
      70             : 
      71           0 :       COMPLEX                 ::  cmt(DIMENSION%neigd, hybrid%maxlmindx, atoms%nat)
      72           0 :       COMPLEX                 ::  exchange(hybrid%nbands(nk), hybrid%nbands(nk))
      73           0 :       COMPLEX, ALLOCATABLE     ::  carr(:, :), carr2(:, :), carr3(:, :)
      74             : 
      75           0 :       LOGICAL                 ::  ldum(hybrid%nbands(nk), hybrid%nbands(nk))
      76             : 
      77           0 :       WRITE (6, '(A)') new_LINE('n')//new_LINE('n')//'### valence-core-core-valence exchange ###'
      78           0 :       WRITE (6, '(A)') new_LINE('n')//'        k-point       band    exchange (core contribution)'
      79             : 
      80             :       ! read in mt wavefunction coefficients from file cmt
      81           0 :       CALL read_cmt(cmt, nk)
      82           0 :       ALLOCATE (fprod(atoms%jmtd, 5), larr(5), parr(5))
      83             : 
      84             :       ! generate ldum(nbands(nk),nbands(nk)), which is true if the corresponding matrix entry is non-zero
      85             :       ic1 = 0
      86           0 :       ldum = .FALSE.
      87           0 :       DO iband1 = 1, hybrid%nbands(nk)
      88           0 :          ndb1 = degenerat(iband1)
      89           0 :          IF (ndb1 >= 1) THEN
      90           0 :             ic1 = ic1 + 1
      91           0 :             ic2 = 0
      92           0 :             DO iband2 = 1, hybrid%nbands(nk)
      93           0 :                ndb2 = degenerat(iband2)
      94           0 :                IF (ndb2 >= 1) THEN
      95           0 :                   ic2 = ic2 + 1
      96           0 :                   IF (symequivalent(ic2, ic1)) THEN
      97           0 :                      IF (ndb1 /= ndb2) STOP 'exchange: failure symequivalent'
      98           0 :                      DO i = 0, ndb1 - 1
      99           0 :                         DO j = 0, ndb2 - 1
     100           0 :                            ldum(iband1 + i, iband2 + j) = .TRUE.
     101             :                         END DO
     102             :                      END DO
     103             : 
     104             :                   END IF
     105             :                END IF
     106             :             END DO
     107             :          END IF
     108             :       END DO
     109             : 
     110           0 :       exchange = 0
     111           0 :       iatom = 0
     112           0 :       rdum = 0
     113           0 :       DO itype = 1, atoms%ntype
     114           0 :          DO ieq = 1, atoms%neq(itype)
     115           0 :             iatom = iatom + 1
     116           0 :             DO l1 = 0, hybdat%lmaxc(itype)
     117           0 :                DO p1 = 1, hybdat%nindxc(l1, itype)
     118             : 
     119           0 :                   DO l = 0, hybrid%lcutm1(itype)
     120             : 
     121             :                      ! Define core-valence product functions
     122             : 
     123           0 :                      n = 0
     124           0 :                      DO l2 = 0, atoms%lmax(itype)
     125           0 :                         IF (l < ABS(l1 - l2) .OR. l > l1 + l2) CYCLE
     126             : 
     127           0 :                         DO p2 = 1, hybrid%nindx(l2, itype)
     128           0 :                            n = n + 1
     129           0 :                            M = SIZE(fprod, 2)
     130           0 :                            IF (n > M) THEN
     131           0 :                               ALLOCATE (fprod2(atoms%jmtd, M), larr2(M), parr2(M))
     132           0 :                               fprod2 = fprod; larr2 = larr; parr2 = parr
     133           0 :                               DEALLOCATE (fprod, larr, parr)
     134           0 :                               ALLOCATE (fprod(atoms%jmtd, M + 5), larr(M + 5), parr(M + 5))
     135           0 :                               fprod(:, :M) = fprod2
     136           0 :                               larr(:M) = larr2
     137           0 :                               parr(:M) = parr2
     138           0 :                               DEALLOCATE (fprod2, larr2, parr2)
     139             :                            END IF
     140             :                            fprod(:atoms%jri(itype), n) = (hybdat%core1(:atoms%jri(itype), p1, l1, itype)*hybdat%bas1(:atoms%jri(itype), p2, l2, itype) &
     141           0 :                                                           + hybdat%core2(:atoms%jri(itype), p1, l1, itype)*hybdat%bas2(:atoms%jri(itype), p2, l2, itype))/atoms%rmsh(:atoms%jri(itype), itype)
     142           0 :                            larr(n) = l2
     143           0 :                            parr(n) = p2
     144             :                         END DO
     145             :                      END DO
     146             : 
     147             :                      ! Evaluate radial integrals (special part of Coulomb matrix : contribution from single MT)
     148             : 
     149           0 :                      ALLOCATE (integral(n, n), carr(n, hybrid%nbands(nk)), carr2(n, lapw%nv(jsp)), carr3(n, lapw%nv(jsp)))
     150             : 
     151           0 :                      DO i = 1, n
     152           0 :                         CALL primitivef(primf1, fprod(:, i)*atoms%rmsh(:, itype)**(l + 1), atoms%rmsh, atoms%dx, atoms%jri, atoms%jmtd, itype, atoms%ntype)
     153           0 :                         CALL primitivef(primf2, fprod(:atoms%jri(itype), i)/atoms%rmsh(:atoms%jri(itype), itype)**l, atoms%rmsh, atoms%dx, atoms%jri, atoms%jmtd, -itype, atoms%ntype)  ! -itype is to enforce inward integration
     154             : 
     155           0 :                         primf1(:atoms%jri(itype)) = primf1(:atoms%jri(itype))/atoms%rmsh(:atoms%jri(itype), itype)**l
     156           0 :                         primf2(:atoms%jri(itype)) = primf2(:atoms%jri(itype))*atoms%rmsh(:atoms%jri(itype), itype)**(l + 1)
     157           0 :                         DO j = 1, n
     158           0 :                            integrand = fprod(:, j)*(primf1 + primf2)
     159             :                            integral(i, j) = fpi_const/(2*l + 1)*intgrf(integrand, atoms%jri, atoms%jmtd, atoms%rmsh, &
     160           0 :                                                                        atoms%dx, atoms%ntype, itype, hybdat%gridf)
     161             :                         END DO
     162             :                      END DO
     163             : 
     164             :                      ! Add everything up
     165             : 
     166           0 :                      DO m1 = -l1, l1
     167           0 :                         DO M = -l, l
     168           0 :                            m2 = m1 + M
     169             : 
     170           0 :                            carr = 0
     171           0 :                            ic = 0
     172           0 :                            DO n1 = 1, hybrid%nbands(nk)
     173             : 
     174           0 :                               DO i = 1, n
     175           0 :                                  ll = larr(i)
     176           0 :                                  IF (ABS(m2) > ll) CYCLE
     177             : 
     178           0 :                                  lm = SUM((/((2*l2 + 1)*hybrid%nindx(l2, itype), l2=0, ll - 1)/)) + (m2 + ll)*hybrid%nindx(ll, itype) + parr(i)
     179             : 
     180           0 :                                  carr(i, n1) = cmt(n1, lm, iatom)*gaunt(l1, ll, l, m1, m2, M, hybdat%maxfac, hybdat%fac, hybdat%sfac)
     181             : 
     182             :                               END DO
     183           0 :                               DO n2 = 1, n1
     184           0 :                                  IF (ldum(n2, n1)) THEN
     185           0 :                                     ic = ic + 1
     186           0 :                                     exchange(n2, n1) = exchange(n2, n1) + dot_PRODUCT(carr(:, n1), MATMUL(integral, carr(:, n2)))
     187             :                                  END IF
     188             :                               END DO
     189             :                            END DO
     190             :                         END DO
     191             :                      END DO
     192             : 
     193           0 :                      DEALLOCATE (integral, carr, carr2, carr3)
     194             : 
     195             :                   END DO
     196             :                END DO
     197             :             END DO
     198             :          END DO
     199             :       END DO
     200             : 
     201           0 :       IF (l_real) THEN
     202           0 :          IF (ANY(ABS(AIMAG(exchange)) > 1e-10)) THEN
     203           0 :             IF (mpi%irank == 0) WRITE (6, '(A)') 'exchangeCore: Warning! Unusually large imaginary component.'
     204           0 :             WRITE (*, *) MAXVAL(ABS(AIMAG(exchange)))
     205           0 :             STOP 'exchangeCore: Unusually large imaginary component.'
     206             :          END IF
     207             :       ENDIF
     208             : 
     209             :       ! add the core-valence contribution to the exchange matrix ex_vv
     210             : 
     211             :       !      ic         = 0
     212           0 :       sum_offdia = 0
     213           0 :       IF (l_real) THEN
     214           0 :          DO n1 = 1, hybrid%nobd(nk)
     215           0 :             DO n2 = 1, hybrid%nbands(nk)
     216           0 :                ex_vv_r(n2, n1, nk) = ex_vv_r(n2, n1, nk) - exchange(n1, n2)
     217           0 :                IF (n1 /= n2) sum_offdia = sum_offdia + 2*ABS(exchange(n1, n2))
     218             :             END DO
     219             :          END DO
     220             :       ELSE
     221           0 :          DO n1 = 1, hybrid%nobd(nk)
     222           0 :             DO n2 = 1, hybrid%nbands(nk)
     223           0 :                ex_vv_c(n2, n1, nk) = ex_vv_c(n2, n1, nk) - exchange(n1, n2)
     224           0 :                IF (n1 /= n2) sum_offdia = sum_offdia + 2*ABS(exchange(n1, n2))
     225             :             END DO
     226             :          END DO
     227             :       END IF
     228             : 
     229           0 :       DO n1 = 1, hybrid%nobd(nk)
     230           0 :          results%te_hfex%core = results%te_hfex%core - results%w_iks(n1, nk, jsp)*exchange(n1, n1)
     231             :       END DO
     232             : 
     233           0 :       WRITE (6, '(A,F20.15)') 'sum of the absolut real part of the non diagonal elements', sum_offdia
     234             : 
     235           0 :    END SUBROUTINE exchange_vccv
     236             : 
     237           0 :    SUBROUTINE exchange_vccv1(nk, atoms, hybrid, hybdat, DIMENSION, jsp, lapw, &
     238           0 :                              nsymop, nsest, indx_sest, mpi, a_ex, results, mat_ex)
     239             : 
     240             :       USE m_constants
     241             :       USE m_util
     242             :       USE m_wrapper
     243             :       USE m_types
     244             :       USE m_io_hybrid
     245             :       IMPLICIT NONE
     246             : 
     247             :       TYPE(t_hybdat), INTENT(IN)   :: hybdat
     248             :       TYPE(t_results), INTENT(INOUT)   :: results
     249             :       TYPE(t_mpi), INTENT(IN)   :: mpi
     250             :       TYPE(t_dimension), INTENT(IN)   :: DIMENSION
     251             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
     252             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     253             :       TYPE(t_lapw), INTENT(IN)   :: lapw
     254             : 
     255             :       !     -scalars -
     256             :       INTEGER, INTENT(IN)      :: jsp
     257             :       INTEGER, INTENT(IN)      :: nk
     258             :       INTEGER, INTENT(IN)      ::  nsymop
     259             :       REAL, INTENT(IN)         ::  a_ex
     260             :       !     - arays -
     261             :       INTEGER, INTENT(IN)      ::  nsest(hybrid%nbands(nk)), indx_sest(hybrid%nbands(nk), hybrid%nbands(nk))
     262             : 
     263             :       TYPE(t_mat), INTENT(INOUT):: mat_ex
     264             :       !     - local scalars -
     265             :       INTEGER                 ::  iatom, ieq, itype, ic, l, l1, l2
     266             :       INTEGER                 ::  ll, lm, m1, m2, p1, p2, n, n1, n2, nn2, i, j
     267             :       INTEGER                 ::  iband1, iband2, ndb1, ndb2, ic1, ic2
     268             :       INTEGER                 ::  m
     269             : 
     270             :       REAL                    ::  time1, time2
     271             :       REAL                    ::  rdum
     272             :       REAL                    ::  sum_offdia
     273             : 
     274             :       COMPLEX                 ::  cdum
     275             :       !     - local arrays -
     276           0 :       INTEGER, ALLOCATABLE     ::  larr(:), larr2(:)
     277           0 :       INTEGER, ALLOCATABLE     ::  parr(:), parr2(:)
     278             : 
     279           0 :       REAL                    ::  integrand(atoms%jmtd)
     280           0 :       REAL                    ::  primf1(atoms%jmtd), primf2(atoms%jmtd)
     281           0 :       REAL, ALLOCATABLE        ::  fprod(:, :), fprod2(:, :)
     282           0 :       REAL, ALLOCATABLE        ::  integral(:, :)
     283             : 
     284           0 :       COMPLEX                 ::  cmt(DIMENSION%neigd, hybrid%maxlmindx, atoms%nat)
     285           0 :       COMPLEX                 ::  exchange(hybrid%nbands(nk), hybrid%nbands(nk))
     286           0 :       COMPLEX, ALLOCATABLE     ::  carr(:, :), carr2(:, :), carr3(:, :)
     287             : 
     288             :       LOGICAL                 ::  ldum(hybrid%nbands(nk), hybrid%nbands(nk))
     289             : 
     290             :       ! read in mt wavefunction coefficients from file cmt
     291             : 
     292           0 :       CALL read_cmt(cmt, nk)
     293             : 
     294           0 :       ALLOCATE (fprod(atoms%jmtd, 5), larr(5), parr(5))
     295             : 
     296           0 :       exchange = 0
     297           0 :       iatom = 0
     298           0 :       rdum = 0
     299           0 :       DO itype = 1, atoms%ntype
     300           0 :          DO ieq = 1, atoms%neq(itype)
     301           0 :             iatom = iatom + 1
     302           0 :             DO l1 = 0, hybdat%lmaxc(itype)
     303           0 :                DO p1 = 1, hybdat%nindxc(l1, itype)
     304             : 
     305           0 :                   DO l = 0, hybrid%lcutm1(itype)
     306             : 
     307             :                      ! Define core-valence product functions
     308             : 
     309           0 :                      n = 0
     310           0 :                      DO l2 = 0, atoms%lmax(itype)
     311           0 :                         IF (l < ABS(l1 - l2) .OR. l > l1 + l2) CYCLE
     312             : 
     313           0 :                         DO p2 = 1, hybrid%nindx(l2, itype)
     314           0 :                            n = n + 1
     315           0 :                            M = SIZE(fprod, 2)
     316           0 :                            IF (n > M) THEN
     317           0 :                               ALLOCATE (fprod2(atoms%jmtd, M), larr2(M), parr2(M))
     318           0 :                               fprod2 = fprod; larr2 = larr; parr2 = parr
     319           0 :                               DEALLOCATE (fprod, larr, parr)
     320           0 :                               ALLOCATE (fprod(atoms%jmtd, M + 5), larr(M + 5), parr(M + 5))
     321           0 :                               fprod(:, :M) = fprod2
     322           0 :                               larr(:M) = larr2
     323           0 :                               parr(:M) = parr2
     324           0 :                               DEALLOCATE (fprod2, larr2, parr2)
     325             :                            END IF
     326             :                            fprod(:atoms%jri(itype), n) = (hybdat%core1(:atoms%jri(itype), p1, l1, itype)*hybdat%bas1(:atoms%jri(itype), p2, l2, itype) &
     327           0 :                                                           + hybdat%core2(:atoms%jri(itype), p1, l1, itype)*hybdat%bas2(:atoms%jri(itype), p2, l2, itype))/atoms%rmsh(:atoms%jri(itype), itype)
     328           0 :                            larr(n) = l2
     329           0 :                            parr(n) = p2
     330             :                         END DO
     331             :                      END DO
     332             : 
     333             :                      ! Evaluate radial integrals (special part of Coulomb matrix : contribution from single MT)
     334             : 
     335           0 :                      ALLOCATE (integral(n, n), carr(n, hybrid%nbands(nk)), carr2(n, lapw%nv(jsp)), carr3(n, lapw%nv(jsp)))
     336             : 
     337           0 :                      DO i = 1, n
     338           0 :                         CALL primitivef(primf1, fprod(:atoms%jri(itype), i)*atoms%rmsh(:atoms%jri(itype), itype)**(l + 1), atoms%rmsh, atoms%dx, atoms%jri, atoms%jmtd, itype, atoms%ntype)
     339           0 :                         CALL primitivef(primf2, fprod(:atoms%jri(itype), i)/atoms%rmsh(:atoms%jri(itype), itype)**l, atoms%rmsh, atoms%dx, atoms%jri, atoms%jmtd, -itype, atoms%ntype)  ! -itype is to enforce inward integration
     340             : 
     341           0 :                         primf1(:atoms%jri(itype)) = primf1(:atoms%jri(itype))/atoms%rmsh(:atoms%jri(itype), itype)**l
     342           0 :                         primf2(:atoms%jri(itype)) = primf2(:atoms%jri(itype))*atoms%rmsh(:atoms%jri(itype), itype)**(l + 1)
     343           0 :                         DO j = 1, n
     344           0 :                            integrand = fprod(:, j)*(primf1 + primf2)
     345             :                            integral(i, j) = fpi_const/(2*l + 1)*intgrf(integrand, atoms%jri, atoms%jmtd, atoms%rmsh, &
     346           0 :                                                                        atoms%dx, atoms%ntype, itype, hybdat%gridf)
     347             :                         END DO
     348             :                      END DO
     349             : 
     350             :                      ! Add everything up
     351             : 
     352           0 :                      DO m1 = -l1, l1
     353           0 :                         DO M = -l, l
     354           0 :                            m2 = m1 + M
     355             : 
     356           0 :                            carr = 0
     357           0 :                            DO n1 = 1, hybrid%nbands(nk)
     358             : 
     359           0 :                               DO i = 1, n
     360           0 :                                  ll = larr(i)
     361           0 :                                  IF (ABS(m2) > ll) CYCLE
     362             : 
     363             :                                  lm = SUM((/((2*l2 + 1)*hybrid%nindx(l2, itype), l2=0, ll - 1)/)) &
     364           0 :                                       + (m2 + ll)*hybrid%nindx(ll, itype) + parr(i)
     365             : 
     366           0 :                                  carr(i, n1) = cmt(n1, lm, iatom)*gaunt(l1, ll, l, m1, m2, M, hybdat%maxfac, hybdat%fac, hybdat%sfac)
     367             : 
     368             :                               END DO
     369           0 :                               DO n2 = 1, nsest(n1)!n1
     370           0 :                                  nn2 = indx_sest(n2, n1)
     371           0 :                                  exchange(nn2, n1) = exchange(nn2, n1) + dot_PRODUCT(carr(:, n1), MATMUL(integral, carr(:, nn2)))
     372             : 
     373             :                               END DO
     374             :                            END DO
     375             :                         END DO
     376             :                      END DO
     377             : 
     378           0 :                      DEALLOCATE (integral, carr, carr2, carr3)
     379             : 
     380             :                   END DO
     381             :                END DO
     382             :             END DO
     383             :          END DO
     384             :       END DO
     385             : 
     386           0 :       IF (mat_ex%l_real) THEN
     387           0 :          IF (ANY(ABS(AIMAG(exchange)) > 1e-10)) THEN
     388           0 :             IF (mpi%irank == 0) WRITE (6, '(A)') 'exchangeCore: Warning! Unusually large imaginary component.'
     389           0 :             WRITE (*, *) MAXVAL(ABS(AIMAG(exchange)))
     390           0 :             STOP 'exchangeCore: Unusually large imaginary component.'
     391             :          END IF
     392             :       ENDIF
     393             : 
     394           0 :       DO n1 = 1, hybrid%nobd(nk)
     395           0 :          results%te_hfex%core = results%te_hfex%Core - a_ex*results%w_iks(n1, nk, jsp)*exchange(n1, n1)
     396             :       END DO
     397             : 
     398             :       ! add the core-valence contribution to the exchange matrix mat_ex
     399             :       ! factor 1/nsymop is needed due to the symmetrization in symmetrizeh
     400             : 
     401           0 :       ic = 0
     402           0 :       sum_offdia = 0
     403           0 :       IF (mat_ex%l_real) THEN
     404           0 :          mat_ex%data_r = mat_ex%data_r + exchange/nsymop
     405             :       ELSE
     406           0 :          mat_ex%data_c = mat_ex%data_c + CONJG(exchange)/nsymop
     407             :       END IF
     408             : 
     409           0 :    END SUBROUTINE exchange_vccv1
     410             : 
     411           0 :    SUBROUTINE exchange_cccc(nk, atoms, hybdat, ncstd, sym, kpts, a_ex, mpi, results)
     412             : 
     413             :       USE m_constants
     414             :       USE m_util
     415             :       USE m_wrapper
     416             :       USE m_gaunt
     417             :       USE m_trafo
     418             :       USE m_types
     419             :       USE m_io_hybrid
     420             :       IMPLICIT NONE
     421             :       TYPE(t_hybdat), INTENT(IN)   :: hybdat
     422             :       TYPE(t_results), INTENT(INOUT)   :: results
     423             :       TYPE(t_mpi), INTENT(IN)   :: mpi
     424             :       TYPE(t_sym), INTENT(IN)   :: sym
     425             :       TYPE(t_kpts), INTENT(IN)   :: kpts
     426             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     427             : 
     428             :       ! - scalars -
     429             :       INTEGER, INTENT(IN)    ::  nk, ncstd
     430             : 
     431             :       REAL, INTENT(IN)    ::  a_ex
     432             : 
     433             :       ! - arays -
     434             : 
     435             :       ! - local scalars -
     436             :       INTEGER               ::  itype, ieq, icst, icst1, icst2, iatom, iatom0
     437             :       INTEGER               ::  l1, l2, l, ll, llmax
     438             :       INTEGER               ::  m1, m2, mm, m
     439             :       INTEGER               ::  n1, n2, n
     440             : 
     441             :       REAL                  ::  rdum, rdum1
     442             :       ! - local arrays -
     443           0 :       INTEGER               ::  point(hybdat%maxindxc, -hybdat%lmaxcd:hybdat%lmaxcd, 0:hybdat%lmaxcd, atoms%nat)
     444           0 :       REAL                  ::  rprod(atoms%jmtd), primf1(atoms%jmtd), primf2(atoms%jmtd), integrand(atoms%jmtd)
     445           0 :       COMPLEX               ::  exch(ncstd, ncstd)
     446             : 
     447             :       !       IF ( irank == 0 ) THEN
     448             :       !         WRITE(6,'(//A)') '### core-core-core-core exchange ###'
     449             :       !         WRITE(6,'(/A)') '        k-point       band    exchange'
     450             :       !       END IF
     451             : 
     452             :       ! set up point
     453           0 :       icst = 0
     454           0 :       iatom = 0
     455           0 :       DO itype = 1, atoms%ntype
     456           0 :          DO ieq = 1, atoms%neq(itype)
     457           0 :             iatom = iatom + 1
     458           0 :             DO l = 0, hybdat%lmaxc(itype)
     459           0 :                DO M = -l, l
     460           0 :                   DO n = 1, hybdat%nindxc(l, itype)
     461           0 :                      icst = icst + 1
     462           0 :                      point(n, M, l, iatom) = icst
     463             :                   END DO
     464             :                END DO
     465             :             END DO
     466             :          END DO
     467             :       END DO
     468             : 
     469           0 :       llmax = 2*hybdat%lmaxcd
     470           0 :       exch = 0
     471           0 :       iatom0 = 0
     472           0 :       DO itype = 1, atoms%ntype
     473             : 
     474           0 :          DO l1 = 0, hybdat%lmaxc(itype)  ! left core state
     475           0 :             DO l2 = 0, hybdat%lmaxc(itype)  ! right core state
     476           0 :                DO l = 0, hybdat%lmaxc(itype)   ! occupied core state
     477             : 
     478           0 :                   DO ll = ABS(l1 - l), l1 + l
     479           0 :                      IF (ll < ABS(l - l2) .OR. ll > l + l2) CYCLE
     480           0 :                      IF (MOD(l + l1 + ll, 2) /= 0) CYCLE
     481           0 :                      IF (MOD(l + l2 + ll, 2) /= 0) CYCLE
     482             : 
     483           0 :                      DO m1 = -l1, l1
     484           0 :                         m2 = m1
     485           0 :                         IF (ABS(m2) > l2) CYCLE
     486           0 :                         DO M = -l, l
     487           0 :                            mm = M - m1
     488           0 :                            IF (ABS(mm) > ll) CYCLE
     489           0 :                            rdum = fpi_const/(2*ll + 1)*gaunt1(l, ll, l1, M, mm, m1, llmax)*gaunt1(l, ll, l2, M, mm, m2, llmax)
     490             : 
     491           0 :                            DO n = 1, hybdat%nindxc(l, itype)
     492           0 :                               DO n2 = 1, hybdat%nindxc(l2, itype)
     493             :                                  rprod(:atoms%jri(itype)) = (hybdat%core1(:atoms%jri(itype), n, l, itype)*hybdat%core1(:atoms%jri(itype), n2, l2, itype) &
     494           0 :                                                              + hybdat%core2(:atoms%jri(itype), n, l, itype)*hybdat%core2(:atoms%jri(itype), n2, l2, itype))/atoms%rmsh(:atoms%jri(itype), itype)
     495             : 
     496           0 :                                  CALL primitivef(primf1, rprod(:)*atoms%rmsh(:, itype)**(ll + 1), atoms%rmsh, atoms%dx, atoms%jri, atoms%jmtd, itype, atoms%ntype)
     497           0 :                                  CALL primitivef(primf2, rprod(:atoms%jri(itype))/atoms%rmsh(:atoms%jri(itype), itype)**ll, atoms%rmsh, atoms%dx, atoms%jri, atoms%jmtd, -itype, atoms%ntype)  ! -itype is to enforce inward integration
     498             : 
     499           0 :                                  primf1 = primf1/atoms%rmsh(:, itype)**ll
     500           0 :                                  primf2 = primf2*atoms%rmsh(:, itype)**(ll + 1)
     501             : 
     502           0 :                                  DO n1 = 1, hybdat%nindxc(l1, itype)
     503             : 
     504             :                                     rprod(:atoms%jri(itype)) = (hybdat%core1(:atoms%jri(itype), n, l, itype)*hybdat%core1(:atoms%jri(itype), n1, l1, itype) &
     505           0 :                                                                 + hybdat%core2(:atoms%jri(itype), n, l, itype)*hybdat%core2(:atoms%jri(itype), n1, l1, itype))/atoms%rmsh(:atoms%jri(itype), itype)
     506             : 
     507           0 :                                     integrand = rprod*(primf1 + primf2)
     508             : 
     509             :                                     rdum1 = rdum*intgrf(integrand, atoms%jri, atoms%jmtd, &
     510           0 :                                                         atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
     511             : 
     512           0 :                                     iatom = iatom0
     513           0 :                                     DO ieq = 1, atoms%neq(itype)
     514           0 :                                        iatom = iatom + 1
     515           0 :                                        icst1 = point(n1, m1, l1, iatom)
     516           0 :                                        icst2 = point(n2, m2, l2, iatom)
     517           0 :                                        exch(icst1, icst2) = exch(icst1, icst2) + rdum1
     518             :                                     END DO
     519             :                                  END DO  !n1
     520             : 
     521             :                               END DO  !n2
     522             :                            END DO  !n
     523             : 
     524             :                         END DO  !M
     525             :                      END DO  !m1
     526             : 
     527             :                   END DO  !ll
     528             : 
     529             :                END DO  !l
     530             :             END DO  !l2
     531             :          END DO  !l1
     532           0 :          iatom0 = iatom0 + atoms%neq(itype)
     533             :       END DO  !itype
     534             : 
     535           0 :       IF (sym%invs) THEN
     536           0 :          CALL symmetrize(exch, ncstd, ncstd, 3, .FALSE., atoms, hybdat%lmaxc, hybdat%lmaxcd, hybdat%nindxc, sym)
     537           0 :          IF (ANY(ABS(AIMAG(exch)) > 1E-6)) STOP 'exchange_cccc: exch possesses significant imaginary part'
     538             :       ENDIF
     539             :       !       DO icst = 1,ncstd
     540             :       !         IF ( irank == 0 )
     541             :       !      &    WRITE(6,'(    ''  ('',F5.3,'','',F5.3,'','',F5.3,'')'',I4,1X,F12.5)')bkpt,icst,REAL(exch(icst,icst))*(-27.211608)
     542             :       !       END DO
     543             : 
     544             :       ! add core exchange contributions to the te_hfex
     545             : 
     546           0 :       DO icst1 = 1, ncstd
     547           0 :          results%te_hfex%core = results%te_hfex%core - a_ex*kpts%wtkpt(nk)*exch(icst1, icst1)
     548             :       END DO
     549             : 
     550           0 :    END SUBROUTINE exchange_cccc
     551             : 
     552           0 :    SUBROUTINE exchange_cccv(nk, atoms, hybdat, hybrid, DIMENSION, maxbands, ncstd, &
     553           0 :                             bkpt, sym, mpi, exch_cv_r, exch_cv_c, l_real)
     554             : 
     555             :       USE m_constants
     556             :       USE m_util
     557             :       USE m_wrapper
     558             :       USE m_gaunt
     559             :       USE m_trafo
     560             :       USE m_io_hybrid
     561             :       USE m_types
     562             :       IMPLICIT NONE
     563             :       TYPE(t_hybdat), INTENT(IN)   :: hybdat
     564             :       TYPE(t_mpi), INTENT(IN)   :: mpi
     565             :       TYPE(t_dimension), INTENT(IN)   :: DIMENSION
     566             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
     567             :       TYPE(t_sym), INTENT(IN)   :: sym
     568             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     569             :       ! - scalars -
     570             :       INTEGER, INTENT(IN)    ::  nk, ncstd
     571             :       INTEGER, INTENT(IN)    :: maxbands
     572             : 
     573             :       ! - arays -
     574             :       REAL, INTENT(IN)    ::  bkpt(3)
     575             :       LOGICAL, INTENT(IN)    :: l_real
     576             :       REAL, INTENT(INOUT) ::  exch_cv_r(:, :, :)!(maxbands,ncstd,nkpti)
     577             :       COMPLEX, INTENT(INOUT) ::  exch_cv_c(:, :, :) !(maxbands,ncstd,nkpti)
     578             :       ! - local scalars -
     579             :       INTEGER               ::  itype, ieq, icst, icst1, icst2, iatom, iatom0
     580             :       INTEGER               ::    iatom1, iband
     581             :       INTEGER               ::  l1, l2, l, ll, llmax
     582             :       INTEGER               ::  lm2, lmp2
     583             :       INTEGER               ::  m1, m2, mm, m
     584             :       INTEGER               ::  n1, n2, n, nn
     585             : 
     586             :       REAL                  ::  rdum0, rdum1, rdum2, rdum3, rdum4
     587             :       COMPLEX               ::  cdum
     588             :       COMPLEX, PARAMETER     ::  img = (0.0, 1.0)
     589             :       ! - local arrays -
     590           0 :       INTEGER               ::  point(hybdat%maxindxc, -hybdat%lmaxcd:hybdat%lmaxcd, 0:hybdat%lmaxcd, atoms%nat)
     591           0 :       INTEGER               ::  lmstart(0:atoms%lmaxd, atoms%ntype)
     592           0 :       REAL                  ::  rprod(atoms%jmtd), primf1(atoms%jmtd), primf2(atoms%jmtd)
     593           0 :       REAL                  ::  integrand(atoms%jmtd)
     594           0 :       COMPLEX               ::  cexp(atoms%nat)
     595           0 :       COMPLEX               ::  exch(hybrid%nbands(nk), ncstd)
     596           0 :       COMPLEX               ::  cmt(DIMENSION%neigd, hybrid%maxlmindx, atoms%nat), carr(hybrid%nbands(nk))
     597             : 
     598           0 :       IF (mpi%irank == 0) THEN
     599           0 :          WRITE (6, '(//A)') '### core-core-core-valence exchange  ###'
     600           0 :          WRITE (6, '(/A)') '        k-point       band    exchange'
     601             :       END IF
     602             : 
     603             :       ! set up point
     604           0 :       icst = 0
     605           0 :       iatom = 0
     606           0 :       DO itype = 1, atoms%ntype
     607           0 :          DO ieq = 1, atoms%neq(itype)
     608           0 :             iatom = iatom + 1
     609           0 :             DO l = 0, hybdat%lmaxc(itype)
     610           0 :                DO M = -l, l
     611           0 :                   DO n = 1, hybdat%nindxc(l, itype)
     612           0 :                      icst = icst + 1
     613           0 :                      point(n, M, l, iatom) = icst
     614             :                   END DO
     615             :                END DO
     616             :             END DO
     617             :          END DO
     618             :       END DO
     619             : 
     620             :       ! lmstart = lm start index for each l-quantum number and atom type (for cmt-coefficients)
     621           0 :       DO itype = 1, atoms%ntype
     622           0 :          DO l = 0, atoms%lmax(itype)
     623           0 :             lmstart(l, itype) = SUM((/(hybrid%nindx(ll, itype)*(2*ll + 1), ll=0, l - 1)/))
     624             :          END DO
     625             :       END DO
     626             : 
     627             :       ! read in cmt coefficient at k-point nk
     628             : 
     629           0 :       CALL read_cmt(cmt, nk)
     630           0 :       iatom = 0
     631           0 :       DO itype = 1, atoms%ntype
     632           0 :          DO ieq = 1, atoms%neq(itype)
     633           0 :             iatom = iatom + 1
     634           0 :             cexp(iatom) = EXP(img*tpi_const*dot_PRODUCT(bkpt(:), atoms%taual(:, iatom)))
     635             :          END DO
     636             :       END DO
     637             : 
     638           0 :       cmt = CONJG(cmt)
     639             : 
     640           0 :       llmax = MAX(2*hybdat%lmaxcd, atoms%lmaxd)
     641             : 
     642           0 :       exch = 0
     643           0 :       iatom0 = 0
     644           0 :       DO itype = 1, atoms%ntype
     645             : 
     646           0 :          DO l1 = 0, hybdat%lmaxc(itype)  ! left core state
     647           0 :             DO l2 = 0, atoms%lmax(itype)  ! right valence state
     648           0 :                DO l = 0, hybdat%lmaxc(itype)   ! occupied core state
     649             : 
     650           0 :                   DO ll = ABS(l1 - l), l1 + l
     651           0 :                      IF (ll < ABS(l - l2) .OR. ll > l + l2) CYCLE
     652           0 :                      IF (MOD(l + l1 + ll, 2) /= 0) CYCLE
     653           0 :                      IF (MOD(l + l2 + ll, 2) /= 0) CYCLE
     654             : 
     655             :                      !                 WRITE(*,*) 'l1,l2,l,ll',l1,l2,l,ll
     656           0 :                      rdum0 = fpi_const/(2*ll + 1)
     657             : 
     658           0 :                      DO m1 = -l1, l1
     659           0 :                         m2 = m1
     660           0 :                         IF (ABS(m2) > l2) CYCLE
     661           0 :                         lm2 = lmstart(l2, itype) + (m2 + l2)*hybrid%nindx(l2, itype)
     662             : 
     663           0 :                         DO M = -l, l
     664           0 :                            mm = M - m1
     665           0 :                            IF (ABS(M - m1) > ll) CYCLE
     666             : 
     667             :                            rdum1 = gaunt1(l, ll, l1, M, mm, m1, llmax) &
     668           0 :                                    *gaunt1(l, ll, l2, M, mm, m1, llmax)*rdum0
     669             : 
     670           0 :                            DO n = 1, hybdat%nindxc(l, itype)
     671           0 :                               DO n2 = 1, hybrid%nindx(l2, itype)
     672           0 :                                  lmp2 = lm2 + n2
     673             : 
     674             :                                  rprod(:atoms%jri(itype)) = (hybdat%core1(:atoms%jri(itype), n, l, itype)*hybdat%bas1(:atoms%jri(itype), n2, l2, itype) &
     675           0 :                                                              + hybdat%core2(:atoms%jri(itype), n, l, itype)*hybdat%bas2(:atoms%jri(itype), n2, l2, itype))/atoms%rmsh(:atoms%jri(itype), itype)
     676             : 
     677           0 :                                  CALL primitivef(primf1, rprod(:atoms%jri(itype))*atoms%rmsh(:atoms%jri(itype), itype)**(ll + 1), atoms%rmsh, atoms%dx, atoms%jri, atoms%jmtd, itype, atoms%ntype)
     678           0 :                                  CALL primitivef(primf2, rprod(:atoms%jri(itype))/atoms%rmsh(:atoms%jri(itype), itype)**ll, atoms%rmsh, atoms%dx, atoms%jri, atoms%jmtd, -itype, atoms%ntype)  ! -itype is to enforce inward integration
     679             : 
     680           0 :                                  primf1(:atoms%jri(itype)) = primf1(:atoms%jri(itype))/atoms%rmsh(:atoms%jri(itype), itype)**ll
     681           0 :                                  primf2 = primf2*atoms%rmsh(:, itype)**(ll + 1)
     682             : 
     683           0 :                                  DO n1 = 1, hybdat%nindxc(l1, itype)
     684             : 
     685             :                                     rprod(:) = (hybdat%core1(:, n, l, itype)*hybdat%core1(:, n1, l1, itype) &
     686           0 :                                                 + hybdat%core2(:, n, l, itype)*hybdat%core2(:, n1, l1, itype))/atoms%rmsh(:atoms%jri(itype), itype)
     687             : 
     688           0 :                                     integrand = rprod*(primf1 + primf2)
     689             : 
     690           0 :                                     rdum2 = rdum1*intgrf(integrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
     691             : 
     692           0 :                                     iatom = iatom0
     693           0 :                                     DO ieq = 1, atoms%neq(itype)
     694           0 :                                        iatom = iatom + 1
     695           0 :                                        icst1 = point(n1, m1, l1, iatom)
     696           0 :                                        cdum = rdum2*cexp(iatom)
     697           0 :                                        DO iband = 1, hybrid%nbands(nk)
     698             : 
     699           0 :                                           exch(iband, icst1) = exch(iband, icst1) + cdum*cmt(iband, lmp2, iatom)
     700             : 
     701             :                                        END DO
     702             :                                     END DO
     703             : 
     704             :                                  END DO  !n1
     705             : 
     706             :                               END DO  !n2
     707             :                            END DO  !n
     708             : 
     709             :                         END DO  !M
     710             :                      END DO  !m1
     711             : 
     712             :                   END DO  !ll
     713             : 
     714             :                END DO  !l
     715             :             END DO  !l2
     716             :          END DO  !l1
     717           0 :          iatom0 = iatom0 + atoms%neq(itype)
     718             :       END DO  !itype
     719             : 
     720           0 :       IF (l_real) THEN
     721             :          !symmetrize core-wavefunctions such that phi(-r) = phi(r)*
     722           0 :          CALL symmetrize(exch, hybrid%nbands(nk), ncstd, 2, .FALSE., atoms, hybdat%lmaxc, hybdat%lmaxcd, hybdat%nindxc, sym)
     723             : 
     724           0 :          IF (ANY(ABS(AIMAG(exch)) > 1E-6)) STOP 'exchange_cccv: exch possesses significant imaginary part'
     725             :       ENDIF
     726             : 
     727           0 :       IF (l_real) THEN
     728           0 :          DO icst = 1, ncstd
     729           0 :             DO iband = 1, hybrid%nbands(nk)
     730           0 :                exch_cv_r(iband, icst, nk) = exch_cv_r(iband, icst, nk) - exch(iband, icst)
     731             :             END DO
     732             :          END DO
     733             :       ELSE
     734           0 :          DO icst = 1, ncstd
     735           0 :             DO iband = 1, hybrid%nbands(nk)
     736           0 :                exch_cv_c(iband, icst, nk) = exch_cv_c(iband, icst, nk) - exch(iband, icst)
     737             :             END DO
     738             :          END DO
     739             :       END IF
     740             : 
     741           0 :    END SUBROUTINE exchange_cccv
     742             : 
     743             : END MODULE m_exchange_core

Generated by: LCOV version 1.13