exchange_val_hf.F90 24.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
!     Calculates the HF exchange term 
!
!                                          s          s*          s            s*
!                                       phi    (r) phi     (r) phi     (r') phi    (r')
!                         occ.             n_1k       n'k+q       n'k+q        n_2k
!     exchange(n,q)  =  - SUM  INT INT  ------------------------------------------- dr dr'
!                         k,n'                           | r - r' |
!
!                         occ                  s          s    ~        ~       s         s
!                    =  - SUM  SUM  v     < phi      | phi     M    > < M    phi     | phi      >
!                         k,n' I,J   k,IJ      n'k+q      n_1k  q,I      q,J    n_2k      n'k+q
!
!     for the different combinations of n_1 and n_2 and where n' runs only over the valence states.     
!     ( n_1,n_2:  valence-valence, core-core,core-valence )
!
!
!     At the Gamma point (k=0) v diverges. After diagonalization of v at k=0 the divergence is
!     restricted to the head element I=1. Furthermore, we expand <...> with kp perturbation theory.
!     As a result, the total I=1 element is given by a sum of a divergent 1/k**2-term and an
!     angular dependent term. The former is separated from the numerical k-summation and treated
!     analytically while the latter is spherically averaged and added to the k=0 contribution of
!     the numerical k-summation. (A better knowledge of the integrand's behavior at the BZ edges
!     might further improve the integration.)
!
!     The divergence at the Gamma point is integrated with one of the following algorithms:
! (1) Switching-Off Function
!     In a sphere of radius k0=radshmin/2 a switching-off function g(k)=1-(k/k0)**n*(n+1-n*k/k0)
!     (n=npot) is defined. The 1/k**2 divergence is subtracted from the BZ integral in the form
!     g(k)/k**2 and integrated analytically. The non-divergent rest is integrated numerically.
! (2) Periodic Function (similar to the one used by Massidda PRB 48, 5058)
!     The function  F(k) = SUM(G) exp(-expo*|k+G|**3) / |k+G|**2  is subtracted from the BZ integral
!     and integrated analytically. The non-divergent rest is integrated numerically.
!     The parameter expo is chosen such that exp(-expo*q**3)=1/2
!     with q = radius of sphere with same volume as BZ.
! (3) Periodic Function (same as Massidda's) with expo->0
!     The function  F(k) = lim(expo->0) SUM(G) exp(-expo*|k+G|**2) / |k+G|**2  is subtracted from
!     the BZ integral and integrated analytically. The contribution to the BZ integral including
!     the "tail" is
!     vol/(8*pi**3) INT F(k) d^3k - P SUM(k) F(k)  ( P = principal value ) .
!     For expo->0 the two terms diverge. Therefore a cutoff radius q0 is introduced and related to
!     expo by exp(-expo*q0**2)=delta  ( delta = small value, e.g., delta = 1d-10 ) .
!     The resulting formula
!     vol/(4*pi**1.5*sqrt(expo)) * erf(sqrt(a)*q0) - sum(q,0<q<q0) exp(-expo*q**2)/q**2
!     converges well with q0. (Should be the default.)
      MODULE m_exchange_valence_hf

Daniel Wortmann's avatar
Daniel Wortmann committed
47 48
        LOGICAL,PARAMETER:: zero_order=.false.,ibs_corr=.false.
        INTEGER,PARAMETER:: maxmem=600
49 50 51 52

      CONTAINS

      SUBROUTINE exchange_valence_hf(&
Daniel Wortmann's avatar
Daniel Wortmann committed
53
                    nk,kpts,nkpt_EIBZ, sym,atoms,hybrid,&
Daniel Wortmann's avatar
Daniel Wortmann committed
54
                    cell, dimension,input,jsp, hybdat, mnobd, lapw,&
55
                    eig_irr,results,parent,pointer_EIBZ,n_q,wl_iks,&
Daniel Wortmann's avatar
Daniel Wortmann committed
56
                    it,xcpot, noco,nsest,indx_sest,&
57
                    mpi,irank2,isize2,comm,mat_ex)
58 59 60 61 62 63 64 65 66


      USE m_wrapper
      USE m_constants   
      USE m_trafo
      USE m_wavefproducts
      USE m_olap
      USE m_spmvec
      USE m_hsefunctional ,ONLY: dynamic_hse_adjustment
67
#if defined(CPP_MPI)&&defined(CPP_NEVER)
68 69 70
      USE m_mpi_work_dist
      USE m_mpi_tags
#endif
Daniel Wortmann's avatar
Daniel Wortmann committed
71
      USE m_io_hybrid
72 73 74
      USE m_kp_perturbation
      USE m_types
      IMPLICIT NONE
Daniel Wortmann's avatar
Daniel Wortmann committed
75
      TYPE(t_hybdat),INTENT(IN)   :: hybdat
76 77 78 79
      TYPE(t_results),INTENT(IN)   :: results
      TYPE(t_xcpot),INTENT(IN)   :: xcpot
      TYPE(t_mpi),INTENT(IN)   :: mpi
      TYPE(t_dimension),INTENT(IN)   :: dimension
Daniel Wortmann's avatar
Daniel Wortmann committed
80
      TYPE(t_hybrid),INTENT(INOUT)   :: hybrid
81 82 83 84 85 86 87 88 89 90 91
      TYPE(t_input),INTENT(IN)   :: input
      TYPE(t_noco),INTENT(IN)   :: noco
      TYPE(t_sym),INTENT(IN)   :: sym
      TYPE(t_cell),INTENT(IN)   :: cell
      TYPE(t_kpts),INTENT(IN)   :: kpts
      TYPE(t_atoms),INTENT(IN)   :: atoms
      TYPE(t_lapw),INTENT(IN)   :: lapw

!     - scalars -
      INTEGER,INTENT(IN)      :: it  ,irank2 ,isize2,comm
      INTEGER,INTENT(IN)      :: jsp
Daniel Wortmann's avatar
Daniel Wortmann committed
92
      INTEGER,INTENT(IN)      ::  nk ,nkpt_EIBZ
Daniel Wortmann's avatar
Daniel Wortmann committed
93
      INTEGER,INTENT(IN)      :: mnobd 
94 95 96 97 98



!     - arrays -
      INTEGER,INTENT(IN)      ::  n_q(nkpt_EIBZ)
Daniel Wortmann's avatar
Daniel Wortmann committed
99

100 101
      INTEGER,INTENT(IN)      ::  parent(kpts%nkptf)
      INTEGER,INTENT(IN)      ::  pointer_EIBZ(nkpt_EIBZ)
Daniel Wortmann's avatar
Daniel Wortmann committed
102
      INTEGER,INTENT(IN)      ::  nsest(hybrid%nbands(nk)),indx_sest(hybrid%nbands(nk),hybrid%nbands(nk))
Daniel Wortmann's avatar
Daniel Wortmann committed
103 104

 
Daniel Wortmann's avatar
Daniel Wortmann committed
105
      REAL   ,INTENT(IN)      ::  eig_irr(dimension%neigd,kpts%nkpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
106

107
      REAL   ,INTENT(IN)      ::  wl_iks(dimension%neigd,kpts%nkptf)
108
    
Daniel Wortmann's avatar
Daniel Wortmann committed
109
      TYPE(t_mat),INTENT(INOUT):: mat_ex
110 111 112 113 114 115

!     - local scalars -
      INTEGER                 ::  iband,iband1,ibando,ikpt,ikpt0
      INTEGER                 ::  i,ic,ix,iy,iz
      INTEGER                 ::  irecl_coulomb,irecl_coulomb1
      INTEGER                 ::  j
Daniel Wortmann's avatar
Daniel Wortmann committed
116
      INTEGER                 ::  m1,m2
117 118 119 120 121
      INTEGER                 ::  n,n1,n2,nn,nn2
      INTEGER                 ::  nkqpt
      INTEGER                 ::  npot
      INTEGER                 ::  ok
      INTEGER                 ::  psize
Daniel Wortmann's avatar
Daniel Wortmann committed
122 123 124
      REAL                    ::  rdum
      REAL                    ::  k0
     
125 126 127 128
      REAL , SAVE             ::  divergence

      COMPLEX                 ::  cdum,cdum1,cdum2 
      COMPLEX                 ::  exch0
Daniel Wortmann's avatar
Daniel Wortmann committed
129
 
130 131 132 133 134
      LOGICAL, SAVE           ::  initialize = .true.

!     - local arrays -
      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/) )
