! 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