LCOV - code coverage report
Current view: top level - hybrid - exchange_val_hf.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 147 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 2 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             : !     Calculates the HF exchange term
       8             : !
       9             : !                                          s          s*          s            s*
      10             : !                                       phi    (r) phi     (r) phi     (r') phi    (r')
      11             : !                         occ.             n_1k       n'k+q       n'k+q        n_2k
      12             : !     exchange(n,q)  =  - SUM  INT INT  ------------------------------------------- dr dr'
      13             : !                         k,n'                           | r - r' |
      14             : !
      15             : !                         occ                  s          s    ~        ~       s         s
      16             : !                    =  - SUM  SUM  v     < phi      | phi     M    > < M    phi     | phi      >
      17             : !                         k,n' I,J   k,IJ      n'k+q      n_1k  q,I      q,J    n_2k      n'k+q
      18             : !
      19             : !     for the different combinations of n_1 and n_2 and where n' runs only over the valence states.
      20             : !     ( n_1,n_2:  valence-valence, core-core,core-valence )
      21             : !
      22             : !
      23             : !     At the Gamma point (k=0) v diverges. After diagonalization of v at k=0 the divergence is
      24             : !     restricted to the head element I=1. Furthermore, we expand <...> with kp perturbation theory.
      25             : !     As a result, the total I=1 element is given by a sum of a divergent 1/k**2-term and an
      26             : !     angular dependent term. The former is separated from the numerical k-summation and treated
      27             : !     analytically while the latter is spherically averaged and added to the k=0 contribution of
      28             : !     the numerical k-summation. (A better knowledge of the integrand's behavior at the BZ edges
      29             : !     might further improve the integration.)
      30             : !
      31             : !     The divergence at the Gamma point is integrated with one of the following algorithms:
      32             : ! (1) Switching-Off Function
      33             : !     In a sphere of radius k0=radshmin/2 a switching-off function g(k)=1-(k/k0)**n*(n+1-n*k/k0)
      34             : !     (n=npot) is defined. The 1/k**2 divergence is subtracted from the BZ integral in the form
      35             : !     g(k)/k**2 and integrated analytically. The non-divergent rest is integrated numerically.
      36             : ! (2) Periodic Function (similar to the one used by Massidda PRB 48, 5058)
      37             : !     The function  F(k) = SUM(G) exp(-expo*|k+G|**3) / |k+G|**2  is subtracted from the BZ integral
      38             : !     and integrated analytically. The non-divergent rest is integrated numerically.
      39             : !     The parameter expo is chosen such that exp(-expo*q**3)=1/2
      40             : !     with q = radius of sphere with same volume as BZ.
      41             : ! (3) Periodic Function (same as Massidda's) with expo->0
      42             : !     The function  F(k) = lim(expo->0) SUM(G) exp(-expo*|k+G|**2) / |k+G|**2  is subtracted from
      43             : !     the BZ integral and integrated analytically. The contribution to the BZ integral including
      44             : !     the "tail" is
      45             : !     vol/(8*pi**3) INT F(k) d^3k - P SUM(k) F(k)  ( P = principal value ) .
      46             : !     For expo->0 the two terms diverge. Therefore a cutoff radius q0 is introduced and related to
      47             : !     expo by exp(-expo*q0**2)=delta  ( delta = small value, e.g., delta = 1*10.0**-10 ) .
      48             : !     The resulting formula
      49             : !     vol/(4*pi**1.5*sqrt(expo)) * erf(sqrt(a)*q0) - sum(q,0<q<q0) exp(-expo*q**2)/q**2
      50             : !     converges well with q0. (Should be the default.)
      51             : 
      52             : MODULE m_exchange_valence_hf
      53             : 
      54             :    LOGICAL, PARAMETER:: zero_order = .false., ibs_corr = .false.
      55             :    INTEGER, PARAMETER:: maxmem = 600
      56             : 
      57             : CONTAINS
      58             : 
      59           0 :    SUBROUTINE exchange_valence_hf(nk, kpts, nkpt_EIBZ, sym, atoms, hybrid, cell, dimension, input, jsp, hybdat, mnobd, lapw, &
      60           0 :                                   eig_irr, results, parent, pointer_EIBZ, n_q, wl_iks, it, xcpot, noco, nsest, indx_sest, &
      61             :                                   mpi, mat_ex)
      62             : 
      63             :       USE m_types
      64             :       USE m_wrapper
      65             :       USE m_constants
      66             :       USE m_trafo
      67             :       USE m_wavefproducts
      68             :       USE m_olap
      69             :       USE m_spmvec
      70             :       USE m_hsefunctional, ONLY: dynamic_hse_adjustment
      71             :       USE m_io_hybrid
      72             :       USE m_kp_perturbation
      73             : 
      74             :       IMPLICIT NONE
      75             : 
      76             :       TYPE(t_results), INTENT(IN)    :: results
      77             :       TYPE(t_xcpot_inbuild), INTENT(IN)    :: xcpot
      78             :       TYPE(t_mpi), INTENT(IN)    :: mpi
      79             :       TYPE(t_dimension), INTENT(IN)    :: dimension
      80             :       TYPE(t_hybrid), INTENT(INOUT) :: hybrid
      81             :       TYPE(t_input), INTENT(IN)    :: input
      82             :       TYPE(t_noco), INTENT(IN)    :: noco
      83             :       TYPE(t_sym), INTENT(IN)    :: sym
      84             :       TYPE(t_cell), INTENT(IN)    :: cell
      85             :       TYPE(t_kpts), INTENT(IN)    :: kpts
      86             :       TYPE(t_atoms), INTENT(IN)    :: atoms
      87             :       TYPE(t_lapw), INTENT(IN)    :: lapw
      88             :       TYPE(t_mat), INTENT(INOUT) :: mat_ex
      89             :       TYPE(t_hybdat), INTENT(INOUT) :: hybdat
      90             : 
      91             :       ! scalars
      92             :       INTEGER, INTENT(IN)    :: it
      93             :       INTEGER, INTENT(IN)    :: jsp
      94             :       INTEGER, INTENT(IN)    :: nk, nkpt_EIBZ
      95             :       INTEGER, INTENT(IN)    :: mnobd
      96             : 
      97             :       ! arrays
      98             :       INTEGER, INTENT(IN)    ::  n_q(nkpt_EIBZ)
      99             : 
     100             :       INTEGER, INTENT(IN)    ::  parent(kpts%nkptf)
     101             :       INTEGER, INTENT(IN)    ::  pointer_EIBZ(nkpt_EIBZ)
     102             :       INTEGER, INTENT(IN)    ::  nsest(hybrid%nbands(nk))
     103             :       INTEGER, INTENT(IN)    ::  indx_sest(hybrid%nbands(nk), hybrid%nbands(nk))
     104             : 
     105             :       REAL, INTENT(IN)    ::  eig_irr(dimension%neigd, kpts%nkpt)
     106             :       REAL, INTENT(IN)    ::  wl_iks(dimension%neigd, kpts%nkptf)
     107             : 
     108             :       ! local scalars
     109             :       INTEGER                 ::  iband, iband1, ibando, ikpt, ikpt0
     110             :       INTEGER                 ::  i, ic, ix, iy, iz
     111             :       INTEGER                 ::  irecl_coulomb, irecl_coulomb1
     112             :       INTEGER                 ::  j
     113             :       INTEGER                 ::  m1, m2
     114             :       INTEGER                 ::  n, n1, n2, nn, nn2
     115             :       INTEGER                 ::  nkqpt
     116             :       INTEGER                 ::  npot
     117             :       INTEGER                 ::  ok
     118             :       INTEGER                 ::  psize
     119             :       REAL                    ::  rdum
     120             :       REAL                    ::  k0
     121             : 
     122             :       REAL, SAVE             ::  divergence
     123             : 
     124             :       COMPLEX                 ::  cdum, cdum1, cdum2
     125             :       COMPLEX                 ::  exch0
     126             : 
     127             :       LOGICAL, SAVE           ::  initialize = .true.
     128             : 
     129             :       ! local arrays
     130             :       INTEGER              :: kcorner(3, 8) = reshape((/0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1/), (/3, 8/))
     131           0 :       COMPLEX              :: exchcorrect(kpts%nkptf)
     132           0 :       COMPLEX              :: dcprod(hybrid%nbands(nk), hybrid%nbands(nk), 3)
     133           0 :       COMPLEX              :: exch_vv(hybrid%nbands(nk), hybrid%nbands(nk))
     134             :       COMPLEX              :: hessian(3, 3)
     135           0 :       COMPLEX              :: proj_ibsc(3, mnobd, hybrid%nbands(nk))
     136           0 :       COMPLEX              :: olap_ibsc(3, 3, mnobd, mnobd)
     137           0 :       REAL                 :: carr1_v_r(hybrid%maxbasm1), carr1_c_r(hybrid%maxbasm1)
     138           0 :       COMPLEX              :: carr1_v_c(hybrid%maxbasm1), carr1_c_c(hybrid%maxbasm1)
     139           0 :       COMPLEX, ALLOCATABLE :: phase_vv(:, :)
     140           0 :       REAL, ALLOCATABLE :: cprod_vv_r(:, :, :), cprod_cv_r(:, :, :), carr3_vv_r(:, :, :), carr3_cv_r(:, :, :)
     141           0 :       COMPLEX, ALLOCATABLE :: cprod_vv_c(:, :, :), cprod_cv_c(:, :, :), carr3_vv_c(:, :, :), carr3_cv_c(:, :, :)
     142             : 
     143             : 
     144             : #if ( !defined CPP_NOSPMVEC && !defined CPP_IRAPPROX )
     145           0 :       REAL                 :: coulomb_mt1(hybrid%maxindxm1 - 1, hybrid%maxindxm1 - 1, 0:hybrid%maxlcutm1, atoms%ntype)
     146           0 :       REAL                 :: coulomb_mt2_r(hybrid%maxindxm1 - 1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1 + 1, atoms%nat)
     147           0 :       REAL                 :: coulomb_mt3_r(hybrid%maxindxm1 - 1, atoms%nat, atoms%nat)
     148           0 :       COMPLEX              :: coulomb_mt2_c(hybrid%maxindxm1 - 1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1 + 1, atoms%nat)
     149           0 :       COMPLEX              :: coulomb_mt3_c(hybrid%maxindxm1 - 1, atoms%nat, atoms%nat)
     150             : #else
     151             :       REAL                 :: coulomb_r(hybrid%maxbasm1*(hybrid%maxbasm1 + 1)/2)
     152             :       COMPLEX              :: coulomb_c(hybrid%maxbasm1*(hybrid%maxbasm1 + 1)/2)
     153             : #endif
     154             : 
     155             : #ifdef CPP_IRCOULOMBAPPROX
     156             :       REAL                 :: coulomb_mtir_r((hybrid%maxlcutm1 + 1)**2*atoms%nat, &
     157             :                                              (hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(hybrid%ngptm))
     158             : #else
     159             :       REAL                 :: coulomb_mtir_r(((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(hybrid%ngptm))* &
     160           0 :                                              ((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(hybrid%ngptm) + 1)/2)
     161             : #endif
     162             : 
     163             : #ifdef CPP_IRCOULOMBAPPROX
     164             :       COMPLEX              :: coulomb_mtir_c((hybrid%maxlcutm1 + 1)**2*atoms%nat, &
     165             :                                              (hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(hybrid%ngptm))
     166             : #else
     167             :       COMPLEX              :: coulomb_mtir_c(((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(hybrid%ngptm))* &
     168           0 :                                              ((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(hybrid%ngptm) + 1)/2)
     169             : #endif
     170             : 
     171           0 :       LOGICAL              :: occup(dimension%neigd)
     172           0 :       CALL timestart("valence exchange calculation")
     173             : 
     174           0 :       IF (initialize) THEN !it .eq. 1 .and. nk .eq. 1) THEN
     175           0 :          call calc_divergence(cell, kpts, divergence)
     176           0 :          PRINT *, "Divergence:", divergence
     177           0 :          initialize = .false.
     178             :       END IF
     179             : 
     180             :       ! calculate valence-valence-valence-valence, core-valence-valence-valence
     181             :       ! and core-valence-valence-core exchange at current k-point
     182             :       ! the sum over the inner occupied valence states is restricted to the EIBZ(k)
     183             :       ! the contribution of the Gamma-point is treated separately (see below)
     184             : 
     185             :       ! determine package size loop over the occupied bands
     186           0 :       rdum = hybrid%maxbasm1*hybrid%nbands(nk)*4/1048576.
     187           0 :       psize = 1
     188           0 :       DO iband = mnobd, 1, -1
     189             :          ! ensure that the packages have equal size
     190           0 :          IF (modulo(mnobd, iband) == 0) THEN
     191             :             ! choose packet size such that cprod is smaller than memory threshold
     192           0 :             IF (rdum*iband <= maxmem) THEN
     193           0 :                psize = iband
     194           0 :                EXIT
     195             :             END IF
     196             :          END IF
     197             :       END DO
     198             : 
     199           0 :       IF (psize /= mnobd) THEN
     200           0 :          WRITE (6, '(A,A,i3,A,f7.2,A)') ' Divide the loop over the occupied hybrid%bands in packages', &
     201           0 :             ' of the size', psize, ' (cprod=', rdum*psize, 'MB)'
     202             :       END IF
     203           0 :       ALLOCATE (phase_vv(psize, hybrid%nbands(nk)), stat=ok)
     204           0 :       IF (ok /= 0) STOP 'exchange_val_hf: error allocation phase'
     205           0 :       phase_vv = 0
     206             :       IF (ok /= 0) STOP 'exchange_val_hf: error allocation phase'
     207             : 
     208           0 :       if (mat_ex%l_real) THEN
     209           0 :          ALLOCATE (cprod_vv_c(hybrid%maxbasm1, 0, 0), carr3_vv_c(hybrid%maxbasm1, 0, 0))
     210           0 :          ALLOCATE (cprod_vv_r(hybrid%maxbasm1, psize, hybrid%nbands(nk)), stat=ok)
     211           0 :          IF (ok /= 0) STOP 'exchange_val_hf: error allocation cprod'
     212           0 :          ALLOCATE (carr3_vv_r(hybrid%maxbasm1, psize, hybrid%nbands(nk)), stat=ok)
     213           0 :          IF (ok /= 0) STOP 'exchange_val_hf: error allocation carr3'
     214           0 :          cprod_vv_r = 0; carr3_vv_r = 0
     215             :       ELSE
     216           0 :          ALLOCATE (cprod_vv_r(hybrid%maxbasm1, 0, 0), carr3_vv_r(hybrid%maxbasm1, 0, 0))
     217           0 :          ALLOCATE (cprod_vv_c(hybrid%maxbasm1, psize, hybrid%nbands(nk)), stat=ok)
     218           0 :          IF (ok /= 0) STOP 'exchange_val_hf: error allocation cprod'
     219           0 :          ALLOCATE (carr3_vv_c(hybrid%maxbasm1, psize, hybrid%nbands(nk)), stat=ok)
     220           0 :          IF (ok /= 0) STOP 'exchange_val_hf: error allocation carr3'
     221           0 :          cprod_vv_c = 0; carr3_vv_c = 0
     222             :       END IF
     223             : 
     224           0 :       exch_vv = 0
     225             : 
     226           0 :       DO ikpt = 1, nkpt_EIBZ
     227             : 
     228           0 :          ikpt0 = pointer_EIBZ(ikpt)
     229             : 
     230           0 :          n = hybrid%nbasp + hybrid%ngptm(ikpt0)
     231           0 :          IF (hybrid%nbasm(ikpt0) /= n) STOP 'error hybrid%nbasm'
     232           0 :          nn = n*(n + 1)/2
     233             : 
     234             :          ! read in coulomb matrix from direct access file coulomb
     235           0 :          IF (mat_ex%l_real) THEN
     236           0 :             CALL read_coulomb_spm_r(kpts%bkp(ikpt0), coulomb_mt1, coulomb_mt2_r, coulomb_mt3_r, coulomb_mtir_r)
     237             :          ELSE
     238           0 :             CALL read_coulomb_spm_c(kpts%bkp(ikpt0), coulomb_mt1, coulomb_mt2_c, coulomb_mt3_c, coulomb_mtir_c)
     239             :          END IF
     240             : 
     241           0 :          IF (kpts%bkp(ikpt0) /= ikpt0) THEN
     242             : #if( !defined CPP_NOSPMVEC && !defined CPP_IRAPPROX )
     243           0 :             IF ((kpts%bksym(ikpt0) > sym%nop) .and. (.not. mat_ex%l_real)) THEN
     244           0 :                coulomb_mt2_c = conjg(coulomb_mt2_c)
     245           0 :                coulomb_mtir_c = conjg(coulomb_mtir_c)
     246             :             END IF
     247             : #else
     248             :             if (.not. mat_ex%l_real) THEN
     249             :                IF (kpts%bksym(ikpt0) > sym%nop) coulomb = conjg(coulomb)
     250             :             endif
     251             : #endif
     252             :          END IF
     253             : 
     254           0 :          DO ibando = 1, mnobd, psize
     255             : 
     256           0 :             IF (mat_ex%l_real) THEN
     257             : #ifdef CPP_IRAPPROX
     258             :                CALL wavefproducts_inv(1, hybdat, dimension, input, jsp, atoms, lapw, obsolete, kpts, nk, ikpt0, &
     259             :                                       mnobd, hybrid, parent, cell, sym, noco, nkqpt, cprod_vv)
     260             : #else
     261             :                CALL wavefproducts_inv5(1, hybrid%nbands(nk), ibando, ibando + psize - 1, dimension, input, jsp, atoms, &
     262             :                                        lapw, kpts, nk, ikpt0, hybdat, mnobd, hybrid, parent, cell, hybrid%nbasp, sym, &
     263           0 :                                        noco, nkqpt, cprod_vv_r)
     264             : #endif
     265             :             ELSE
     266             : #ifdef CPP_IRAPPROX
     267             :                CALL wavefproducts_noinv(1, hybdat, nk, ikpt0, dimension, input, jsp, cell, atoms, hybrid,
     268             :                kpts, mnobd, lapw, sym, noco, nkqpt, cprod_vv)
     269             : #else
     270             :                CALL wavefproducts_noinv5(1, hybrid%nbands(nk), ibando, ibando + psize - 1, nk, ikpt0, dimension, input, jsp, &!jsp,&
     271           0 :                                          cell, atoms, hybrid, hybdat, kpts, mnobd, lapw, sym, hybrid%nbasp, noco, nkqpt, cprod_vv_c)
     272             : #endif
     273             :             END IF
     274             : 
     275             :             ! The sparse matrix technique is not feasible for the HSE
     276             :             ! functional. Thus, a dynamic adjustment is implemented
     277             :             ! The mixed basis functions and the potential difference
     278             :             ! are Fourier transformed, so that the exchange can be calculated
     279             :             ! in Fourier space
     280             : #ifndef CPP_NOSPMVEC
     281           0 :             IF (xcpot%is_name("hse") .OR. xcpot%is_name("vhse")) THEN
     282           0 :                iband1 = hybrid%nobd(nkqpt)
     283             : 
     284             :                exch_vv = exch_vv + &
     285             :                          dynamic_hse_adjustment(atoms%rmsh, atoms%rmt, atoms%dx, atoms%jri, atoms%jmtd, kpts%bkf(:, ikpt0), ikpt0, &
     286             :                                                 kpts%nkptf, cell%bmat, cell%omtil, atoms%ntype, atoms%neq, atoms%nat, atoms%taual, &
     287             :                                                 hybrid%lcutm1, hybrid%maxlcutm1, hybrid%nindxm1, hybrid%maxindxm1, hybrid%gptm, &
     288             :                                                 hybrid%ngptm(ikpt0), hybrid%pgptm(:, ikpt0), hybrid%gptmd, hybrid%basm1, &
     289             :                                                 hybrid%nbasm(ikpt0), iband1, hybrid%nbands(nk), nsest, ibando, psize, indx_sest, &
     290             :                                                 atoms%invsat, sym%invsatnr, mpi%irank, cprod_vv_r(:hybrid%nbasm(ikpt0), :, :), &
     291           0 :                                                 cprod_vv_c(:hybrid%nbasm(ikpt0), :, :), mat_ex%l_real, wl_iks(:iband1, nkqpt), n_q(ikpt))
     292             :             END IF
     293             : #endif
     294             : 
     295             :             ! the Coulomb matrix is only evaluated at the irrecuible k-points
     296             :             ! bra_trafo transforms cprod instead of rotating the Coulomb matrix
     297             :             ! from IBZ to current k-point
     298           0 :             IF (kpts%bkp(ikpt0) /= ikpt0) THEN
     299             :                CALL bra_trafo2(mat_ex%l_real, carr3_vv_r(:hybrid%nbasm(ikpt0), :, :), cprod_vv_r(:hybrid%nbasm(ikpt0), :, :), &
     300             :                                carr3_vv_c(:hybrid%nbasm(ikpt0), :, :), cprod_vv_c(:hybrid%nbasm(ikpt0), :, :), &
     301             :                                hybrid%nbasm(ikpt0), psize, hybrid%nbands(nk), kpts%bkp(ikpt0), ikpt0, kpts%bksym(ikpt0), sym, &
     302           0 :                                hybrid, kpts, cell, atoms, phase_vv)
     303           0 :                IF (mat_ex%l_real) THEN
     304           0 :                   cprod_vv_r(:hybrid%nbasm(ikpt0), :, :) = carr3_vv_r(:hybrid%nbasm(ikpt0), :, :)
     305             :                ELSE
     306           0 :                   cprod_vv_c(:hybrid%nbasm(ikpt0), :, :) = carr3_vv_c(:hybrid%nbasm(ikpt0), :, :)
     307             :                ENDIF
     308             :             ELSE
     309           0 :                phase_vv(:, :) = (1.0, 0.0)
     310             :             END IF
     311             : 
     312             :             ! calculate exchange matrix at ikpt0
     313             : 
     314           0 :             call timestart("exchange matrix")
     315           0 :             DO n1 = 1, hybrid%nbands(nk)
     316           0 :                DO iband = 1, psize
     317           0 :                   IF ((ibando + iband - 1) > hybrid%nobd(nkqpt)) CYCLE
     318             : 
     319           0 :                   cdum = wl_iks(ibando + iband - 1, nkqpt)*conjg(phase_vv(iband, n1))/n_q(ikpt)
     320             : 
     321           0 :                   IF (mat_ex%l_real) THEN
     322           0 :                      carr1_v_r(:n) = 0
     323             :                      CALL spmvec_invs(atoms, hybrid, hybdat, ikpt0, kpts, cell, coulomb_mt1, coulomb_mt2_r, coulomb_mt3_r, &
     324           0 :                                       coulomb_mtir_r, cprod_vv_r(:n, iband, n1), carr1_v_r(:n))
     325             :                   ELSE
     326           0 :                      carr1_v_c(:n) = 0
     327             :                      CALL spmvec_noinvs(atoms, hybrid, hybdat, ikpt0, kpts, cell, coulomb_mt1, coulomb_mt2_c, coulomb_mt3_c, &
     328           0 :                                         coulomb_mtir_c, cprod_vv_c(:n, iband, n1), carr1_v_c(:n))
     329             :                   END IF
     330             : 
     331           0 :                   IF (mat_ex%l_real) THEN
     332           0 :                      DO n2 = 1, nsest(n1)!n1
     333           0 :                         nn2 = indx_sest(n2, n1)
     334             :                         exch_vv(nn2, n1) = exch_vv(nn2, n1) + cdum*phase_vv(iband, nn2)* &
     335           0 :                                            dotprod(carr1_v_r(:n), cprod_vv_r(:n, iband, nn2))
     336             :                      END DO !n2
     337             :                   ELSE
     338           0 :                      DO n2 = 1, nsest(n1)!n1
     339           0 :                         nn2 = indx_sest(n2, n1)
     340             :                         exch_vv(nn2, n1) = exch_vv(nn2, n1) + cdum*phase_vv(iband, nn2)* &
     341           0 :                                            dotprod(carr1_v_c(:n), cprod_vv_c(:n, iband, nn2))
     342             :                      END DO !n2
     343             :                   END IF
     344             :                END DO
     345             :             END DO  !n1
     346           0 :             call timestop("exchange matrix")
     347             :          END DO !ibando
     348             :       END DO  !ikpt
     349             : 
     350             : !   WRITE(7001,'(a,i7)') 'nk: ', nk
     351             : !   DO n1=1,hybrid%nbands(nk)
     352             : !      DO n2=1,n1
     353             : !         WRITE(7001,'(2i7,2f15.8)') n2, n1, exch_vv(n2,n1)
     354             : !     END DO
     355             : !   END DO
     356             : 
     357             :       ! add contribution of the gamma point to the different cases (exch_vv,exch_cv,exch_cc)
     358             : 
     359             :       ! valence-valence-valence-valence exchange
     360             : 
     361           0 :       IF ((.not. xcpot%is_name("hse")) .AND. (.not. xcpot%is_name("vhse"))) THEN ! no gamma point correction needed for HSE functional
     362             :          IF (zero_order .and. .not. ibs_corr) THEN
     363             :             WRITE (6, '(A)') ' Take zero order terms into account.'
     364             :          ELSE IF (zero_order .and. ibs_corr) THEN
     365             :             WRITE (6, '(A)') ' Take zero order terms and ibs-correction into account.'
     366             :          END IF
     367             : 
     368             :          IF (zero_order) THEN
     369             :             CALL dwavefproducts(dcprod, nk, 1, hybrid%nbands(nk), 1, hybrid%nbands(nk), .false., atoms, hybrid, &
     370             :                                 cell, hybdat, kpts, kpts%nkpt, lapw, dimension, jsp, eig_irr)
     371             : 
     372             :             ! make dcprod hermitian
     373             :             DO n1 = 1, hybrid%nbands(nk)
     374             :                DO n2 = 1, n1
     375             :                   dcprod(n1, n2, :) = (dcprod(n1, n2, :) - conjg(dcprod(n2, n1, :)))/2
     376             :                   dcprod(n2, n1, :) = -conjg(dcprod(n1, n2, :))
     377             :                END DO
     378             :             END DO
     379             : 
     380             :             IF (ibs_corr) THEN
     381             :                CALL ibs_correction(nk, atoms, dimension, input, jsp, hybdat, hybrid, lapw, kpts, kpts%nkpt, cell, mnobd, &
     382             :                                    sym, proj_ibsc, olap_ibsc)
     383             :             END IF
     384             :          END IF
     385             : 
     386             :          !This should be done with w_iks I guess!TODO
     387           0 :          occup = .false.
     388           0 :          DO i = 1, hybrid%ne_eig(nk)
     389           0 :             IF (results%ef >= eig_irr(i, nk)) THEN
     390           0 :                occup(i) = .true.
     391           0 :             ELSE IF ((eig_irr(i, nk) - results%ef) <= 1E-06) THEN
     392           0 :                occup(i) = .true.
     393             :             END IF
     394             :          END DO
     395             : 
     396           0 :          DO n1 = 1, hybrid%nbands(nk)
     397           0 :             DO n2 = 1, nsest(n1)!n1
     398           0 :                nn2 = indx_sest(n2, n1)
     399           0 :                exchcorrect = 0
     400           0 :                exch0 = 0
     401             : 
     402             :                ! if zero_order = .true. add averaged k-dependent term to the numerical integration at Gamma-point contribution
     403             : 
     404             :                ! if we start with a system with a small DFT band gap (like GaAs), the contribution
     405             :                ! of the highest occupied and lowest unoccupied state in Hessian is typically
     406             :                ! large; a correct numerical integration requires a dense k-point mesh, so
     407             :                ! we don't add the contribution exchcorrect for such materials
     408             : 
     409             :                IF (zero_order) THEN
     410             :                   hessian = 0
     411             :                   IF (occup(n1) .and. occup(nn2)) THEN
     412             :                      DO i = 1, 3
     413             :                         j = i
     414             :                         DO iband = 1, hybrid%nbands(nk)
     415             :                            IF (occup(iband)) THEN
     416             :                               hessian(i, j) = hessian(i, j) + conjg(dcprod(iband, n1, i))*dcprod(iband, nn2, j)
     417             :                            END IF
     418             :                            hessian(i, j) = hessian(i, j) - dcprod(iband, nn2, i)*conjg(dcprod(iband, n1, j))
     419             :                         END DO
     420             : 
     421             :                         ! ibs correction
     422             :                         IF (ibs_corr) THEN
     423             :                            hessian(i, j) = hessian(i, j) - olap_ibsc(i, j, n1, nn2)/cell%omtil
     424             :                            DO iband = 1, hybrid%nbands(nk)
     425             :                               hessian(i, j) = hessian(i, j) + conjg(proj_ibsc(i, nn2, iband))*proj_ibsc(j, n1, iband)/cell%omtil
     426             :                            END DO
     427             :                         END IF
     428             :                      END DO
     429             :                   ELSE
     430             :                      DO i = 1, 3
     431             :                         j = i
     432             :                         DO iband = 1, hybrid%nbands(nk)
     433             :                            IF (occup(iband)) THEN
     434             :                               hessian(i, j) = hessian(i, j) + conjg(dcprod(iband, n1, i))*dcprod(iband, nn2, j)
     435             :                            END IF
     436             :                         END DO
     437             :                      END DO
     438             :                   END IF
     439             : 
     440             :                   exchcorrect(1) = fpi_const/3*(hessian(1, 1) + hessian(2, 2) + hessian(3, 3))
     441             :                   exch0 = exchcorrect(1)/kpts%nkptf
     442             :                END IF
     443             : 
     444             :                ! tail correction/contribution from all other k-points (it  goes into exchcorrect )
     445             : 
     446             :                ! Analytic contribution
     447             : 
     448           0 :                cdum2 = 0
     449             :                !multiply divergent contribution with occupation number;
     450             :                !this only affects metals
     451           0 :                IF (n1 == nn2) THEN
     452           0 :                   cdum2 = fpi_const/cell%omtil*divergence*wl_iks(n1, nk)*kpts%nkptf
     453             :                END IF
     454             : 
     455             :                ! due to the symmetrization afterwards the factor 1/n_q(1) must be added
     456             : 
     457           0 :                IF (n1 == nn2) hybrid%div_vv(n1, nk, jsp) = REAL(cdum2)
     458           0 :                exch_vv(nn2, n1) = exch_vv(nn2, n1) + (exch0 + cdum2)/n_q(1)
     459             : 
     460             :             END DO !n2
     461             :          END DO !n1
     462             :       END IF ! xcpot%icorr .ne. icorr_hse
     463             : 
     464           0 :       IF (mat_ex%l_real) THEN
     465           0 :          IF (any(abs(aimag(exch_vv)) > 1E-08)) CALL judft_warn('unusally large imaginary part of exch_vv', &
     466           0 :                                                                calledby='exchange_val_hf.F90')
     467             :       END IF
     468             : 
     469             : !   WRITE(7000,'(a,i7)') 'nk: ', nk
     470             : !   DO n1=1,hybrid%nbands(nk)
     471             : !      DO n2=1,n1
     472             : !         WRITE(7000,'(2i7,2f15.8)') n2, n1, exch_vv(n2,n1)
     473             : !      END DO
     474             : !   END DO
     475             : 
     476             :       ! write exch_vv in mat_ex
     477           0 :       CALL mat_ex%alloc(matsize1=hybrid%nbands(nk))
     478           0 :       IF (mat_ex%l_real) THEN
     479           0 :          mat_ex%data_r = exch_vv
     480             :       ELSE
     481           0 :          mat_ex%data_c = exch_vv
     482             :       END IF
     483           0 :       CALL timestop("valence exchange calculation")
     484             : 
     485           0 :    END SUBROUTINE exchange_valence_hf
     486             : 
     487           0 :    SUBROUTINE calc_divergence(cell, kpts, divergence)
     488             : 
     489             :       USE m_util, ONLY: cerf
     490             :       USE m_types
     491             :       USE m_constants
     492             : 
     493             :       IMPLICIT NONE
     494             : 
     495             :       TYPE(t_cell), INTENT(IN)  :: cell
     496             :       TYPE(t_kpts), INTENT(IN)  :: kpts
     497             :       REAL, INTENT(OUT) :: divergence
     498             : 
     499             :       INTEGER :: ix, iy, iz, sign, n
     500             :       logical :: found
     501             :       REAL    :: expo, rrad, k(3), kv1(3), kv2(3), kv3(3), knorm2
     502             :       COMPLEX :: cdum
     503             : 
     504           0 :       expo = 5e-3
     505           0 :       rrad = sqrt(-log(5e-3)/expo)
     506           0 :       cdum = sqrt(expo)*rrad
     507           0 :       divergence = cell%omtil/(tpi_const**2)*sqrt(pi_const/expo)*cerf(cdum)
     508           0 :       rrad = rrad**2
     509           0 :       kv1 = cell%bmat(1, :)/kpts%nkpt3(1)
     510           0 :       kv2 = cell%bmat(2, :)/kpts%nkpt3(2)
     511           0 :       kv3 = cell%bmat(3, :)/kpts%nkpt3(3)
     512             :       n = 1
     513             :       found = .true.
     514             : 
     515           0 :       DO WHILE (found)
     516           0 :          found = .false.
     517           0 :          DO ix = -n, n
     518           0 :             DO iy = -(n - abs(ix)), n - abs(ix)
     519           0 :                iz = n - abs(ix) - abs(iy)
     520           0 :                DO sign = -1, 1, 2
     521           0 :                   iz = sign*iz
     522           0 :                   k(1) = ix*kv1(1) + iy*kv2(1) + iz*kv3(1)
     523           0 :                   k(2) = ix*kv1(2) + iy*kv2(2) + iz*kv3(2)
     524           0 :                   k(3) = ix*kv1(3) + iy*kv2(3) + iz*kv3(3)
     525           0 :                   knorm2 = k(1)**2 + k(2)**2 + k(3)**2
     526           0 :                   IF (knorm2 < rrad) THEN
     527           0 :                      found = .true.
     528           0 :                      divergence = divergence - exp(-expo*knorm2)/knorm2/kpts%nkptf
     529             :                   END IF
     530           0 :                   IF (iz == 0) exit
     531             :                END DO
     532             :             END DO
     533             :          END DO
     534           0 :          n = n + 1
     535             :       END DO
     536             : 
     537           0 :    END SUBROUTINE calc_divergence
     538             : 
     539             : END MODULE m_exchange_valence_hf

Generated by: LCOV version 1.13