Daniel Wortmann's avatar
Daniel Wortmann committed
135
   
136 137
      COMPLEX,ALLOCATABLE     ::  phase_vv(:,:)
      COMPLEX                 ::  exchcorrect(kpts%nkptf)
Daniel Wortmann's avatar
Daniel Wortmann committed
138
      COMPLEX                 ::  dcprod(hybrid%nbands(nk),hybrid%nbands(nk),3) 
139

Daniel Wortmann's avatar
Daniel Wortmann committed
140
      COMPLEX(8)              ::  exch_vv(hybrid%nbands(nk),hybrid%nbands(nk))
141
#if defined(CPP_MPI)&&defined(CPP_NEVER)
Daniel Wortmann's avatar
Daniel Wortmann committed
142
      COMPLEX(8)              ::  buf_vv(hybrid%nbands(nk),nbands(nk))
143 144
#endif
      COMPLEX                 ::  hessian(3,3)
Daniel Wortmann's avatar
Daniel Wortmann committed
145
      COMPLEX                 ::  proj_ibsc(3,mnobd,hybrid%nbands(nk))
146 147
      COMPLEX                 ::  olap_ibsc(3,3,mnobd,mnobd)
#if ( !defined CPP_NOSPMVEC && !defined CPP_IRAPPROX )
Daniel Wortmann's avatar
Daniel Wortmann committed
148
      REAL                    ::  coulomb_mt1(hybrid%maxindxm1-1,hybrid%maxindxm1-1, 0:hybrid%maxlcutm1,atoms%ntype)       
Daniel Wortmann's avatar
Daniel Wortmann committed
149 150 151 152
      REAL                    ::  coulomb_mt2_r(hybrid%maxindxm1-1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1+1,atoms%nat)
      REAL                    ::  coulomb_mt3_r(hybrid%maxindxm1-1,atoms%nat,atoms%nat)
      COMPLEX                 ::  coulomb_mt2_c(hybrid%maxindxm1-1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1+1,atoms%nat)
      COMPLEX                 ::  coulomb_mt3_c(hybrid%maxindxm1-1,atoms%nat,atoms%nat)
153 154 155

#else

Daniel Wortmann's avatar
Daniel Wortmann committed
156 157
      REAL                    ::  coulomb_r(hybrid%maxbasm1*(hybrid%maxbasm1+1)/2)
      COMPLEX                 ::  coulomb_c(hybrid%maxbasm1*(hybrid%maxbasm1+1)/2)
158 159 160

#endif

Daniel Wortmann's avatar
Daniel Wortmann committed
161 162
      REAL   ,ALLOCATABLE     ::  cprod_vv_r(:,:,:),cprod_cv_r(:,:,:), carr3_vv_r(:,:,:),carr3_cv_r(:,:,:)
      REAL                    ::  carr1_v_r(hybrid%maxbasm1),carr1_c_r(hybrid%maxbasm1)
163
#ifdef CPP_IRCOULOMBAPPROX
Daniel Wortmann's avatar
Daniel Wortmann committed
164
      REAL                    ::  coulomb_mtir_r((hybrid%maxlcutm1+1)**2*atoms%nat , (hybrid%maxlcutm1+1)**2*atoms%nat +maxval(hybrid%ngptm) )
165
#else
Daniel Wortmann's avatar
Daniel Wortmann committed
166
      REAL                    ::  coulomb_mtir_r(((hybrid%maxlcutm1+1)**2*atoms%nat +maxval(hybrid%ngptm))* ((hybrid%maxlcutm1+1)**2*atoms%nat +maxval(hybrid%ngptm)+1)/2 )
167
#endif
Daniel Wortmann's avatar
Daniel Wortmann committed
168 169
      COMPLEX,ALLOCATABLE     ::  cprod_vv_c(:,:,:),cprod_cv_c(:,:,:), carr3_vv_c(:,:,:),carr3_cv_c(:,:,:)
      COMPLEX                 ::  carr1_v_c(hybrid%maxbasm1),carr1_c_c(hybrid%maxbasm1)
170 171

#ifdef CPP_IRCOULOMBAPPROX
Daniel Wortmann's avatar
Daniel Wortmann committed
172
      COMPLEX                 ::  coulomb_mtir_c((hybrid%maxlcutm1+1)**2*atoms%nat , (hybrid%maxlcutm1+1)**2*atoms%nat +maxval(hybrid%ngptm) )
173
#else
Daniel Wortmann's avatar
Daniel Wortmann committed
174
      COMPLEX                 ::  coulomb_mtir_c(((hybrid%maxlcutm1+1)**2*atoms%nat +maxval(hybrid%ngptm))* ((hybrid%maxlcutm1+1)**2*atoms%nat +maxval(hybrid%ngptm)+1)/2 )
175 176 177
#endif

      LOGICAL                 ::  occup(dimension%neigd)
178
#if defined(CPP_MPI)&&defined(CPP_NEVER)
179 180 181 182
      INCLUDE "mpif.h"
      INTEGER                 :: ierr,ierr2,length,rank
      CHARACTER(LEN=MPI_MAX_ERROR_STRING) :: errmsg
#endif
Daniel Wortmann's avatar
Daniel Wortmann committed
183 184
     
     
185
      IF( initialize ) THEN !it .eq. 1 .and. nk .eq. 1) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
186
         call calc_divergence(cell,kpts,divergence)
187
         PRINT *,"Divergence:",divergence
Daniel Wortmann's avatar
Daniel Wortmann committed
188
         initialize = .false.
189
      END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
190
   
191 192 193 194 195 196 197 198

      ! calculate valence-valence-valence-valence, core-valence-valence-valence
      ! and core-valence-valence-core exchange at current k-point
      ! the sum over the inner occupied valence states is restricted to the EIBZ(k)
      ! the contribution of the Gamma-point is treated separately (see below)


      ! determine package size loop over the occupied bands
Daniel Wortmann's avatar
Daniel Wortmann committed
199
      if (mat_ex%l_real) THEn
Daniel Wortmann's avatar
Daniel Wortmann committed
200
         rdum  = hybrid%maxbasm1*hybrid%nbands(nk)*4/1048576.
Daniel Wortmann's avatar
Daniel Wortmann committed
201
      else
Daniel Wortmann's avatar
Daniel Wortmann committed
202
         rdum  = hybrid%maxbasm1*hybrid%nbands(nk)*4/1048576.
Daniel Wortmann's avatar
Daniel Wortmann committed
203
      endif
204 205 206 207 208 209 210 211 212 213 214 215 216
      psize = 1
      DO iband = mnobd,1,-1
        ! ensure that the packages have equal size
        IF( modulo(mnobd,iband) .eq. 0 ) THEN
          ! choose packet size such that cprod is smaller than memory threshold
          IF( rdum*iband .le. maxmem ) THEN
            psize = iband
            EXIT
          END IF
        END IF
      END DO

      IF( psize .ne. mnobd ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
217
        WRITE(6,'(A,A,i3,A,f7.2,A)') ' Divide the loop over the occupied hybrid%bands in packages', ' of the size',psize,' (cprod=',rdum*psize,'MB)'
218
      END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
219
      ALLOCATE( phase_vv(psize,hybrid%nbands(nk)),stat=ok )
220
      IF( ok .ne. 0 ) STOP 'exchange_val_hf: error allocation phase'
Daniel Wortmann's avatar
Daniel Wortmann committed
221 222 223
      phase_vv=0
      IF( ok .ne. 0 ) STOP 'exchange_val_hf: error allocation phase'
      if (mat_ex%l_real) THEN
224
         ALLOCATE( cprod_vv_c(0,0,0), carr3_vv_c(0,0,0))
Daniel Wortmann's avatar
Daniel Wortmann committed
225
         ALLOCATE( cprod_vv_r(hybrid%maxbasm1,psize,hybrid%nbands(nk)),stat=ok )
Daniel Wortmann's avatar
Daniel Wortmann committed
226
         IF( ok .ne. 0 ) STOP 'exchange_val_hf: error allocation cprod'
Daniel Wortmann's avatar
Daniel Wortmann committed
227
         ALLOCATE( carr3_vv_r(hybrid%maxbasm1,psize,hybrid%nbands(nk)),stat=ok )
Daniel Wortmann's avatar
Daniel Wortmann committed
228 229 230
         IF( ok .ne. 0 ) STOP 'exchange_val_hf: error allocation carr3'
         cprod_vv_r = 0 ; carr3_vv_r = 0 
      ELSE
231
         ALLOCATE( cprod_vv_r(0,0,0), carr3_vv_r(0,0,0))
Daniel Wortmann's avatar
Daniel Wortmann committed
232
         ALLOCATE( cprod_vv_c(hybrid%maxbasm1,psize,hybrid%nbands(nk)),stat=ok )
Daniel Wortmann's avatar
Daniel Wortmann committed
233
         IF( ok .ne. 0 ) STOP 'exchange_val_hf: error allocation cprod'
Daniel Wortmann's avatar
Daniel Wortmann committed
234
         ALLOCATE( carr3_vv_c(hybrid%maxbasm1,psize,hybrid%nbands(nk)),stat=ok )
Daniel Wortmann's avatar
Daniel Wortmann committed
235 236 237 238
         IF( ok .ne. 0 ) STOP 'exchange_val_hf: error allocation carr3'
         cprod_vv_c = 0 ; carr3_vv_c = 0
      endif
         
239 240
      exch_vv = 0

Daniel Wortmann's avatar
Daniel Wortmann committed
241
      DO ikpt = 1,nkpt_EIBZ
242 243
        ikpt0 = pointer_EIBZ(ikpt)

Daniel Wortmann's avatar
Daniel Wortmann committed
244
        n  = hybrid%nbasp + hybrid%ngptm(ikpt0)
Daniel Wortmann's avatar
Daniel Wortmann committed
245
        IF( hybrid%nbasm(ikpt0) .ne. n ) STOP 'error hybrid%nbasm'
246 247 248 249
        nn = n*(n+1)/2

        ! read in coulomb matrix from direct access file coulomb
#if( !defined CPP_NOSPMVEC && !defined CPP_IRAPPROX )
250 251 252 253 254
        IF (mat_ex%l_real) THEN
	   CALL read_coulomb_spm_r(kpts%bkp(ikpt0),coulomb_mt1,coulomb_mt2_r,coulomb_mt3_r,coulomb_mtir_r)
        ELSE
           CALL read_coulomb_spm_c(kpts%bkp(ikpt0),coulomb_mt1,coulomb_mt2_c,coulomb_mt3_c,coulomb_mtir_c)
        END IF
255
#else
Daniel Wortmann's avatar
Daniel Wortmann committed
256
	call read_coulomb(kpts%bkp(ikpt0),coulomb)
257 258 259 260
#endif

        IF( kpts%bkp(ikpt0) .ne. ikpt0 ) THEN
#if( !defined CPP_NOSPMVEC && !defined CPP_IRAPPROX )
Daniel Wortmann's avatar
Daniel Wortmann committed
261 262 263
          IF( kpts%bksym(ikpt0) .gt. sym%nop.and..not.mat_ex%l_real ) THEN
            coulomb_mt2_c = conjg(coulomb_mt2_c)
            coulomb_mtir_c= conjg(coulomb_mtir_c)
264 265 266 267
          END IF

#else

Daniel Wortmann's avatar
Daniel Wortmann committed
268
          if (.not.mat_ex%l_real) THEN
269
          IF( kpts%bksym(ikpt0) .gt. sym%nop ) coulomb = conjg(coulomb)
Daniel Wortmann's avatar
Daniel Wortmann committed
270
       endif
271 272 273 274 275

#endif
        END IF

        DO ibando = 1,mnobd,psize
Daniel Wortmann's avatar
Daniel Wortmann committed
276
        if (mat_ex%l_real) THEN
277 278 279

#ifdef CPP_IRAPPROX
          CALL wavefproducts_inv(&
Daniel Wortmann's avatar
Daniel Wortmann committed
280
                         1,hybdat,dimension,jsp,atoms,&
Daniel Wortmann's avatar
Daniel Wortmann committed
281
                         lapw,obsolete,kpts,&
Daniel Wortmann's avatar
Daniel Wortmann committed
282
                         nk,ikpt0,mnobd,hybrid, parent,cell, sym,&
Daniel Wortmann's avatar
Daniel Wortmann committed
283
                         nkqpt,cprod_vv)
284 285
#else
          CALL wavefproducts_inv5(&
Daniel Wortmann's avatar
Daniel Wortmann committed
286
                         1,hybrid%nbands(nk),ibando,ibando+psize-1,&
287
                         dimension,input,jsp,atoms,&
Daniel Wortmann's avatar
Daniel Wortmann committed
288 289 290 291 292
                         lapw,kpts,&
                         nk,ikpt0,hybdat,mnobd,hybrid,&
                         parent,cell,hybrid%nbasp,sym,&
                         noco,&
                         nkqpt,cprod_vv_r)
293 294
#endif

Daniel Wortmann's avatar
Daniel Wortmann committed
295
       else
296 297
#ifdef CPP_IRAPPROX
          CALL wavefproducts_noinv(&
Daniel Wortmann's avatar
Daniel Wortmann committed
298 299 300 301
                         1,hybdat,nk,ikpt0,dimension,jsp,&
                         cell,atoms,hybrid, 
                         kpts,&
                         mnobd,&
302
                         lapw,sym,&
Daniel Wortmann's avatar
Daniel Wortmann committed
303
                         nkqpt,&
304 305 306
                         cprod_vv)
#else
          CALL wavefproducts_noinv5(&
Daniel Wortmann's avatar
Daniel Wortmann committed
307
                         1,hybrid%nbands(nk),ibando,ibando+psize-1,&
308
                         nk,ikpt0,dimension,input,jsp, &!jsp,&
Daniel Wortmann's avatar
Daniel Wortmann committed
309 310 311 312
                         cell,atoms,hybrid,hybdat, &
                         kpts,&
                         mnobd,&
                         lapw,sym, &
Daniel Wortmann's avatar
Daniel Wortmann committed
313 314
                         hybrid%nbasp,noco,&
                         nkqpt,cprod_vv_c)
315 316
#endif

Daniel Wortmann's avatar
Daniel Wortmann committed
317
endif
318 319 320 321 322 323 324

          ! The sparse matrix technique is not feasible for the HSE
          ! functional. Thus, a dynamic adjustment is implemented
          ! The mixed basis functions and the potential difference
          ! are Fourier transformed, so that the exchange can be calculated
          ! in Fourier space
#ifndef CPP_NOSPMVEC
325
          IF ( xcpot%is_name("hse") .OR. xcpot%is_name("vhse") ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
326
            iband1  = hybrid%nobd(nkqpt)
327
            exch_vv = exch_vv + dynamic_hse_adjustment(&
Daniel Wortmann's avatar
Daniel Wortmann committed
328
                       atoms%rmsh,atoms%rmt,atoms%dx,atoms%jri,atoms%jmtd,kpts%bkf(:,ikpt0),ikpt0,kpts%nkptf,&
Daniel Wortmann's avatar
Daniel Wortmann committed
329
                       cell%bmat,cell%omtil,atoms%ntype,atoms%neq,atoms%nat,atoms%taual,hybrid%lcutm1,hybrid%maxlcutm1,&
Daniel Wortmann's avatar
Daniel Wortmann committed
330
                       hybrid%nindxm1,hybrid%maxindxm1,hybrid%gptm,hybrid%ngptm(ikpt0),hybrid%pgptm(:,ikpt0),&
Daniel Wortmann's avatar
Daniel Wortmann committed
331
                       hybrid%gptmd,hybrid%basm1,hybrid%nbasm(ikpt0),iband1,hybrid%nbands(nk),nsest,&
332
                       ibando,psize,indx_sest,atoms%invsat,sym%invsatnr,mpi%irank,&
Daniel Wortmann's avatar
Daniel Wortmann committed
333 334
                       cprod_vv_r(:hybrid%nbasm(ikpt0),:,:),&
                       cprod_vv_c(:hybrid%nbasm(ikpt0),:,:),&
Daniel Wortmann's avatar
Daniel Wortmann committed
335
                       mat_ex%l_real,wl_iks(:iband1,nkqpt),n_q(ikpt))
336 337 338 339 340 341 342
          END IF
#endif

          ! the Coulomb matrix is only evaluated at the irrecuible k-points
          ! bra_trafo transforms cprod instead of rotating the Coulomb matrix
          ! from IBZ to current k-point
          IF( kpts%bkp(ikpt0) .ne. ikpt0 ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
343
             CALL bra_trafo2(&
344 345 346 347 348 349
                  mat_ex%l_real,carr3_vv_r(:hybrid%nbasm(ikpt0),:,:),cprod_vv_r(:hybrid%nbasm(ikpt0),:,:),&
                  carr3_vv_c(:hybrid%nbasm(ikpt0),:,:),cprod_vv_c(:hybrid%nbasm(ikpt0),:,:),&
                  hybrid%nbasm(ikpt0),psize,hybrid%nbands(nk),&
                  kpts%bkp(ikpt0),ikpt0,kpts%bksym(ikpt0),sym,&
                  hybrid,kpts,cell,atoms,&
                  phase_vv)
Daniel Wortmann's avatar
Daniel Wortmann committed
350 351

             IF (mat_ex%l_real) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
352
                cprod_vv_r(:hybrid%nbasm(ikpt0),:,:) = carr3_vv_r(:hybrid%nbasm(ikpt0),:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
353
             ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
354
                cprod_vv_c(:hybrid%nbasm(ikpt0),:,:) = carr3_vv_c(:hybrid%nbasm(ikpt0),:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
355
             ENDIF
356 357 358 359 360 361
          ELSE
            phase_vv(:,:) = (1d0,0d0)
          END IF

          ! calculate exchange matrix at ikpt0

Daniel Wortmann's avatar
Daniel Wortmann committed
362
          DO n1=1,hybrid%nbands(nk)
363
            DO iband = 1,psize
Daniel Wortmann's avatar
Daniel Wortmann committed
364
              IF( ibando + iband - 1 .gt. hybrid%nobd(nkqpt) ) CYCLE
365 366 367 368
              cdum  = wl_iks(ibando+iband-1,nkqpt)&
                    * conjg(phase_vv(iband,n1))/n_q(ikpt)

#if( !defined CPP_NOSPMVEC && !defined CPP_IRAPPROX )
Daniel Wortmann's avatar
Daniel Wortmann committed
369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385
              if (mat_ex%l_real) THEN
                 carr1_v_r(:n) = 0 
                 CALL spmvec_invs(atoms,hybrid,&
                      hybdat,ikpt0,kpts,&
                      cell,&
                      coulomb_mt1,coulomb_mt2_r,coulomb_mt3_r,&
                      coulomb_mtir_r,cprod_vv_r(:n,iband,n1),&
                      carr1_v_r(:n))
              ELSE
                 carr1_v_c(:n) = 0 
                 CALL spmvec_noinvs(atoms,hybrid,&
                      hybdat,ikpt0,kpts,&
                      cell,&
                      coulomb_mt1,coulomb_mt2_c,coulomb_mt3_c,&
                      coulomb_mtir_c,cprod_vv_c(:n,iband,n1),&
                      carr1_v_c(:n))
              endif
386
#else
Daniel Wortmann's avatar
Daniel Wortmann committed
387 388 389 390 391
              if (mat_ex%l_real) THEN
                 carr1_v_r(:n) = matvec( coulomb_r(:nn),cprod_vv_r(:n,iband,n1) )
              ELSE
                 carr1_v_r(:n) = matvec( coulomb_c(:nn),cprod_vv_c(:n,iband,n1) )
              endif
392 393
#endif

Daniel Wortmann's avatar
Daniel Wortmann committed
394 395 396 397 398 399 400 401 402 403 404 405 406
              if (mat_ex%l_real) THEN
                 DO n2=1,nsest(n1)!n1
                    nn2 = indx_sest(n2,n1)
                    exch_vv(nn2,n1) = exch_vv(nn2,n1) + cdum*phase_vv(iband,nn2)&
                         *dotprod( carr1_v_r(:n), cprod_vv_r(:n,iband,nn2) )
                 END DO !n2
              else
                 DO n2=1,nsest(n1)!n1
                    nn2 = indx_sest(n2,n1)
                    exch_vv(nn2,n1) = exch_vv(nn2,n1) + cdum*phase_vv(iband,nn2)&
                         *dotprod( carr1_v_c(:n), cprod_vv_c(:n,iband,nn2) )
                 END DO !n2
              end if
407 408 409 410 411 412 413 414 415
            END DO
          END DO  !n1
        END DO !ibando
      END DO  !ikpt


      !
      ! add contribution of the gamma point to the different cases (exch_vv,exch_cv,exch_cc)
      !
Daniel Wortmann's avatar
Daniel Wortmann committed
416
     
417 418
      ! valence-valence-valence-valence exchange

419
      IF ( (.not.xcpot%is_name("hse")) .AND. (.not.xcpot%is_name("vhse")) ) THEN ! no gamma point correction needed for HSE functional
420 421 422 423 424 425 426
        IF( zero_order .and. .not. ibs_corr ) THEN
          WRITE(6,'(A)') ' Take zero order terms into account.'
        ELSE IF( zero_order .and.  ibs_corr ) THEN
          WRITE(6,'(A)') ' Take zero order terms and ibs-correction into account.'
        END IF
        IF( zero_order ) THEN
          CALL dwavefproducts(  &
Daniel Wortmann's avatar
Daniel Wortmann committed
427
                            dcprod,nk,1,hybrid%nbands(nk),1,hybrid%nbands(nk),.false., atoms,hybrid,&
Daniel Wortmann's avatar
Daniel Wortmann committed
428
                            cell,hybdat, kpts,kpts%nkpt,lapw,&
Daniel Wortmann's avatar
Daniel Wortmann committed
429 430
                            dimension,jsp,&
                            eig_irr )
431 432

          ! make dcprod hermitian
Daniel Wortmann's avatar
Daniel Wortmann committed
433
          DO n1 = 1,hybrid%nbands(nk)
434 435 436 437 438 439 440 441 442 443 444
            DO n2 = 1,n1
              dcprod(n1,n2,:) = (dcprod(n1,n2,:) &
                              - conjg(dcprod(n2,n1,:)))/2   
              dcprod(n2,n1,:) = -conjg(dcprod(n1,n2,:))
            END DO
          END DO

          IF( ibs_corr ) THEN
            CALL ibs_correction(&
                        nk,atoms,&
                        dimension,input,jsp,&
Daniel Wortmann's avatar
Daniel Wortmann committed
445
                        hybdat,hybrid,&
Daniel Wortmann's avatar
Daniel Wortmann committed
446
                        lapw,kpts,kpts%nkpt,&
Daniel Wortmann's avatar
Daniel Wortmann committed
447 448
                        cell,mnobd,&
                        sym,&
449 450 451 452
                        proj_ibsc,olap_ibsc)
          END IF

        END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
453 454
        
        !This should be done with w_iks I guess!TODO
455
        occup = .false.
Daniel Wortmann's avatar
Daniel Wortmann committed
456
        DO i=1,hybrid%ne_eig(nk)
457 458 459 460 461 462 463 464
          IF ( results%ef  .ge. eig_irr(i,nk) ) THEN
            occup(i) = .true.
          ELSE IF ( eig_irr(i,nk) - results%ef .le. 1E-06) THEN
             occup(i) = .true.
          END IF
        END DO


Daniel Wortmann's avatar
Daniel Wortmann committed
465
        DO n1 = 1,hybrid%nbands(nk)
466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484
          DO n2 = 1,nsest(n1)!n1
            nn2 = indx_sest(n2,n1)
            exchcorrect = 0
            exch0       = 0

            ! if zero_order = .true. add averaged k-dependent term to the numerical
            ! integration at Gamma-point contribution
            !
            ! if we start with a system with a small DFT band gap (like GaAs), the contribution
            ! of the highest occupied and lowest unoccupied state in Hessian is typically
            ! large; a correct numerical integration requires a dense k-point mesh, so
            ! we don't add the contribution exchcorrect for such materials 

            IF( zero_order ) THEN
              hessian = 0
              IF( occup(n1) .and. occup(nn2) ) THEN
                DO i = 1,3
                  j = i

Daniel Wortmann's avatar
Daniel Wortmann committed
485
                  DO iband = 1,hybrid%nbands(nk)
486 487 488 489 490 491 492 493
                    IF( occup(iband) ) THEN
                      hessian(i,j) = hessian(i,j) + conjg(dcprod(iband,n1,i)) *dcprod(iband,nn2,j)
                    END IF
                    hessian(i,j) = hessian(i,j) - dcprod(iband,nn2,i) * conjg(dcprod(iband,n1,j))
                  END DO

                  ! ibs correction
                  IF( ibs_corr ) THEN 
Daniel Wortmann's avatar
Daniel Wortmann committed
494
                    hessian(i,j) = hessian(i,j) - olap_ibsc(i,j,n1,nn2)/cell%omtil
Daniel Wortmann's avatar
Daniel Wortmann committed
495
                    DO iband = 1,hybrid%nbands(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
496
                      hessian(i,j) = hessian(i,j) + conjg(proj_ibsc(i,nn2,iband)) * proj_ibsc(j,n1,iband)/cell%omtil
497 498 499 500 501 502 503 504
                    END DO
                  END IF

                END DO
              ELSE

                DO i = 1,3
                  j = i 
Daniel Wortmann's avatar
Daniel Wortmann committed
505
                  DO iband = 1,hybrid%nbands(nk)
506 507 508 509 510 511 512
                    IF( occup(iband) ) THEN
                      hessian(i,j) = hessian(i,j) + conjg(dcprod(iband,n1,i)) * dcprod(iband,nn2,j)
                    END IF
                  END DO
                END DO

              END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
513 514
 
              exchcorrect(1) = fpi_const/3 * (hessian(1,1)+hessian(2,2)+hessian(3,3))
515 516 517 518 519 520 521 522 523 524 525 526
              exch0          = exchcorrect(1)/kpts%nkptf
            END IF


            ! tail correction/contribution from all other k-points (it  goes into exchcorrect )

            ! Analytic contribution

            cdum2 = 0
            !multiply divergent contribution with occupation number;
            !this only affects metals 
            IF ( n1 .eq. nn2 ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
527
               cdum2 = fpi_const/cell%omtil * divergence * wl_iks(n1,nk)*kpts%nkptf
528 529 530 531
            END IF

            ! due to the symmetrization afterwards the factor 1/n_q(1) must be added

532
            IF( n1 .EQ. nn2 ) hybrid%div_vv(n1,nk,jsp) = REAL(cdum2) 
533 534 535 536 537

            exch_vv(nn2,n1)  = exch_vv(nn2,n1) + (exch0 + cdum2)/n_q(1)

          END DO !n2
        END DO !n1
538
       END IF ! xcpot%icorr .ne. icorr_hse
539 540


Daniel Wortmann's avatar
Daniel Wortmann committed
541 542 543
      if (.not.mat_ex%l_real) THEN
         IF(any( abs(aimag(exch_vv)) .gt. 1E-08)) CALL judft_warn('unusally large imaginary part of exch_vv',calledby='exchange_val_hf.F90')
      ENDIF
544 545

      ! write exch_vv in mat_ex
546 547 548 549 550 551
      CALL mat_ex%alloc(matsize1=hybrid%nbands(nk))
      IF (mat_ex%l_real) THEN
         mat_ex%data_r=exch_vv
      ELSE
         mat_ex%data_c=exch_vv
      END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
552
     
553 554
      END SUBROUTINE exchange_valence_hf

Daniel Wortmann's avatar
Daniel Wortmann committed
555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602
      subroutine calc_divergence(cell,kpts,divergence)
        USE m_util          ,ONLY: cerf
        use m_types
        use m_constants
        implicit none
        TYPE(t_cell),INTENT(IN) :: cell
        TYPE(t_kpts),INTENT(IN) :: kpts
        REAL,INTENT(OUT)        :: divergence
        
        INTEGER :: ix,iy,iz,sign,n
        logical :: found
        REAL    :: expo,rrad,k(3),kv1(3),kv2(3),kv3(3),knorm2
        COMPLEX :: cdum
        
        expo       = 5d-3
        rrad       = sqrt(-log(5d-3)/expo)
        cdum       = sqrt(expo)*rrad
        divergence = cell%omtil / (tpi_const**2) * sqrt(pi_const/expo) * cerf(cdum)
        rrad       = rrad**2
        kv1        = cell%bmat(1,:)/kpts%nkpt3(1)
        kv2        = cell%bmat(2,:)/kpts%nkpt3(2)
        kv3        = cell%bmat(3,:)/kpts%nkpt3(3)
        n          = 1
        found      = .true.
        DO WHILE(found)
           found = .false.
           DO ix = -n,n
              DO iy = -(n-abs(ix)),n-abs(ix)
                 iz     = n - abs(ix) - abs(iy)
                 DO sign=-1,1,2
                    iz=sign*iz
                    k(1)   = ix*kv1(1) + iy*kv2(1) + iz*kv3(1)
                    k(2)   = ix*kv1(2) + iy*kv2(2) + iz*kv3(2)
                    k(3)   = ix*kv1(3) + iy*kv2(3) + iz*kv3(3)
                    knorm2 = k(1)**2   + k(2)**2   + k(3)**2
                    IF(knorm2.lt.rrad) THEN
                       found      = .true.
                       divergence = divergence - exp(-expo*knorm2)/knorm2 / kpts%nkptf
                    END IF
                    IF(iz==0) exit
                 enddo
                 
              END DO
           END DO
           n = n + 1
        END DO
      end subroutine calc_divergence
    END MODULE m_exchange_valence_hf