LCOV - code coverage report
Current view: top level - hybrid - hsefunctional.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 670 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 25 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             : ! Provides all procedures needed to implement the HSE-exchange functional
       8             : ! within a FLAPW framework
       9             : ! Author: M. Schlipf 2009
      10             : MODULE m_hsefunctional
      11             :    USE m_judft
      12             :    IMPLICIT NONE
      13             : 
      14             : #ifdef __PGI
      15             :    REAL, EXTERNAL ::erfc
      16             : #endif
      17             :    ! Constant omega of the HSE exchange functional
      18             :    REAL, PARAMETER :: omega_HSE = 0.11
      19             : 
      20             :    ! Constant for the maximum number of G points
      21             :    INTEGER, PARAMETER :: maxNoGPts = 50
      22             : 
      23             :    ! these arrays are calculated once and then reused
      24             :    LOGICAL, ALLOCATABLE   :: already_known(:)
      25             :    REAL, ALLOCATABLE      :: known_potential(:, :)
      26             : #ifdef CPP_INVERSION
      27             :    REAL, ALLOCATABLE      :: known_fourier_trafo(:, :, :)
      28             : #else
      29             :    COMPLEX, ALLOCATABLE   :: known_fourier_trafo(:, :, :)
      30             : #endif
      31             : 
      32             :    ! variables needed for fourier transformation
      33             :    PRIVATE already_known, known_potential, known_fourier_trafo
      34             : 
      35             :    ! functions needed internally
      36             :    PRIVATE calcYlm, calcSphBes, my_dot_product3x1, my_dot_product3x3, my_dot_product4x1, &
      37             :       my_dot_product4x4, my_sum3d, my_sum4d, gPtsSummation, approximateIntegral, generateIntegrals, &
      38             :       integrateY1ExpErfc, integrateY3ExpErfc, integrateY5ExpErfc, integrateY7ExpErfc, &
      39             :       integrateY9ExpErfc, calculateF, calculateG, calculateH
      40             : 
      41             :    INTERFACE my_sum
      42             :       MODULE PROCEDURE my_sum3d, my_sum4d
      43             :    END INTERFACE
      44             : 
      45             :    INTERFACE my_dot_product
      46             :       MODULE PROCEDURE my_dot_product4x1, my_dot_product4x4, my_dot_product3x1, my_dot_product3x3
      47             :    END INTERFACE
      48             : 
      49             :    INTERFACE gPtsSummation
      50             :       MODULE PROCEDURE crc_gPtsSummation, rrr_gPtsSummation
      51             :    END INTERFACE
      52             : 
      53             : CONTAINS
      54             : 
      55             :    ! Calculate the enhancement factor of the HSE03 functional
      56             :    ! References:
      57             :    ! [1] Heyd, Scuseria, Ernzerhof: Hybrid functionals based on a screened Coulomb potential,
      58             :    !     J. Chem. Phys. 118 (2003) 8207-8215
      59             :    ! [2] Heyd, Scuseria: Assessment and validation of a screened Coulomb hybrid density
      60             :    !     functional, J. Chem. Phys. 120 (2004) 7274-7280
      61             :    ! [3] Ernzerhof, Perdew: Generalized gradient approximation to the angle- and system-
      62             :    !     averaged exchange hole, J. Chem. Phys. 109 (1998) 3313-3319
      63             :    ! Input:  kF       - kF = (3 pi^2 rho) ** 1/3
      64             :    !         s_inp    - reduced gradient = |grad rho| / 2 kF rho
      65             :    ! Output: F_x      - HSE03 enhancement factor
      66             :    !         dFx_ds   - derivative of this factor with respect to s
      67             :    !         d2Fx_ds2 - second derivative with respect to s
      68           0 :    SUBROUTINE calculateEnhancementFactor(kF, s_inp, F_x, dFx_Ds, d2Fx_Ds2, dFx_dkF, d2Fx_dsdkF)
      69             : 
      70             :       IMPLICIT NONE
      71             : 
      72             :       REAL, INTENT(IN)  :: kF, s_inp
      73             :       REAL, INTENT(OUT) :: F_x, dFx_ds, d2Fx_ds2
      74             :       REAL, INTENT(OUT) :: dFx_dkF, d2Fx_dsdkF
      75             : 
      76             :       ! Helper variables
      77             :       REAL :: r1_kF, &                   ! 1 / kF
      78             :               omega_kF, omega_kF_Sqr, &      ! omega_HSE / kF, (omega/kF)^2
      79             :               s, s2, &                      ! reduced gradient s and s^2
      80             :               s_si2, r2s_si3, r6s_si4, &    ! quotient s_chg / s_inp^n (n = 2,3,4)
      81             :               F, G, H, &                    ! results of the functions F, G, and H defined in [3]
      82             :               Hs2, dHs2_ds, d2Hs2_ds2, &    ! H * s^2 and derivatives with respect to s
      83             :               Fs2, dFs2_ds, d2Fs2_ds2, &    ! F * s^2 and derivatives with respect to s
      84             :               dGs2_ds, d2Gs2_ds2, &         ! derivatives with respect to s of G s^2
      85             :               C_term, dCt_ds, d2Ct_ds2, &   ! C * (1 + F(s) s^2) and derivatives with respect to s
      86             :               E_term, dEt_ds, d2Et_ds2, &   ! E * (1 + G(s) s^2) and derivatives with respect to s
      87             :               D_term                        ! D + H s^2
      88             : 
      89             :       ! integrals of y^n Exp(-ay^2) Erfc(w/k y) for n = 1, 3, 5, 7, 9
      90             :       REAL :: intY1ExpErfc, intY3ExpErfc, intY5ExpErfc, intY7ExpErfc, intY9ExpErfc
      91             :       ! integrals of 2/sqrt(pi) y^n Exp(-arg*y^2) for n = 2, 4, 6, 8
      92             :       REAL :: r1_arg, intY2Gauss, intY4Gauss, intY6Gauss, intY8Gauss
      93             : 
      94             :       ! approximation of the not analytical part of the integral and derivatives
      95             :       REAL :: appInt, dAppInt_ds, d2AppInt_ds2
      96             :       REAL :: dAppInt_dkF, d2AppInt_dsdkF
      97             : 
      98             :       REAL, PARAMETER :: r8_9 = 8.0/9.0  ! 8/9
      99             :       REAL, PARAMETER :: &               ! Parameters of the exchange hole as defined in [3]
     100             :          B = -0.37170836, C = -0.077215461, &
     101             :          D = 0.57786348, E = -0.051955731
     102             : 
     103             :       REAL, PARAMETER :: &               ! Correction of the reduced gradient to ensure Lieb-Oxford bound [2]
     104             :          s_thresh = 8.3, s_max = 8.572844, s_chg = 18.796223
     105             :       LOGICAL :: correction              ! Has the value of s been corrected
     106             : 
     107             :       ! If a large value of s would violate the Lieb-Oxford bound, the value of s is reduced,
     108             :       ! so that this condition is fullfilled
     109           0 :       correction = s_inp > s_thresh
     110           0 :       IF (correction) THEN
     111           0 :          s_si2 = s_chg/(s_inp*s_inp)
     112           0 :          s = s_max - s_si2
     113             :       ELSE
     114           0 :          s = s_inp
     115             :       END IF
     116             : 
     117             :       ! Calculate different helper variables
     118           0 :       r1_kF = 1.0/kF
     119           0 :       omega_kF = omega_hse*r1_kF !was omega_VHSE()
     120           0 :       omega_kF_Sqr = omega_kF*omega_kF
     121             : 
     122             :       ! calculate the functions H and F in [3] and its derivatives
     123           0 :       CALL calculateH(s, H, dHs2_ds, d2Hs2_ds2)
     124           0 :       CALL calculateF(s, H, dHs2_ds, d2Hs2_ds2, F, dFs2_ds, d2Fs2_ds2)
     125           0 :       s2 = s*s
     126           0 :       Hs2 = H*s2
     127           0 :       Fs2 = F*s2
     128             :       CALL calculateG(s2, Fs2, dFs2_ds, d2Fs2_ds2, Hs2, dHs2_ds, d2Hs2_ds2, & !Input
     129           0 :                       G, dGs2_ds, d2Gs2_ds2)                                 !Output
     130             : 
     131           0 :       C_term = C*(1 + s2*F)
     132           0 :       dCt_ds = C*dFs2_ds
     133           0 :       d2Ct_ds2 = C*d2Fs2_ds2
     134           0 :       E_term = E*(1 + s2*G)
     135           0 :       dEt_ds = E*dGs2_ds
     136           0 :       d2Et_ds2 = E*d2Gs2_ds2
     137           0 :       D_term = D + Hs2
     138             : 
     139             :       ! approximate the integral using an expansion of the error function (cf. [2])
     140             :       CALL approximateIntegral(omega_kF, Hs2, D_term, dHs2_ds, d2Hs2_ds2, & !Input
     141           0 :                                appInt, dAppInt_ds, d2AppInt_ds2, dAppInt_dkF, d2AppInt_dsdkF)  !Output
     142             : 
     143             :       ! Calculate the integrals
     144             :       !
     145             :       ! inf
     146             :       !   /                2   2      /       \
     147             :       !  |     n  -(D+H(s)s ) y      | omega   |
     148             :       !  | dy y  e               Erfc| ----- y |
     149             :       !  |                           |   k     |
     150             :       ! /                             \   F   /
     151             :       !  0
     152             :       ! for n = 1, 3, 5, 7, and 9
     153             :       !
     154           0 :       intY1ExpErfc = integrateY1ExpErfc(omega_kf, omega_kf_Sqr, D_term)
     155           0 :       intY3ExpErfc = integrateY3ExpErfc(omega_kf, omega_kf_Sqr, D_term)
     156           0 :       intY5ExpErfc = integrateY5ExpErfc(omega_kf, omega_kf_Sqr, D_term)
     157           0 :       intY7ExpErfc = integrateY7ExpErfc(omega_kf, omega_kf_Sqr, D_term)
     158           0 :       intY9ExpErfc = integrateY9ExpErfc(omega_kf, omega_kf_Sqr, D_term)
     159             : 
     160             :       ! Calculate the integrals
     161             :       !
     162             :       !       inf
     163             :       !         /          /                      2    \
     164             :       !   2    |     n    |          2   2   omega   2  |
     165             :       ! -----  | dy y  Exp| -(D+H(s)s ) y  - ------ y   |
     166             :       !  ____  |          |                   k ^2      |
     167             :       ! V Pi  /            \                   F       /
     168             :       !        0
     169             :       ! for n = 2, 4, 6, 8
     170             :       !
     171           0 :       r1_arg = 1.0/(D_term + omega_kF_Sqr)
     172           0 :       intY2Gauss = 0.5*SQRT(r1_arg)*r1_arg
     173           0 :       intY4Gauss = 1.5*intY2Gauss*r1_arg
     174           0 :       intY6Gauss = 2.5*intY4Gauss*r1_arg
     175           0 :       intY8Gauss = 3.5*intY6Gauss*r1_arg
     176             : 
     177             :       ! Calculate the integral
     178             :       !  inf
     179             :       !    /                 /       \
     180             :       !   |                 | omega   |
     181             :       !   | dy y J(s,y) Erfc| ----- y |
     182             :       !   |                 |   k     |
     183             :       !  /                   \   F   /
     184             :       !   0
     185             :       ! where J(s, y) is the exchange hole defined in [3]
     186             :       ! the exchange factor is proportional to this integral
     187           0 :       F_x = -r8_9*(appInt + B*intY1ExpErfc + C_term*intY3ExpErfc + E_term*intY5ExpErfc)
     188             : 
     189             :       ! Calculate the derivatives with respect to s using that the derivatative of the integral
     190             :       ! yields higher orders of the same kind of integral intY1 -> -intY3 -> intY5 ... times
     191             :       ! the derivative of the exponent
     192             :       dFx_ds = -r8_9*(dAppInt_ds - (B*intY3ExpErfc + C_term*intY5ExpErfc + E_term*intY7ExpErfc)*dHs2_ds &
     193           0 :                       + dCt_ds*intY3ExpErfc + dEt_ds*intY5ExpErfc)
     194             :       d2Fx_ds2 = -r8_9*(d2AppInt_ds2 + (B*intY5ExpErfc + C_term*intY7ExpErfc + E_term*intY9ExpErfc)*dHs2_ds**2 &
     195             :                         - (B*intY3ExpErfc + C_term*intY5ExpErfc + E_term*intY7ExpErfc)*d2Hs2_ds2 &
     196             :                         - 2.0*(dCt_ds*intY5ExpErfc + dEt_ds*intY7ExpErfc)*dHs2_ds &
     197           0 :                         + d2Ct_ds2*intY3ExpErfc + d2Et_ds2*intY5ExpErfc)
     198             : 
     199           0 :       dFx_dkF = -r8_9*r1_kF*(omega_kF*(B*intY2Gauss + C_term*intY4Gauss + E_term*intY6Gauss) + dAppInt_dkF)
     200             :       d2Fx_dsdkF = -r8_9*r1_kF*(d2AppInt_dsdkF + omega_kF*(dCt_ds*intY4Gauss + dEt_ds*intY6Gauss &
     201           0 :                                                            - (B*intY4Gauss + C_term*intY6Gauss + E_term*intY8Gauss)*dHs2_ds))
     202             : 
     203             :       ! Correction to the derivatives, if s_inp > s_thresh
     204           0 :       IF (correction) THEN
     205           0 :          r2s_si3 = 2.0*s_si2/s_inp
     206           0 :          r6s_si4 = 3.0*r2s_si3/s_inp
     207           0 :          d2Fx_ds2 = d2Fx_ds2*r2s_si3**2 - dFx_ds*r6s_si4
     208           0 :          d2Fx_dsdkF = d2Fx_dsdkF*r2s_si3
     209           0 :          dFx_ds = dFx_ds*r2s_si3
     210             :       END IF
     211             : 
     212           0 :    END SUBROUTINE calculateEnhancementFactor
     213             : 
     214             :    ! Approximation for the first part of the exchange hole integral
     215             :    ! Calculated using the algorithm described in [2].
     216             :    ! Additionally the first and second derivative of this approximation with
     217             :    ! respect to s are calculated
     218             :    ! Input:  omega_kF       - omega / kF
     219             :    !         Hs2            - H s^2
     220             :    !         D_Hs2          - D + H s^2
     221             :    !         dHs2_ds        - d (H s^2) / ds
     222             :    !         d2Hs2_ds2      - d^2 (H s^2) / ds^2
     223             :    ! Output: appInt         - approximation for one part of the exchange hole integral
     224             :    !         dAppInt_ds     - first derivative with respect to s
     225             :    !         d2AppInt_ds2   - second derivative with respect to s
     226             :    !         dAppInt_dkF    - first derivative with respect to kF
     227             :    !         d2AppInt_dsdkF - mixed derivative with respect to s and kF
     228           0 :    SUBROUTINE approximateIntegral(omega_kF, Hs2, D_Hs2, dHs2_ds, d2Hs2_ds2, &
     229             :                                   appInt, dAppInt_ds, d2AppInt_ds2, dAppInt_dkF, d2AppInt_dsdkF)
     230             : 
     231             :       USE m_exponential_integral, ONLY: calculateExponentialIntegral, gauss_laguerre
     232             : 
     233             :       IMPLICIT NONE
     234             : 
     235             :       REAL, INTENT(IN)  :: omega_kF, Hs2, D_Hs2, dHs2_ds, d2Hs2_ds2
     236             :       REAL, INTENT(OUT) :: appInt, dAppInt_ds, d2AppInt_ds2
     237             :       REAL, INTENT(OUT) :: dAppInt_dkF, d2AppInt_dsdkF
     238             : 
     239             :       REAL    :: w2, bw, r2bw, bw_Hs2, bw_D_Hs2
     240             :       ! variables for temporary storage of the integrals, the prefactors and the dot_product
     241             :       REAL    :: integral(0:12), a_omegaI(0:8), aI_omegaI(0:8), dotpr, dotpr2
     242             :       INTEGER :: i
     243             : 
     244             :       REAL    :: r2w2, r4w2, r2w2_Hs2, r2w2_D_Hs2, arg, exp_e1, dAppInt_dh, d2AppInt_dh2
     245             : 
     246             :       ! parameters of the erfc fit as given in [2]
     247             :       REAL, PARAMETER :: A_2 = 0.5080572, scale = 1.125/A_2, &
     248             :                          a(1:8) = (/-1.128223946706117, 1.452736265762971, &
     249             :                                     -1.243162299390327, 0.971824836115601, &
     250             :                                     -0.568861079687373, 0.246880514820192, &
     251             :                                     -0.065032363850763, 0.008401793031216/), &
     252             :                          b = 1.455915450052607, cutoff = 14.0
     253             : 
     254             :       ! Calculate helper variables
     255           0 :       w2 = omega_kF**2
     256           0 :       bw = b*w2
     257           0 :       r2bw = 2.0*bw
     258           0 :       bw_Hs2 = bw + Hs2
     259           0 :       bw_D_Hs2 = bw + D_Hs2
     260             : 
     261           0 :       IF (bw_Hs2 < cutoff) THEN
     262             :          ! integrate the erfc approximation times the exchange hole
     263           0 :          CALL generateIntegrals(bw_Hs2, bw_D_Hs2, integral)
     264             : 
     265             :          ! combine the solutions of the integrals with the appropriate prefactors
     266           0 :          a_omegaI(0) = 1.0
     267           0 :          a_omegaI(1:8) = a*(/(omega_kF**i, i=1, 8)/)
     268           0 :          aI_omegaI = (/(i*a_omegaI(i), i=0, 8)/)
     269           0 :          appInt = DOT_PRODUCT(a_omegaI, integral(0:8))
     270           0 :          dotpr = DOT_PRODUCT(a_omegaI, integral(2:10))
     271           0 :          dAppInt_ds = -dotpr*dHs2_ds
     272           0 :          dAppInt_dkF = r2bw*dotpr - DOT_PRODUCT(aI_omegaI, integral(0:8))
     273           0 :          dotpr2 = DOT_PRODUCT(a_omegaI, integral(4:12))
     274           0 :          d2AppInt_ds2 = dotpr2*dHs2_ds**2 - dotpr*d2Hs2_ds2
     275           0 :          d2AppInt_dsdkF = -(r2bw*dotpr2 - DOT_PRODUCT(aI_omegaI, integral(2:10)))*dHs2_ds
     276             :       ELSE
     277           0 :          r2w2 = 2.0*w2
     278           0 :          r4w2 = 4.0*w2
     279           0 :          r2w2_Hs2 = r2w2 + Hs2
     280           0 :          r2w2_D_Hs2 = r2w2 + D_Hs2
     281           0 :          arg = scale*r2w2_Hs2
     282           0 :          exp_e1 = gauss_laguerre(arg)
     283           0 :          appInt = A_2*(LOG(r2w2_Hs2/r2w2_D_Hs2) + exp_e1)
     284           0 :          dAppInt_dh = -A_2/r2w2_D_Hs2 + 1.125*exp_e1
     285           0 :          d2AppInt_dh2 = 1.125/r2w2_Hs2 - A_2/r2w2_D_Hs2**2 - 1.265625/A_2*exp_e1
     286           0 :          dAppInt_ds = dAppInt_dh*dHs2_ds
     287           0 :          dAppInt_dkF = -dAppInt_dh*r4w2
     288           0 :          d2AppInt_ds2 = d2AppInt_dh2*dHs2_ds**2 + dAppInt_dh*d2Hs2_ds2
     289           0 :          d2AppInt_dsdkF = d2AppInt_dh2*dHs2_ds*r4w2
     290             :       END IF
     291             : 
     292           0 :    END SUBROUTINE approximateIntegral
     293             : 
     294             :    ! Generate the integrals for n = 0, 1, ..., 12
     295             :    !
     296             :    ! inf
     297             :    !   /        /        2                  \      /                        \
     298             :    !  |     n  |  A   -Dy           A        |    |    /   omega      2\   2 |
     299             :    !  | dy y   | --- e     - --------------  | Exp| - (  b ----- + H s  ) y  |
     300             :    !  |        |  y            /    4   2\   |    |    \     k         /     |
     301             :    ! /          \            y( 1 + - Ay  ) /      \          F             /
     302             :    !  0                        \    9    /
     303             :    !
     304             :    ! Input:  bw_Hs2   - b (omega/kF)^2 + H s^2
     305             :    !         bw_D_Hs2 - b (omega/kF)^2 + D + H s^2
     306             :    ! Output: integral - array with the calculated integrals
     307             :    ! To simplify the calculation use integral(n+2) = - d(integral(n))/d(b omega/kF)
     308           0 :    SUBROUTINE generateIntegrals(bw_Hs2, bw_D_Hs2, integral)
     309             :       USE m_exponential_integral, ONLY: calculateExponentialIntegral, gauss_laguerre, series_laguerre
     310             :       USE m_constants
     311             :       IMPLICIT NONE
     312             : 
     313             :       REAL, INTENT(IN)  :: bw_Hs2, bw_D_Hs2
     314             :       REAL, INTENT(OUT) :: integral(0:12)
     315             : 
     316             :       ! Helper variables
     317             :       REAL :: bw_Hs2_Sqr, bw_Hs2_Cub, sqrt_bw_Hs2, &
     318             :               bw_D_Hs2_Sqr, bw_D_Hs2_Cub, bw_D_Hs2_Tet, sqrt_bw_D_Hs2, &
     319             :               arg, sqrt_arg, r1_arg, e1_arg, exp_arg, exp_erfc, &
     320             :               term1, factor2, term2, arg_n, sum_term, add_term, half_i2_1
     321             :       INTEGER :: i
     322             : 
     323             :       ! A is an exchange hole parameter [3] and b the fit parameter for erfc [2]
     324             :       REAL, PARAMETER :: &
     325             :          A = 1.0161144, A_2 = A/2.0, scale = 2.25/A, sqrtA = 1.008025, & !sqrt(A)
     326             :          b = 1.455915450052607
     327             : 
     328             :       ! Calculate many helper variables
     329           0 :       bw_Hs2_Sqr = bw_Hs2*bw_Hs2
     330           0 :       bw_Hs2_Cub = bw_Hs2*bw_Hs2_Sqr
     331           0 :       sqrt_bw_Hs2 = SQRT(bw_Hs2)
     332           0 :       bw_D_Hs2_Sqr = bw_D_Hs2*bw_D_Hs2
     333           0 :       bw_D_Hs2_Cub = bw_D_Hs2*bw_D_Hs2_Sqr
     334           0 :       bw_D_Hs2_Tet = bw_D_Hs2_Sqr*bw_D_Hs2_Sqr
     335           0 :       sqrt_bw_D_Hs2 = SQRT(bw_D_Hs2)
     336             : 
     337             :       ! arg = 9/4 * (b omega/kF + H s^2) / A
     338           0 :       arg = scale*bw_Hs2
     339           0 :       sqrt_arg = SQRT(arg)
     340           0 :       r1_arg = 1.0/arg
     341             : 
     342             :       ! calculate e^(arg), E1(arg), and erfc(sqrt(arg))
     343           0 :       exp_arg = EXP(arg)
     344           0 :       exp_erfc = exp_arg*erfc(sqrt_arg)
     345             : 
     346           0 :       IF (arg > series_laguerre) THEN
     347           0 :          term2 = gauss_laguerre(arg)
     348             :       ELSE
     349           0 :          CALL calculateExponentialIntegral(arg, e1_arg)
     350           0 :          term2 = exp_arg*e1_arg
     351             :       END IF
     352             : 
     353             :       ! The n = 0 integral is
     354             :       ! A/2 ( ln((b omega/kF + H s^2) / (b omega/kF + D + H s^2))
     355             :       !     + e^(arg) E1(arg) )
     356           0 :       integral(0) = A_2*(LOG(bw_Hs2/bw_D_Hs2) + term2)
     357             : 
     358             :       ! Calculate now all even n's by successive derivation
     359             :       ! The log(...) term gives term proportional to 1/(b omega/kF + D + H s^2)^i
     360             :       ! The e^(arg) E1(arg) term reproduces itself with a prefactor and
     361             :       ! generates an additional 1/arg term which produces higher 1/arg^i terms
     362             :       ! when further deriviated
     363           0 :       term1 = A_2/bw_D_Hs2
     364           0 :       factor2 = -1.125
     365           0 :       arg_n = -1.0/arg
     366           0 :       integral(2) = term1 + factor2*term2
     367             : 
     368           0 :       DO i = 1, 5
     369           0 :          term1 = term1/bw_D_Hs2*REAL(i)
     370           0 :          factor2 = -factor2*scale
     371           0 :          term2 = term2 + arg_n
     372             : 
     373           0 :          integral(2*(i + 1)) = term1 + factor2*term2
     374             : 
     375           0 :          arg_n = -arg_n*REAL(i)/arg
     376             :       END DO
     377             : 
     378             :       ! The n = 1 integral is
     379             :       ! A/2 ( sqrt(pi) / sqrt( b omega/kF + D + H s2 )
     380             :       ! - 3/4 sqrt(A) pi e^(arg) erfc(sqrt(arg))
     381           0 :       term1 = A_2*SQRT(PI_const)/sqrt_bw_D_Hs2
     382           0 :       term2 = PI_const*exp_erfc
     383           0 :       factor2 = -0.75*sqrtA
     384             : 
     385           0 :       integral(1) = term1 + factor2*term2
     386             : 
     387             :       ! Calculate now all uneven n's by successive derivation
     388             :       ! The 1 / sqrt(...) term gives higher orders of 1 / (...)^((2i+1)/2)
     389             :       ! The e^(arg) erfc(sqrt(arg)) term reproduces itself with a prefactor
     390             :       ! and generates an additional 1/sqrt(arg) term which produces higher
     391             :       ! 1/(arg)^((2i+1)/2) terms when further deriviated
     392           0 :       sum_term = -1.125*SQRT(pi_const)/sqrt_bw_Hs2
     393           0 :       add_term = sum_term
     394           0 :       half_i2_1 = -0.5
     395             : 
     396           0 :       DO i = 3, 11, 2
     397           0 :          factor2 = -factor2*scale
     398           0 :          term1 = -term1*half_i2_1/bw_D_Hs2
     399           0 :          integral(i) = term1 + term2*factor2 + sum_term
     400             : 
     401           0 :          add_term = -add_term*half_i2_1/bw_Hs2
     402           0 :          sum_term = -sum_term*scale + add_term
     403           0 :          half_i2_1 = half_i2_1 - 1.0
     404             :       ENDDO
     405             : 
     406           0 :    END SUBROUTINE generateIntegrals
     407             : 
     408             :    ! Calculate the integral
     409             :    !
     410             :    ! inf
     411             :    !   /                  /       \         /                 ________________ \
     412             :    !  |       -a y^2     | omega   |    1  |     omega   /   /     / omega \2   |
     413             :    !  | dy y e       Erfc| ----- y | = --- | 1 - -----  /   / a + (  -----  )   |
     414             :    !  |                  |   k     |   2 a |       k   /   V       \   k   /    |
     415             :    ! /                    \   F   /         \       F                   F      /
     416             :    !  0
     417             :    !
     418             :    ! Input:  omega_kF     - omega / kF
     419             :    !         omega_kF_Sqr - (omega / kF)^2
     420             :    !         a            - factor of the gaussian function
     421             :    ! Return: result of the integration
     422             :    REAL FUNCTION integrateY1ExpErfc(omega_kF, omega_kF_Sqr, a)
     423             :       IMPLICIT NONE
     424             : 
     425             :       REAL, INTENT(IN)  :: omega_kF, omega_kF_Sqr, a
     426             : 
     427             :       ! helper variable sqrt(a + (omega/kF)^2)
     428             :       REAL :: sqrt_term
     429             : 
     430             :       ! calculate helper variable
     431           0 :       sqrt_term = SQRT(a + omega_kF_Sqr)
     432             : 
     433             :       ! calculate the integral
     434           0 :       integrateY1ExpErfc = 0.5*(1 - omega_kF/sqrt_term)/a
     435             : 
     436             :    END FUNCTION integrateY1ExpErfc
     437             : 
     438             :    ! Calculate the integral
     439             :    !
     440             :    ! inf
     441             :    !   /                   /       \           /                   ________________                   ________________3 \
     442             :    !  |     3  -a y^2     | omega   |     1   |       omega   /   /     / omega \2       omega   /   /     / omega \2    |
     443             :    !  | dy y  e       Erfc| ----- y | = ----- | 2 - 2 -----  /   / a + (  -----  )  - a  -----  /   / a + (  -----  )    |
     444             :    !  |                   |   k     |   4 a^2 |         k   /   V       \   k   /          k   /   V       \   k   /     |
     445             :    ! /                     \   F   /           \         F                   F              F                   F       /
     446             :    !  0
     447             :    !
     448             :    ! Input:  omega_kF     - omega / kF
     449             :    !         omega_kF_Sqr - (omega / kF)^2
     450             :    !         a            - factor of the gaussian function
     451             :    ! Return: result of the integration
     452           0 :    REAL FUNCTION integrateY3ExpErfc(omega_kF, omega_kF_Sqr, a)
     453             :       IMPLICIT NONE
     454             : 
     455             :       REAL, INTENT(IN)  :: omega_kF, omega_kF_Sqr, a
     456             : 
     457             :       ! helper variables
     458             :       REAL :: term, sqrt_term, sqrt_term_3
     459             : 
     460             :       ! calculate helper variables
     461           0 :       term = a + omega_kF_Sqr
     462           0 :       sqrt_term = SQRT(term)
     463           0 :       sqrt_term_3 = term*sqrt_term
     464             : 
     465             :       ! calculate the integral
     466           0 :       integrateY3ExpErfc = 0.25*(2.0 - a*omega_kF/sqrt_term_3 - 2.0*omega_kF/sqrt_term)/(a*a)
     467             : 
     468           0 :    END FUNCTION integrateY3ExpErfc
     469             : 
     470             :    ! Calculate the integral
     471             :    !
     472             :    ! inf                                     /                                   \
     473             :    !   /                   /       \        |                  /              2\  |
     474             :    !  |     5  -a y^2     | omega   |    1  |       omega     |      a     3 a  | |
     475             :    !  | dy y  e       Erfc| ----- y | = --- | 1 - ----------  | 1 + --- + ----- | |
     476             :    !  |                   |   k     |   a^3 |     k  sqrt(x)  |     2 x   8 x^2 | |
     477             :    ! /                     \   F   /        |      F           \               /  |
     478             :    !  0                                      \                                   /
     479             :    ! where
     480             :    !          /     \ 2
     481             :    !         | omega |
     482             :    ! x = a + | ----- |
     483             :    !         |   k   |
     484             :    !          \   F /
     485             :    !
     486             :    ! Input:  omega_kF     - omega / kF
     487             :    !         omega_kF_Sqr - (omega / kF)^2
     488             :    !         a            - factor of the gaussian function
     489             :    ! Return: result of the integration
     490           0 :    REAL FUNCTION integrateY5ExpErfc(omega_kF, omega_kF_Sqr, a)
     491             :       IMPLICIT NONE
     492             : 
     493             :       REAL, INTENT(IN)  :: omega_kF, omega_kF_Sqr, a
     494             : 
     495             :       ! helper variables
     496             :       REAL :: term, sqrt_term, a_Sqr
     497             : 
     498             :       ! calculate helper variables
     499           0 :       term = a + omega_kF_Sqr
     500           0 :       sqrt_term = SQRT(term)
     501           0 :       a_Sqr = a*a
     502             : 
     503             :       ! calculate the integral
     504           0 :       integrateY5ExpErfc = (1.0 - (0.375*a_Sqr/(term*term) + a/(2.0*term) + 1.0)*omega_kF/sqrt_term)/(a_Sqr*a)
     505             : 
     506           0 :    END FUNCTION integrateY5ExpErfc
     507             : 
     508             :    ! Calculate the integral
     509             :    !
     510             :    ! inf                                     /                                            \
     511             :    !   /                   /       \        |                  /              2       3 \  |
     512             :    !  |     7  -a y^2     | omega   |    3  |       omega     |      a     3 a     5 a   | |
     513             :    !  | dy y  e       Erfc| ----- y | = --- | 1 - ----------  | 1 + --- + ----- + ------ | |
     514             :    !  |                   |   k     |   a^4 |     k  sqrt(x)  |     2 x   8 x^2   16 x^3 | |
     515             :    ! /                     \   F   /        |      F           \                        /  |
     516             :    !  0                                      \                                            /
     517             :    ! where
     518             :    !          /     \ 2
     519             :    !         | omega |
     520             :    ! x = a + | ----- |
     521             :    !         |   k   |
     522             :    !          \   F /
     523             :    !
     524             :    ! Input:  omega_kF     - omega / kF
     525             :    !         omega_kF_Sqr - (omega / kF)^2
     526             :    !         a            - factor of the gaussian function
     527             :    ! Return: result of the integration
     528           0 :    REAL FUNCTION integrateY7ExpErfc(omega_kF, omega_kF_Sqr, a)
     529             :       IMPLICIT NONE
     530             : 
     531             :       REAL, INTENT(IN)  :: omega_kF, omega_kF_Sqr, a
     532             : 
     533             :       ! helper variables
     534             :       REAL :: term2, term, sqrt_term, a_Sqr
     535             : 
     536             :       ! calculate helper variables
     537           0 :       term = a + omega_kF_Sqr
     538           0 :       term2 = term*term
     539           0 :       sqrt_term = SQRT(term)
     540           0 :       a_Sqr = a*a
     541             : 
     542             :       ! calculate the integral
     543             :       integrateY7ExpErfc = 3.0/(a_Sqr*a_Sqr)*(1.0 - omega_kF/sqrt_term* &
     544           0 :                                               (0.3125*(a_Sqr*a)/(term2*term) + 0.375*a_Sqr/term2 + a/(2.0*term) + 1.0))
     545             : 
     546           0 :    END FUNCTION integrateY7ExpErfc
     547             : 
     548             :    ! Calculate the integral
     549             :    !
     550             :    ! inf                                     /                                                      \
     551             :    !   /                   /       \        |                  /              2       3           \  |
     552             :    !  |     9  -a y^2     | omega   |    12 |       omega     |      a     3 a     5 a      35 a^4 | |
     553             :    !  | dy y  e       Erfc| ----- y | = --- | 1 - ----------  | 1 + --- + ----- + ------ + ------- | |
     554             :    !  |                   |   k     |   a^5 |     k  sqrt(x)  |     2 x   8 x^2   16 x^3   128 x^4 | |
     555             :    ! /                     \   F   /        |      F           \                                  /  |
     556             :    !  0                                      \                                                      /
     557             :    ! where
     558             :    !          /     \ 2
     559             :    !         | omega |
     560             :    ! x = a + | ----- |
     561             :    !         |   k   |
     562             :    !          \   F /
     563             :    !
     564             :    ! Input:  omega_kF     - omega / kF
     565             :    !         omega_kF_Sqr - (omega / kF)^2
     566             :    !         a            - factor of the gaussian function
     567             :    ! Return: result of the integration
     568           0 :    REAL FUNCTION integrateY9ExpErfc(omega_kF, omega_kF_Sqr, a)
     569             :       IMPLICIT NONE
     570             : 
     571             :       REAL, INTENT(IN)  :: omega_kF, omega_kF_Sqr, a
     572             : 
     573             :       ! helper variables
     574             :       REAL :: term2, term, sqrt_term, a_Sqr, a_Tet
     575             : 
     576             :       ! calculate helper variables
     577           0 :       term = a + omega_kF_Sqr
     578           0 :       term2 = term*term
     579           0 :       sqrt_term = SQRT(term)
     580           0 :       a_Sqr = a*a
     581           0 :       a_Tet = a_Sqr*a_Sqr
     582             : 
     583             :       ! calculate the integral
     584             :       integrateY9ExpErfc = 12.0*(1.0 - (0.2734375*a_Tet/(term2*term2) + 0.3125*(a_Sqr*a)/(term2*term) &
     585           0 :                                         + 0.375*a_Sqr/term2 + a/(2.0*term) + 1.0)*omega_kF/sqrt_term)/(a_Tet*a)
     586             : 
     587           0 :    END FUNCTION integrateY9ExpErfc
     588             : 
     589             :    ! Calculate the function F(s) given in [3]
     590             :    ! Input:  s         - reduced gradient
     591             :    !         H         - value of the function H(s)
     592             :    !         dHs2_ds   - first derivative of H(s)s^2
     593             :    !         dsHs2_ds2 - second derivative of H(s)s^2
     594             :    ! Output: F         - value of the function F
     595             :    !         dFs2_ds   - first derivative of F(s)s^2
     596             :    !         d2Fs2_ds2 - second derivative of F(s)s^2
     597             :    SUBROUTINE calculateF(s, H, dHs2_ds, d2Hs2_ds2, F, dFs2_ds, d2Fs2_ds2)
     598             :       IMPLICIT NONE
     599             : 
     600             :       REAL, INTENT(IN)  :: s, H, dHs2_ds, d2Hs2_ds2
     601             :       REAL, INTENT(OUT) :: F, dFs2_ds, d2Fs2_ds2
     602             : 
     603             :       REAL, PARAMETER   :: slope = 6.4753871
     604             :       REAL, PARAMETER   :: shift = 0.4796583
     605             : 
     606             :       ! calculate the function F(s), d(s^2F(s))/ds, and d^2(s^2F(s))/ds^2
     607           0 :       F = slope*H + shift
     608           0 :       dFs2_ds = slope*dHs2_ds + 2.0*s*shift
     609           0 :       d2Fs2_ds2 = slope*d2Hs2_ds2 + 2.0*shift
     610             :    END SUBROUTINE calculateF
     611             : 
     612             :    ! Calculate the function G(s) given in [3]
     613             :    ! Input:  s2        - s^2 where s is the reduced gradient
     614             :    !         Fs2       - F(s) s^2
     615             :    !         dFs2_ds   - first derivative of F(s)s^2 with respect to s
     616             :    !         d2Fs2_ds2 - second derivative of F(s)s^2
     617             :    !         Hs2       - H(s) s^2
     618             :    !         dHs2_ds   - first derivative of H(s)s^2 with respect to s
     619             :    !         d2Hs2_ds2 - second derivative of H(s)s^2
     620             :    ! Output: G         - value of the function G
     621             :    !         dGs2_ds   - first derivative of G(s)s^2 with respect to s
     622             :    !         d2Gs2_ds2 - second derivative of G(s)s^2
     623           0 :    SUBROUTINE calculateG(s2, Fs2, dFs2_ds, d2Fs2_ds2, Hs2, dHs2_ds, d2Hs2_ds2, G, dGs2_ds, d2Gs2_ds2)
     624             :       IMPLICIT NONE
     625             : 
     626             :       REAL, INTENT(IN)  :: s2, Fs2, dFs2_ds, d2Fs2_ds2, Hs2, dHs2_ds, d2Hs2_ds2
     627             :       REAL, INTENT(OUT) :: G, dGs2_ds, d2Gs2_ds2
     628             : 
     629             :       ! helper variables
     630             :       REAL :: AHs2_1_2, AHs2_3_2, r1_Fs2, D_Hs2, D_Hs2Sqr, D_Hs2Cub, &
     631             :               D_Hs2_5_2, D_Hs2_7_2, D_Hs2_9_2, D_Hs2_11_2
     632             :       REAL :: part1, dpart1_dh, d2part1_dh2
     633             :       REAL :: arg1, arg2, exp_erfc, &
     634             :               part2, dpart2_dh, d2part2_dh2
     635             :       REAL :: alpha, Ebeta, r3Pi_4_alpha, beta_s2, Ebeta_s2, &
     636             :               dalpha_dh, d2alpha_dh2, dalpha_df, d2alpha_dfdh, dbeta_dh, d2beta_dh2
     637             : 
     638             :       ! parameters of the exchange hole given in [3]
     639             :       REAL, PARAMETER :: PI_75 = 2.356194490192344929, SQRT_PI = 1.77245385090551602729816748334, &
     640             :                          A = 1.0161144, sqrtA = 1.008025, r9_4A = 2.25/A, & !sqrt(A), 9/4A
     641             :                          B = -0.37170836, &
     642             :                          C = -0.077215461, &
     643             :                          D = 0.57786348, &
     644             :                          E = -0.051955731
     645             : 
     646             :       ! calculate the helper variables
     647           0 :       AHs2_1_2 = sqrtA*SQRT(Hs2)
     648           0 :       AHs2_3_2 = AHs2_1_2*A*Hs2
     649           0 :       r1_Fs2 = 1.0 + Fs2
     650           0 :       D_Hs2 = D + Hs2
     651           0 :       D_Hs2Sqr = D_Hs2*D_Hs2
     652           0 :       D_Hs2Cub = D_Hs2*D_Hs2Sqr
     653           0 :       D_Hs2_5_2 = D_Hs2Sqr*SQRT(D_Hs2)
     654           0 :       D_Hs2_7_2 = D_Hs2_5_2*D_Hs2
     655           0 :       D_Hs2_9_2 = D_Hs2_5_2*D_Hs2Sqr
     656           0 :       D_Hs2_11_2 = D_Hs2_5_2*D_Hs2Cub
     657             : 
     658             :       ! calculate first part of the term called 'a' in [3] eq. (A2) and its derivatives with respect to H(s)s^2 and F(s)s^2
     659             :       !
     660             :       !                          2      2           2 2          2 3
     661             :       !          __ 15E + 6C(1+Fs )(D+Hs ) + 4B(D+Hs )  + 8A(D+Hs )
     662             :       ! part1 = VPi ------------------------------------------------
     663             :       !                                    2  7/2
     664             :       !                        16 ( D + H s  )
     665             :       !
     666           0 :       part1 = SQRT_PI*(15.0*E + 6.0*C*r1_Fs2*D_Hs2 + 4.0*B*D_Hs2Sqr + 8.0*A*D_Hs2Cub)/(16.0*D_Hs2_7_2)
     667             :       !
     668             :       !                               2      2            2 2          2 3
     669             :       ! d part1     __ 105E + 30C(1+Fs )(D+Hs ) + 12B(D+Hs )  + 8A(D+Hs )
     670             :       ! ------- = -VPi ---------------------------------------------------
     671             :       ! d(Hs^2)                                  2  9/2
     672             :       !                              32 ( D + H s  )
     673             :       !
     674           0 :       dpart1_dh = -SQRT_PI*(105.0*E + 30.0*C*r1_Fs2*D_Hs2 + 12.0*B*D_Hs2Sqr + 8.0*A*D_Hs2Cub)/(32.0*D_Hs2_9_2)
     675             :       !
     676             :       !  2                             2      2            2 2           2 3
     677             :       ! d  part1    __ 945E + 210C(1+Fs )(D+Hs ) + 60B(D+Hs )  + 24A(D+Hs )
     678             :       ! -------- = VPi -----------------------------------------------------
     679             :       !        2                                   2  11/2
     680             :       ! d(Hs^2)                        64 ( D + H s  )
     681             :       !
     682           0 :       d2part1_dh2 = SQRT_PI*(945.0*E + 210.0*C*r1_Fs2*D_Hs2 + 60.0*B*D_Hs2Sqr + 24.0*A*D_Hs2Cub)/(64.0*D_Hs2_11_2)
     683             :       !
     684             :       ! d part1    __       3 C
     685             :       ! ------- = VPi -----------------
     686             :       ! d(Fs^2)                  2  5/2
     687             :       !               8 ( D + H s  )
     688             :       !
     689           0 :       dalpha_df = SQRT_PI*0.375*C/D_Hs2_5_2
     690             :       !
     691             :       !      2
     692             :       !     d part1         __       15 C
     693             :       ! --------------- = -VPi ------------------
     694             :       ! d(Fs^2) d(Hs^2)                    2  7/2
     695             :       !                        16 ( D + H s  )
     696             :       !
     697           0 :       d2alpha_dfdh = -2.5*dalpha_df/D_Hs2
     698             : 
     699             :       ! calculate second part of the term called 'a' in [3] eq. (A2) and its derivatives
     700             :       !                                            ________
     701             :       !                       /     2 \       /   /     2   \
     702             :       !         3 Pi  ___    | 9 H s   |     |   / 9 H s     |
     703             :       ! part2 = ---- V A  Exp| ------- | Erfc|  /  -------   |
     704             :       !          4           |   4 A   |     | V     4 A     |
     705             :       !                       \       /       \             /
     706             :       !
     707           0 :       arg1 = r9_4A*Hs2
     708           0 :       arg2 = SQRT(arg1)
     709           0 :       exp_erfc = EXP(arg1)*erfc(arg2)
     710           0 :       part2 = PI_75*sqrtA*exp_erfc
     711             :       !
     712             :       !                      /                                \
     713             :       ! d part2   3 Pi  ___ |  9                    3          |
     714             :       ! ------- = ---- V A  | --- exp_erfc - ----------------  |
     715             :       ! d(Hs^2)    4        | 4 A                          2   |
     716             :       !                      \               2 sqrt(Pi A Hs ) /
     717             :       !
     718           0 :       dpart2_dh = PI_75*sqrtA*(r9_4A*exp_erfc - 1.5/(SQRT_PI*AHs2_1_2))
     719             :       !
     720             :       !  2                    /                                 2    \
     721             :       ! d  part2   3 Pi  ___ |   81                  6 A - 27 Hs      |
     722             :       ! -------- = ---- V A  | ------ exp_erfc + -------------------- |
     723             :       !        2    4        | 16 A^2                           2 3/2 |
     724             :       ! d(Hs^2)               \                  8 sqrt(Pi)(A Hs )   /
     725             :       !
     726           0 :       d2part2_dh2 = PI_75*sqrtA*(r9_4A**2*exp_erfc + (6.0*A - 27.0*Hs2)/(8.0*SQRT_PI*AHs2_3_2))
     727             : 
     728             :       ! calculate the 'a' and 'b' terms ([3] eq. (A2) and (A3) )
     729             :       !
     730             :       ! a = part1 - part2
     731             :       !
     732             :       !             __  2
     733             :       !         15 VPi s
     734             :       ! b = -----------------
     735             :       !                2  7/2
     736             :       !     16 (D + H s  )
     737             :       !
     738           0 :       alpha = part1 - part2
     739           0 :       beta_s2 = 0.9375*SQRT_PI/D_Hs2_7_2
     740           0 :       Ebeta = E*beta_s2*s2
     741             : 
     742             :       ! the derivatives of a with respect to Fs^2 and Hs^2
     743             :       !
     744             :       !    da      d part1   d part2
     745             :       ! -------- = ------- - -------
     746             :       ! d (Hs^2)   d(Hs^2)   d(Hs^2)
     747             :       !
     748           0 :       dalpha_dh = dpart1_dh - dpart2_dh
     749             :       !
     750             :       !    2        2          2
     751             :       !   d a      d part1    d part2
     752             :       ! -------- = -------- - --------
     753             :       !        2          2          2
     754             :       ! d(Hs^2)    d(Hs^2)    d(Hs^2)
     755             :       !
     756           0 :       d2alpha_dh2 = d2part1_dh2 - d2part2_dh2
     757             :       !
     758             :       !    da      d part1
     759             :       ! -------- = -------
     760             :       ! d (Fs^2)   d(Fs^2)
     761             :       !
     762             :       !        2               2
     763             :       !       d a             d part1
     764             :       ! --------------- = ---------------
     765             :       ! d(Fs^2) d(Hs^2)   d(Fs^2) d(Hs^2)
     766             :       !
     767             :       ! (this expressions are already defined earlier)
     768             :       !
     769             :       ! and the derivatives of b/s^2 with respect to Hs^2
     770             :       !                        __
     771             :       ! d(b/s^2)          105 VPi          7  b / s^2
     772             :       ! -------- = - ----------------- = - -  --------
     773             :       ! d (Hs^2)                2  9/2     2         2
     774             :       !              32 (D + H s  )           D + H s
     775             :       !
     776           0 :       dbeta_dh = -3.5*beta_s2/D_Hs2
     777             :       !
     778             :       !  2                     __
     779             :       ! d (b/s^2)         945 VPi          9  d(b/s^2)
     780             :       ! --------- = ------------------ = - -  --------
     781             :       !         2              2  11/2     2  d (Hs^2)
     782             :       ! d (Hs^2)    64 (D + H s  )
     783             :       !
     784           0 :       d2beta_dh2 = -4.5*dbeta_dh/D_Hs2
     785             : 
     786             :       ! calculate the function G(s) and its derivatives
     787             :       !
     788             :       !       3Pi/4 + a
     789             :       ! G = - ---------
     790             :       !          E b
     791             :       !
     792           0 :       r3Pi_4_alpha = PI_75 + alpha
     793           0 :       Ebeta_s2 = E*beta_s2
     794           0 :       G = -r3Pi_4_alpha/Ebeta
     795             :       !
     796             :       !             /  /3 Pi    \  d(b/s^2)   /           d a    \  d(Hs^2)       d a   d(Fs^2)
     797             :       !            <  ( ---- + a ) -------   /  b/s^2  - -------  > -------  -  ------- -------
     798             :       ! d (Gs^2)    \  \ 4      /  d(Hs^2)  /            d(Hs^2) /    d s       d(Fs^2)   ds
     799             :       ! -------- = ----------------------------------------------------------------------------
     800             :       !    ds                                     E b/s^2
     801             :       !
     802           0 :       dGs2_ds = ((r3Pi_4_alpha*dbeta_dh/beta_s2 - dalpha_dh)*dHs2_ds - dalpha_df*dFs2_ds)/Ebeta_s2
     803             :       !
     804             :       !                                                     /3Pi    \  /             2\           /           \
     805             :       !                                                    ( --- + a )( b  h" + b   h' ) + 2b  h'( f'a  + h'a  ) + ...
     806             :       !                                            2        \ 4     /  \ h       hh   /      h    \   f      h/
     807             :       !  2         -f" a  - h" a  - 2 f' h' a   - h' a   + -----------------------------------------------------------
     808             :       ! d (Gs^2)        f       h            fh       hh                                b/s^2
     809             :       ! -------- = ---------------------------------------------------------------------------------------------------
     810             :       !   ds^2                                               E b/s^2
     811             :       !
     812             :       !            /3Pi    \   2  2
     813             :       !           ( --- + a ) b  h'
     814             :       !            \ 4     /   h
     815             :       ! ... = - 2 -----------------
     816             :       !                  b/s^2
     817             :       !
     818             :       ! where the primes indicate derivative with respect to s and the subscripts derivatives with respect to the subscript
     819             :       ! and f and h are abbrevations for f = Fs^2 and h = Hs^2
     820             :       !
     821             :       d2Gs2_ds2 = (-d2Fs2_ds2*dalpha_df - d2Hs2_ds2*dalpha_dh &
     822             :                    - 2.0*dFs2_ds*dHs2_ds*d2alpha_dfdh - dHs2_ds**2*d2alpha_dh2 &
     823             :                    + (r3Pi_4_alpha*(dbeta_dh*d2Hs2_ds2 + d2beta_dh2*dHs2_ds**2) &
     824             :                       + 2.0*dbeta_dh*dHs2_ds*(dFs2_ds*dalpha_df + dHs2_ds*dalpha_dh) &
     825           0 :                       - 2.0*r3Pi_4_alpha*(dbeta_dh*dHs2_ds)**2/beta_s2)/beta_s2)/Ebeta_s2
     826           0 :    END SUBROUTINE calculateG
     827             : 
     828             :    ! Calculate the function H(s) given in [3]
     829             :    ! Input:  s         - reduced gradient
     830             :    ! Output: H         - value of the function H
     831             :    !         dHs2_ds   - first derivative d(s^2*H(s))/ds
     832             :    !         d2Hs2_ds2 - second derivative d^2(s^2H(s))/ds^2
     833           0 :    SUBROUTINE calculateH(s, H, dHs2_ds, d2Hs2_ds2)
     834             :       IMPLICIT NONE
     835             : 
     836             :       REAL, INTENT(IN)  :: s
     837             :       REAL, INTENT(OUT) :: H, dHs2_ds, d2Hs2_ds2
     838             : 
     839             :       ! helper variables
     840             :       REAL :: s2, s3, s4, s5, s6
     841             :       REAL :: numer, denom
     842             :       REAL :: dnum_ds, dden_ds
     843             :       REAL :: d2num_ds2, d2den_ds2
     844             : 
     845             :       ! parameters given in [3]
     846             :       REAL, PARAMETER :: &
     847             :          a1 = 0.00979681, &
     848             :          a2 = 0.0410834, &
     849             :          a3 = 0.187440, &
     850             :          a4 = 0.00120824, &
     851             :          a5 = 0.0347188
     852             : 
     853             :       ! calculate helper variables
     854           0 :       s2 = s*s
     855           0 :       s3 = s2*s
     856           0 :       s4 = s2*s2
     857           0 :       s5 = s3*s2
     858           0 :       s6 = s4*s2
     859             : 
     860             :       ! calculate function H(s) with [3] eq. (A5)
     861             :       !
     862             :       !            2         4
     863             :       !        a  s   +  a  s
     864             :       !         1         2
     865             :       ! H = -------------------------
     866             :       !             4       5       6
     867             :       !     1 + a  s  + a  s  + a  s
     868             :       !          3       4       5
     869             :       !
     870           0 :       numer = a1*s2 + a2*s4
     871           0 :       denom = 1.0 + a3*s4 + a4*s5 + a5*s6
     872           0 :       H = numer/denom
     873             : 
     874             :       ! calculate the first derivative of s^2 H(s)
     875             :       !
     876             :       !     /    \
     877             :       !  d | f(x) |   f'(x) - f(x)g'(x) / g(x)
     878             :       ! -- | ---- | = ------------------------
     879             :       ! dx | g(x) |            g(x)
     880             :       !     \    /
     881             :       !
     882           0 :       numer = numer*s2
     883           0 :       dnum_ds = 4.0*s3*a1 + 6.0*s5*a2
     884           0 :       dden_ds = 4.0*s3*a3 + 5.0*s4*a4 + 6.0*s5*a5
     885           0 :       dHs2_ds = (dnum_ds - (dden_ds*numer)/denom)/denom
     886             : 
     887             :       ! calculate the second derivative of s^2 H(s)
     888             :       !
     889             :       !                                                               2
     890             :       !                         2 f'(x)g'(x) + f(x)g"(x) - 2 f(x)g'(x) /g(x)
     891             :       !    2  /    \    f"(x) - --------------------------------------------
     892             :       !   d  | f(x) |                              g(x)
     893             :       ! ---- | ---- | = ----------------------------------------------------
     894             :       ! dx^2 | g(x) |                          g(x)
     895             :       !       \    /
     896             :       !
     897           0 :       d2num_ds2 = 12.0*s2*a1 + 30.0*s4*a2
     898           0 :       d2den_ds2 = 12.0*s2*a3 + 20.0*s3*a4 + 30.0*s4*a5
     899           0 :       d2Hs2_ds2 = (d2num_ds2 - (2.0*dnum_ds*dden_ds + numer*d2den_ds2 - 2.0*numer*dden_ds**2/denom)/denom)/denom
     900           0 :    END SUBROUTINE calculateH
     901             : 
     902             :    ! Calculate several Fourier transformations
     903             :    ! a) the Fourier transformed potential
     904             :    !                                                 /          \
     905             :    !          /      | erf |       \     4 Pi       |   |k+G|^2  |
     906             :    ! V (k) = ( k + G | --- | k + G  ) = -------  Exp| - -------  |                             [1]
     907             :    !  G       \      |  r  |       /    |k+G|^2     |    4 w^2   |
     908             :    !                                                 \          /
     909             :    !
     910             :    ! b) muffin-tin basis function
     911             :    !                                                            R
     912             :    !                                     -L                      /
     913             :    !           /      |      \     4 Pi i         ^     -i G R  |     2
     914             :    ! MT (k) = ( k + G | M     ) = ---------- Y ( k+G ) e        | dr r  Phi(r) j ( |k+G| r )  [2]
     915             :    !   G,I     \      |  k,I /     Sqrt(Om)   LM                |               L
     916             :    !                                                           /
     917             :    !                                                            0
     918             :    !
     919             :    ! c) interstitial basis function
     920             :    !                                                      -----
     921             :    !          /      |      \                       4 Pi   \     -i(G - G_I)R_a  Sin(|G_I - G| R_a) - |G_I - G| R_a Cos(|G_I - G| R_a)
     922             :    ! IN    = ( k + G | M     ) = kronecker(G,G ) - ------   )   e               ------------------------------------------------------- [3]
     923             :    !   G,I    \      |  k,I /                 I      Om    /                                          |G_I - G|^3
     924             :    !                                                      -----
     925             :    !                                                        a
     926             :    !                                                                            \_________________________ ___________________________/
     927             :    !                                                                                                      V
     928             :    !                                                                                                     I_a
     929             :    !
     930             :    ! In the code:
     931             :    ! V_G(k):  potential
     932             :    ! MT_G(k): muffintin
     933             :    ! IN_G:    interstitial
     934             :    ! I_a:     inter_atom
     935             :    !
     936             :    ! Input:
     937             :    ! rmsh        - array of radial mesh points
     938             :    ! rmt         - radius of the muffin tin
     939             :    ! dx          - logarithmic increment of radial mesh
     940             :    ! jri         - number of radial mesh points
     941             :    ! jmtd        - dimension of radial mesh points
     942             :    ! bk          - one(!) k-vector for which FT is calculated
     943             :    ! bmat        - reciprocal lattice vector for all directions
     944             :    ! vol         - volume of the unit cell
     945             :    ! ntype       - number of atom types
     946             :    ! neq         - number of symmetry-equivalent atoms of atom type i
     947             :    ! natd        - dimesion of atoms in unit cell
     948             :    ! taual       - vector of atom positions (internal coordinates)
     949             :    ! lcutm       - l cutoff for mixed basis
     950             :    ! maxlcutm    - maximum of all these l cutoffs
     951             :    ! nindxm      - number of radial functions of mixed basis
     952             :    ! maxindxm    - maximum of these numbers
     953             :    ! gptm        - reciprocal lattice vectors of the mixed basis (internal coord.)
     954             :    ! ngptm       - number of vectors (for treated k-vector)
     955             :    ! pgptm       - pointer to the appropriate g-vector (for treated k-vector)
     956             :    ! gptmd       - dimension of gptm
     957             :    ! basm        - mixed basis functions (mt + inter) for treated k-vector
     958             :    ! lexp        - cutoff of spherical harmonics expansion of plane wave
     959             :    ! noGPts      - no g-vectors used for Fourier trafo
     960             :    ! Output:
     961             :    ! potential   - Fourier transformation of the potential
     962             :    ! muffintin   - Fourier transformation of all MT functions
     963             :    ! interstital - Fourier transformation of all IR functions
     964           0 :    SUBROUTINE calculate_fourier_transform( &
     965             :       ! Input
     966           0 :       rmsh, rmt, dx, jri, jmtd, bk, &
     967           0 :       bmat, vol, ntype, neq, natd, taual, lcutm, maxlcutm, &
     968           0 :       nindxm, maxindxm, gptm, ngptm, pgptm, gptmd, &
     969           0 :       basm, noGPts, irank, &
     970             :       ! Output
     971           0 :       potential, muffintin, interstitial)
     972             : 
     973             :       USE m_constants
     974             :       USE m_olap, ONLY: gptnorm
     975             :       USE m_util, ONLY: sphbessel, pure_intgrf, intgrf_init, intgrf_out, NEGATIVE_EXPONENT_WARNING, NEGATIVE_EXPONENT_ERROR
     976             : 
     977             :       IMPLICIT NONE
     978             : 
     979             :       ! scalar input
     980             :       INTEGER, INTENT(IN)    :: natd, ntype, maxlcutm
     981             :       INTEGER, INTENT(IN)    :: jmtd, irank
     982             :       INTEGER, INTENT(IN)    :: maxindxm
     983             :       INTEGER, INTENT(IN)    :: gptmd, noGPts
     984             :       REAL, INTENT(IN)       :: vol
     985             : 
     986             :       ! array input
     987             :       INTEGER, INTENT(IN)    :: lcutm(ntype)
     988             :       INTEGER, INTENT(IN)    :: nindxm(0:maxlcutm, ntype), neq(ntype)
     989             :       INTEGER, INTENT(IN)    :: jri(ntype)
     990             :       INTEGER, INTENT(IN)    :: gptm(3, gptmd)
     991             :       INTEGER, INTENT(IN)    :: ngptm
     992             :       INTEGER, INTENT(IN)    :: pgptm(ngptm)
     993             : 
     994             :       REAL, INTENT(IN)       :: bk(3)
     995             :       REAL, INTENT(IN)       :: rmsh(jmtd, ntype), rmt(ntype), dx(ntype)
     996             :       REAL, INTENT(IN)       :: basm(jmtd, maxindxm, 0:maxlcutm, ntype)
     997             :       REAL, INTENT(IN)       :: bmat(3, 3)!,amat(3,3)
     998             :       REAL, INTENT(IN)       :: taual(3, natd)
     999             : 
    1000             :       ! array output
    1001             :       REAL, INTENT(OUT)   :: potential(noGPts)                           ! Fourier transformed potential
    1002             :       COMPLEX, INTENT(OUT)   :: muffintin(noGPts, maxindxm, &                 ! muffin-tin overlap integral
    1003             :                                           (maxlcutm + 1)**2, ntype, MAXVAL(neq))
    1004             :       COMPLEX, INTENT(OUT)   :: interstitial(noGPts, gptmd)                  ! interstistial overlap intergral
    1005             : 
    1006             :       ! private scalars
    1007             :       INTEGER                :: cg, cg2, ci, cl, cn, cr                          ! counter variables
    1008             :       REAL                   ::  r2Pi, r4Pi, pi_omega2                   !  2*Pi, 4*Pi, Pi/omega^2
    1009             :       REAL                   :: sVol, r4Pi_sVol, r4Pi_Vol                   ! sqrt(vol), 4*Pi/sqrt(Vol), 4*Pi/Vol
    1010             :       REAL                   :: omega, r1_omega2, r1_4omega2
    1011             : !     REAL,    PARAMETER     :: omega = omega_VHSE()                        ! omega of the HSE functional
    1012             : !     REAL,    PARAMETER     :: r1_omega2  = 1.0 / omega**2                 ! 1/omega^2
    1013             : !     REAL,    PARAMETER     :: r1_4omega2 = 0.25 * r1_omega2               ! 1/(4 omega^2)
    1014             :       COMPLEX, PARAMETER     :: img = (0.0, 1.0)                             ! i
    1015             : 
    1016             :       ! private arrays
    1017           0 :       INTEGER                :: gPts(3, noGPts)                              ! g vectors (internal units)
    1018           0 :       INTEGER                :: gPts_gptm(3, noGpts, gptmd)                   ! gPts - gptm
    1019           0 :       INTEGER                :: natdPtr(ntype + 1)                            ! pointer to all atoms of one type
    1020           0 :       REAL, ALLOCATABLE   :: gridf(:, :)                                  ! grid for radial integration
    1021           0 :       REAL                   :: k_G(3, noGPts)                               ! k + G
    1022           0 :       REAL                   :: AbsK_G(noGPts), AbsK_G2(noGPts)              ! abs(k+G), abs(k+G)^2
    1023           0 :       REAL                   :: arg(noGPts)                                 ! abs(k+G)^2 / (4*omega^2)
    1024           0 :       REAL                   :: sphbesK_Gr(noGPts, jmtd, 0:maxlcutm, ntype)    ! spherical bessel function of abs(k+G)r
    1025           0 :       TYPE(intgrf_out)       :: intgrMT(noGPts, maxindxm, 0:maxlcutm, ntype)   ! integration in muffin-tin
    1026           0 :       REAL                   :: abs_dg(noGpts, gptmd)                        ! abs(gPts - gptm)
    1027           0 :       COMPLEX                :: imgl(0:maxlcutm)                            ! i^l
    1028           0 :       COMPLEX                :: Ylm(noGPts, (maxlcutm + 1)**2)                 ! spherical harmonics for k+G and all lm
    1029           0 :       COMPLEX                :: expIGR(noGPts, ntype, MAXVAL(neq))            ! exp(-iGR) for all atom types
    1030           0 :       COMPLEX                :: sumInter_atom(noGpts, gptmd)                 ! sum over inter-atom factors
    1031             : 
    1032             :       ! Calculate helper variables
    1033           0 :       r2Pi = 2.0*pi_const
    1034           0 :       r4Pi = 4.0*pi_const
    1035           0 :       sVol = SQRT(vol)
    1036           0 :       r4Pi_sVol = r4Pi/sVol
    1037           0 :       r4Pi_Vol = r4Pi/Vol
    1038           0 :       omega = omega_hse!omega_VHSE()
    1039           0 :       r1_omega2 = 1.0/omega**2
    1040           0 :       r1_4omega2 = 0.25*r1_omega2
    1041           0 :       pi_omega2 = pi_const*r1_omega2
    1042             : 
    1043             :       ! calculate pointers for all atom-types to all atoms
    1044           0 :       natdPtr(ntype + 1) = natd
    1045           0 :       DO cn = ntype, 1, -1
    1046           0 :          natdPtr(cn) = natdPtr(cn + 1) - neq(cn)
    1047             :       END DO
    1048             : 
    1049             :       ! Define imgl(l) = img**l
    1050           0 :       imgl(0) = 1.0
    1051           0 :       DO ci = 1, maxlcutm
    1052           0 :          imgl(ci) = imgl(ci - 1)*img
    1053             :       END DO
    1054             : 
    1055             :       ! generate grid for fast radial integration
    1056             :       CALL intgrf_init(ntype, jmtd, jri, dx, rmsh, gridf)
    1057             : 
    1058             :       ! Set all arrays to 0
    1059           0 :       gPts = 0
    1060           0 :       k_G = 0
    1061           0 :       AbsK_G = 0
    1062           0 :       arg = 0
    1063           0 :       sphbesK_Gr = 0
    1064           0 :       intgrMT%value = 0
    1065           0 :       intgrMT%ierror = 0
    1066           0 :       Ylm = 0
    1067           0 :       expIGR = 0
    1068           0 :       gPts_gptm = 0
    1069           0 :       abs_dg = 0
    1070           0 :       sumInter_atom = 0
    1071           0 :       potential = 0
    1072           0 :       muffintin = 0
    1073           0 :       interstitial = 0
    1074             : 
    1075             : ! Calculate the muffin-tin basis function overlap integral
    1076             : !                                                            R
    1077             : !                                     -L                      /
    1078             : !           /      |      \     4 Pi i         ^     -i G R  |     2
    1079             : ! MT (k) = ( k + G | M     ) = ---------- Y ( k+G ) e        | dr r  Phi(r) j ( |k+G| r )  [2]
    1080             : !   G,I     \      |  k,I /     Sqrt(Om)   LM                |               L
    1081             : !                                                           /
    1082             : !                                                            0
    1083           0 :       IF (ngptm < noGpts) STOP 'hsefunctional: error calculating Fourier coefficients, noGpts too large'
    1084             : 
    1085           0 :       gPts(:, :) = gptm(:, pgptm(1:noGPts))
    1086             : #ifndef __PGI
    1087             : 
    1088           0 :       gpoints:FORALL (cg=1:noGPts)
    1089             :       ntypesA:FORALL (cn=1:ntype)
    1090             :       ! Calculate the phase factor exp(-iGR) for all atoms
    1091             :       FORALL (ci=1:neq(cn))
    1092             :       expIGR(cg, cn, ci) = EXP(-r2Pi*img*DOT_PRODUCT(taual(:, natdPtr(cn) + ci), gPts(:, cg)))
    1093             :       muffintin(cg, :, :, cn, ci) = r4Pi_sVol*expIGR(cg, cn, ci)
    1094             :       END FORALL
    1095             : 
    1096             :       ! Multiplication with i^-L
    1097             :       FORALL (cl=0:lcutm(cn))
    1098             :       muffintin(cg, :, cl*cl + 1:(cl + 1)*(cl + 1), cn, :) = muffintin(cg, :, cl*cl + 1:(cl + 1)*(cl + 1), cn, :)/imgl(cl)
    1099             :       END FORALL
    1100             : 
    1101             :       END FORALL ntypesA
    1102             : 
    1103             :       ! Calculate the k+G, abs(k+G), and abs(k+G)^2
    1104             :       k_G(:, cg) = MATMUL(gPts(:, cg) + bk(:), bmat(:, :))
    1105             :       AbsK_G2(cg) = SUM(k_G(:, cg)**2)
    1106             :       AbsK_G(cg) = SQRT(AbsK_G2(cg))
    1107             : 
    1108             :       ! Spherical harmonics are calculated for all lm's and k+G's
    1109             :       Ylm(cg, :) = calcYlm(k_G(:, cg), maxlcutm)
    1110             : 
    1111             :       ! Perform the integration in eq.[2] for all muffin-tins
    1112             :       ntypesB:FORALL (cn=1:ntype)
    1113             : 
    1114             :       ! Multiplication with the spherical harmonic
    1115             :       FORALL (cl=1:(lcutm(cn) + 1)**2)
    1116             :       muffintin(cg, :, cl, cn, :) = muffintin(cg, :, cl, cn, :)*Ylm(cg, cl)
    1117             :       END FORALL
    1118             : 
    1119             :       ! Calculate the spherical bessel function
    1120             :       FORALL (cr=1:jri(cn))
    1121             :       sphbesK_Gr(cg, cr, :, cn) = calcSphBes(AbsK_G(cg)*rmsh(cr, cn), maxlcutm)
    1122             :       END FORALL
    1123             :       ! integrate the function and multiplicate to the result
    1124             :       FORALL (cl=0:lcutm(cn))
    1125             :       FORALL (ci=1:nindxm(cl, cn))
    1126             :       intgrMT(cg, ci, cl, cn) = pure_intgrf(rmsh(:, cn)*basm(:, ci, cl, cn)* &
    1127             :                                             sphbesK_Gr(cg, :, cl, cn), jri, jmtd, rmsh, dx, ntype, cn, gridf)
    1128             :       muffintin(cg, ci, cl*cl + 1:(cl + 1)*(cl + 1), cn, :) = &
    1129             :          muffintin(cg, ci, cl*cl + 1:(cl + 1)*(cl + 1), cn, :)*intgrMT(cg, ci, cl, cn)%value
    1130             :       END FORALL
    1131             :       END FORALL
    1132             : 
    1133             :       END FORALL ntypesB
    1134             : 
    1135             : ! calculate the overlap with the interstitial basis function
    1136             : !                                                      -----
    1137             : !          /      |      \                       4 Pi   \     -i(G - G_I)R_a  Sin(|G_I - G| R_a) - |G_I - G| R_a Cos(|G_I - G| R_a)
    1138             : ! IN    = ( k + G | M     ) = kronecker(G,G ) - ------   )   e               ------------------------------------------------------- [3]
    1139             : !   G,I    \      |  k,I /                 I      Om    /                                          |G_I - G|^3
    1140             : !                                                      -----
    1141             : !                                                        a
    1142             : !                                                                            \_________________________ ___________________________/
    1143             : !                                                                                                      V
    1144             : !                                                                                                     I_a
    1145             :       ! Calculate the difference of the G vectors and its absolute value
    1146             :       FORALL (cg2=1:gptmd)
    1147             : 
    1148             :       gPts_gptm(:, cg, cg2) = gPts(:, cg) - gptm(:, cg2)
    1149             :       abs_dg(cg, cg2) = gptnorm(gPts_gptm(:, cg, cg2), bmat)
    1150             : 
    1151             :       END FORALL
    1152             : 
    1153             :       END FORALL gpoints
    1154             : 
    1155             :       ! Check if any of the integrations failed and abort if one did
    1156           0 :       IF (ANY(intgrMT%ierror == NEGATIVE_EXPONENT_ERROR)) THEN
    1157           0 :          IF (irank == 0) WRITE (6, *) 'intgrf: Warning! Negative exponent x in extrapolation a+c*r**x'
    1158           0 :       ELSEIF (ANY(intgrMT%ierror == NEGATIVE_EXPONENT_WARNING)) THEN
    1159           0 :          IF (irank == 0) WRITE (6, *) 'intgrf: Negative exponent x in extrapolation a+c*r**x'
    1160           0 :          STOP 'intgrf: Negative exponent x in extrapolation a+c*r**x'
    1161             :       END IF
    1162             : 
    1163             :       ! Calculate the interstitial value with eq.[3] using the limit
    1164             :       !           3
    1165             :       !        R_a
    1166             :       ! I_a = ------
    1167             :       !         3
    1168             :       ! if there's no difference of the two G vectors
    1169           0 :       WHERE (abs_dg == 0)
    1170             :       sumInter_atom = DOT_PRODUCT(rmt**3, neq)/3.0
    1171             :       interstitial = 1.0 - r4Pi_Vol*sumInter_atom
    1172             :       ELSEWHERE
    1173             :       sumInter_atom = calculateSummation(abs_dg, gPts_gptm)
    1174             :       interstitial = -r4Pi_Vol*sumInter_atom
    1175             :       END WHERE
    1176             : 
    1177             : ! Calculate the Fourier transformed potential
    1178             : !                                                 /          \
    1179             : !          /      | erf |       \      4 Pi      |   |k+G|^2  |
    1180             : ! V (k) = ( k + G | --- | k + G' ) = -------  Exp| - -------  | kronecker(G,G')             [1]
    1181             : !  G       \      |  r  |       /    |k+G|^2     |    4 w^2   |
    1182             : !                                                 \          /
    1183           0 :       WHERE (absK_G /= 0)
    1184             :       ! Calculate (k+G)^2 / 4*w^2
    1185             :       arg = -r1_4omega2*AbsK_G2
    1186             : 
    1187             :       ! The potential is calculated using eq.[1]
    1188             :       potential = r4Pi/AbsK_G2*EXP(arg)
    1189             : 
    1190             : ! Calculate the Fourier transformed potential for the full(!) potential
    1191             : !
    1192             : !          /          | erfc |          \      Pi
    1193             : ! V (k) = ( k + G = 0 | ---- | k + G = 0 ) = -----
    1194             : !  G       \          |   r  |          /     w^2
    1195             : !
    1196             :       ELSEWHERE
    1197             :       ! the negative sign is added to compensate the sign when the
    1198             :       ! contribution is added to the Coulomb matrix
    1199             :       potential = -pi_omega2
    1200             :       endwhere
    1201             : 
    1202           0 :       DEALLOCATE (gridf)
    1203             : #else
    1204             :       call judft_error("hsefunctional not implemented for PGI")
    1205             : #endif
    1206             :    CONTAINS
    1207             : 
    1208             :       ! Calculates the inter-atom parts and the summation over all atoms:
    1209             :       ! -----
    1210             :       !  \     -i(G - G_I)R_a  Sin(|G_I - G| R_a) - |G_I - G| R_a Cos(|G_I - G| R_a)
    1211             :       !   )   e               -------------------------------------------------------
    1212             :       !  /                                          |G_I - G|^3
    1213             :       ! -----
    1214             :       !   a
    1215             :       ! Input:  abs_dg    - absolute value |G_I - G|
    1216             :       !         gPts_gptm - vector G - G_I
    1217           0 :       PURE FUNCTION calculateSummation(abs_dg, gPts_gptm)
    1218             :          IMPLICIT NONE
    1219             :          INTEGER, INTENT(IN) :: gPts_gptm(3, noGPts, gptmd)
    1220             :          REAL, INTENT(IN)    :: abs_dg(noGPts, gptmd)
    1221             :          COMPLEX             :: calculateSummation(noGPts, gptmd)
    1222             :          INTEGER             :: cn, ci                                     ! counter variables
    1223           0 :          REAL                :: abs_dgR(noGPts, gptmd, ntype)               ! abs(gPts - gptm)*R (R: radius MT)
    1224           0 :          REAL                :: inter_atom(noGPts, gptmd, ntype)            ! inter-atom interaction for interstitial
    1225           0 :          COMPLEX             :: expIdGR(noGPts, gptmd, ntype, MAXVAL(neq))   ! exp(-i(gPts-gptm)R)
    1226           0 :          COMPLEX             :: sumExpIdGR(noGPts, gptmd, ntype)            ! sum over atom of same type
    1227             : 
    1228           0 :          atoms:FORALL (cn=1:ntype)
    1229             :          !                                                         -i(G-G_I)R_a
    1230             :          ! Calculate for all similar atoms (same type) the factor e
    1231             :          ! use that the G vectors and atomic positions are given in internal units
    1232             :          FORALL (ci=1:neq(cn))
    1233             :          expIdGR(:, :, cn, ci) = EXP(-r2pi*img*my_dot_product(gPts_gptm(:, :, :), taual(:, natdPtr(cn) + ci)))
    1234             :          END FORALL
    1235             :          ! Add all factors for atoms of the same type
    1236             :          sumExpIdGR(:, :, cn) = my_sum(expIdGR(:, :, cn, :))
    1237             : 
    1238             :          ! Calculate the inter-atom factor which is the same for all atoms of the same type
    1239             :          abs_dgR(:, :, cn) = abs_dg*rmt(cn)
    1240             :          inter_atom(:, :, cn) = (SIN(abs_dgR(:, :, cn)) - abs_dgR(:, :, cn)*COS(abs_dgR(:, :, cn)))/abs_dg**3
    1241             :          END FORALL atoms
    1242             : 
    1243             :          ! Add the factors of all atoms together
    1244           0 :          calculateSummation = my_dot_product(sumExpIdGR, inter_atom)
    1245             :       END FUNCTION calculateSummation
    1246             : 
    1247             :    END SUBROUTINE calculate_fourier_transform
    1248             : 
    1249             :    ! Calculate several Fourier transformations
    1250             :    ! a) the Fourier transformed potential
    1251             :    !                                                 /          \
    1252             :    !          /      | erf |       \     4 Pi       |   |k+G|^2  |
    1253             :    ! V (k) = ( k + G | --- | k + G  ) = -------  Exp| - -------  |                             [1]
    1254             :    !  G       \      |  r  |       /    |k+G|^2     |    4 w^2   |
    1255             :    !                                                 \          /
    1256             :    !
    1257             :    ! b) muffin-tin basis function
    1258             :    !                                                            R
    1259             :    !                                     -L                      /
    1260             :    !           /      |      \     4 Pi i         ^     -i G R  |     2
    1261             :    ! MT (k) = ( k + G | M     ) = ---------- Y ( k+G ) e        | dr r  Phi(r) j ( |k+G| r )  [2]
    1262             :    !   G,I     \      |  k,I /     Sqrt(Om)   LM                |               L
    1263             :    !                                                           /
    1264             :    !                                                            0
    1265             :    !
    1266             :    ! In the code:
    1267             :    ! V_G(k):  potential
    1268             :    ! MT_G(k): muffintin
    1269             :    !
    1270             :    ! Input:
    1271             :    ! rmsh        - array of radial mesh points
    1272             :    ! rmt         - radius of the muffin tin
    1273             :    ! dx          - logarithmic increment of radial mesh
    1274             :    ! jri         - number of radial mesh points
    1275             :    ! jmtd        - dimension of radial mesh points
    1276             :    ! bk          - one(!) k-vector for which FT is calculated
    1277             :    ! bmat        - reciprocal lattice vector for all directions
    1278             :    ! vol         - volume of the unit cell
    1279             :    ! ntype       - number of atom types
    1280             :    ! neq         - number of symmetry-equivalent atoms of atom type i
    1281             :    ! natd        - dimesion of atoms in unit cell
    1282             :    ! taual       - vector of atom positions (internal coordinates)
    1283             :    ! lcutm       - l cutoff for mixed basis
    1284             :    ! maxlcutm    - maximum of all these l cutoffs
    1285             :    ! nindxm      - number of radial functions of mixed basis
    1286             :    ! maxindxm    - maximum of these numbers
    1287             :    ! gptm        - reciprocal lattice vectors of the mixed basis (internal coord.)
    1288             :    ! ngptm       - number of vectors (for treated k-vector)
    1289             :    ! pgptm       - pointer to the appropriate g-vector (for treated k-vector)
    1290             :    ! gptmd       - dimension of gptm
    1291             :    ! basm        - mixed basis functions (mt + inter) for treated k-vector
    1292             :    ! lexp        - cutoff of spherical harmonics expansion of plane wave
    1293             :    ! noGPts      - no g-vectors used for Fourier trafo
    1294             :    ! Output:
    1295             :    ! potential   - Fourier transformation of the potential
    1296             :    ! muffintin   - Fourier transformation of all MT functions
    1297           0 :    SUBROUTINE calculate_fourier_transform_once( &
    1298             :       ! Input
    1299           0 :       rmsh, rmt, dx, jri, jmtd, bk, ikpt, nkptf, &
    1300           0 :       bmat, vol, ntype, neq, natd, taual, lcutm, maxlcutm, &
    1301           0 :       nindxm, maxindxm, gptm, ngptm, pgptm, gptmd, &
    1302           0 :       nbasp, basm, noGPts, invsat, invsatnr, irank, &
    1303             :       ! Output
    1304           0 :       potential, fourier_trafo)
    1305             : 
    1306             :       USE m_constants
    1307             :       USE m_util, ONLY: sphbessel, pure_intgrf, intgrf_init, intgrf_out, NEGATIVE_EXPONENT_WARNING, NEGATIVE_EXPONENT_ERROR
    1308             :       USE m_trafo, ONLY: symmetrize
    1309             : 
    1310             :       IMPLICIT NONE
    1311             : 
    1312             :       ! scalar input
    1313             :       INTEGER, INTENT(IN)    :: natd, ntype, maxlcutm
    1314             :       INTEGER, INTENT(IN)    :: ikpt, nkptf, jmtd
    1315             :       INTEGER, INTENT(IN)    :: maxindxm, irank
    1316             :       INTEGER, INTENT(IN)    :: gptmd, noGPts
    1317             :       REAL, INTENT(IN)       :: vol
    1318             : 
    1319             :       ! array input
    1320             :       INTEGER, INTENT(IN)    :: lcutm(ntype)
    1321             :       INTEGER, INTENT(IN)    :: nindxm(0:maxlcutm, ntype), neq(ntype)
    1322             :       INTEGER, INTENT(IN)    :: jri(ntype)
    1323             :       INTEGER, INTENT(IN)    :: gptm(3, gptmd)
    1324             :       INTEGER, INTENT(IN)    :: ngptm, nbasp
    1325             :       INTEGER, INTENT(IN)    :: pgptm(ngptm)
    1326             :       INTEGER, INTENT(IN)    :: invsat(natd), invsatnr(natd)
    1327             : 
    1328             :       REAL, INTENT(IN)       :: bk(3)
    1329             :       REAL, INTENT(IN)       :: rmsh(jmtd, ntype), rmt(ntype), dx(ntype)
    1330             :       REAL, INTENT(IN)       :: basm(jmtd, maxindxm, 0:maxlcutm, ntype)
    1331             :       REAL, INTENT(IN)       :: bmat(3, 3)
    1332             :       REAL, INTENT(IN)       :: taual(3, natd)
    1333             : 
    1334             :       ! array output
    1335             :       REAL, INTENT(OUT)   :: potential(noGPts)                           ! Fourier transformed potential
    1336             : #ifdef CPP_INVERSION
    1337             :       REAL, INTENT(OUT)   :: fourier_trafo(nbasp, noGPts) !muffintin_out(nbasp,noGPts)
    1338             : #else
    1339             :       COMPLEX, INTENT(OUT)   :: fourier_trafo(nbasp, noGPts) !muffintin_out(nbasp,noGPts)
    1340             : #endif
    1341             : 
    1342             :       ! private scalars
    1343             :       INTEGER                :: cg, cg2, ci, cl, cm, cn, cr                       ! counter variables
    1344             :       REAL                   :: r2Pi, r4Pi, pi_omega2                   ! Pi, 2*Pi, 4*Pi, Pi/omega^2
    1345             :       REAL                   :: sVol, r4Pi_sVol, r4Pi_Vol                   ! sqrt(vol), 4*Pi/sqrt(Vol), 4*Pi/Vol
    1346             :       REAL                   :: omega, r1_omega2, r1_4omega2
    1347             : !     REAL,    PARAMETER     :: omega = omega_VHSE()                        ! omega of the HSE functional
    1348             : !     REAL,    PARAMETER     :: r1_omega2  = 1.0 / omega**2                 ! 1/omega^2
    1349             : !     REAL,    PARAMETER     :: r1_4omega2 = 0.25 * r1_omega2               ! 1/(4 omega^2)
    1350             :       COMPLEX, PARAMETER     :: img = (0.0, 1.0)                             ! i
    1351             : 
    1352             :       ! private arrays
    1353           0 :       INTEGER                :: gPts(3, noGPts)                              ! g vectors (internal units)
    1354           0 :       INTEGER                :: natdPtr(ntype + 1)                            ! pointer to all atoms of one type
    1355           0 :       INTEGER                :: ptrType(nbasp), ptrEq(nbasp), &               ! pointer from ibasp to corresponding
    1356           0 :                                 ptrLM(nbasp), ptrN(nbasp)                    ! type, atom, l and m, and n
    1357           0 :       REAL, ALLOCATABLE   :: gridf(:, :)                                  ! grid for radial integration
    1358           0 :       REAL                   :: k_G(3, noGPts)                               ! k + G
    1359           0 :       REAL                   :: AbsK_G(noGPts), AbsK_G2(noGPts)              ! abs(k+G), abs(k+G)^2
    1360           0 :       REAL                   :: arg(noGPts)                                 ! abs(k+G)^2 / (4*omega^2)
    1361           0 :       REAL                   :: sphbesK_Gr(noGPts, jmtd, 0:maxlcutm, ntype)    ! spherical bessel function of abs(k+G)r
    1362           0 :       TYPE(intgrf_out)       :: intgrMT(noGPts, maxindxm, 0:maxlcutm, ntype)   ! integration in muffin-tin
    1363           0 :       COMPLEX                :: imgl(0:maxlcutm)                            ! i^l
    1364           0 :       COMPLEX                :: Ylm(noGPts, (maxlcutm + 1)**2)                 ! spherical harmonics for k+G and all lm
    1365           0 :       COMPLEX                :: expIGR(noGPts, ntype, MAXVAL(neq))            ! exp(-iGR) for all atom types
    1366             :       COMPLEX                :: muffintin(noGPts, maxindxm, &                 ! muffin-tin overlap integral
    1367           0 :                                           (maxlcutm + 1)**2, ntype, MAXVAL(neq))
    1368             : #ifdef CPP_INVERSION
    1369             :       COMPLEX                :: sym_muffintin(nbasp, noGPts)                 ! symmetrized muffin tin
    1370             : #endif
    1371             : 
    1372             :       LOGICAL, SAVE          :: first_entry = .TRUE.                        ! allocate arrays in first entry
    1373             : 
    1374             :       ! allocate arrays in first entry reuse later
    1375           0 :       IF (first_entry) THEN
    1376           0 :          ALLOCATE (already_known(nkptf), &            ! stores which elements are known
    1377           0 :                    known_potential(maxNoGPts, nkptf), &            ! stores the potential for all k-points
    1378           0 :                    known_fourier_trafo(nbasp, maxNoGPts, nkptf))            ! stores the fourier transform of the mixed basis
    1379             :          ! initialization
    1380           0 :          already_known = .FALSE.
    1381           0 :          known_potential = 0.0
    1382           0 :          known_fourier_trafo = 0.0
    1383             :          ! unset flag as arrays are allocated
    1384           0 :          first_entry = .FALSE.
    1385             :       ELSE
    1386             :          ! check if size of arrays has changed and stop with error if they did
    1387           0 :          IF (SIZE(already_known) /= nkptf) STOP 'hsefunctional: Array size changed!'
    1388           0 :          IF (SIZE(known_fourier_trafo, 1) /= nbasp) STOP 'hsefunctional: Array size changed!'
    1389             :       END IF
    1390             : 
    1391             :       ! if the current k-point was not calculated yet
    1392           0 :       IF (.NOT. (already_known(ikpt))) THEN
    1393             : 
    1394             :          ! Calculate helper variables
    1395           0 :          r2Pi = 2.0*pi_const
    1396           0 :          r4Pi = 4.0*pi_const
    1397           0 :          sVol = SQRT(vol)
    1398           0 :          r4Pi_sVol = r4Pi/sVol
    1399           0 :          r4Pi_Vol = r4Pi/Vol
    1400           0 :          omega = omega_hse!omega_VHSE()
    1401           0 :          r1_omega2 = 1.0/omega**2
    1402           0 :          r1_4omega2 = 0.25*r1_omega2
    1403           0 :          pi_omega2 = pi_const*r1_omega2
    1404             : 
    1405             :          ! calculate pointers for all atom-types to all atoms
    1406           0 :          natdPtr(ntype + 1) = natd
    1407           0 :          DO cn = ntype, 1, -1
    1408           0 :             natdPtr(cn) = natdPtr(cn + 1) - neq(cn)
    1409             :          END DO
    1410             : 
    1411             :          ! Define imgl(l) = img**l
    1412           0 :          imgl(0) = 1.0
    1413           0 :          DO ci = 1, maxlcutm
    1414           0 :             imgl(ci) = imgl(ci - 1)*img
    1415             :          END DO
    1416             : 
    1417             :          ! generate grid for fast radial integration
    1418             :          CALL intgrf_init(ntype, jmtd, jri, dx, rmsh, gridf)
    1419             : 
    1420             :          ! Set all arrays to 0
    1421           0 :          gPts = 0
    1422           0 :          k_G = 0
    1423           0 :          AbsK_G = 0
    1424           0 :          arg = 0
    1425           0 :          sphbesK_Gr = 0
    1426           0 :          intgrMT%value = 0
    1427           0 :          intgrMT%ierror = 0
    1428           0 :          Ylm = 0
    1429           0 :          expIGR = 0
    1430           0 :          potential = 0
    1431           0 :          muffintin = 0
    1432             : 
    1433             : ! Calculate the muffin-tin basis function overlap integral
    1434             : !                                                            R
    1435             : !                                     -L                      /
    1436             : !           /      |      \     4 Pi i         ^     -i G R  |     2
    1437             : ! MT (k) = ( k + G | M     ) = ---------- Y ( k+G ) e        | dr r  Phi(r) j ( |k+G| r )  [2]
    1438             : !   G,I     \      |  k,I /     Sqrt(Om)   LM                |               L
    1439             : !                                                           /
    1440             : !                                                            0
    1441           0 :          IF (ngptm < noGpts) STOP 'hsefunctional: error calculating Fourier coefficients, noGpts too large'
    1442             : 
    1443           0 :          gPts(:, :) = gptm(:, pgptm(1:noGPts))
    1444             : 
    1445           0 :          gpoints:FORALL (cg=1:noGPts)
    1446             :          ntypesA:FORALL (cn=1:ntype)
    1447             : 
    1448             :          ! Calculate the phase factor exp(-iGR) for all atoms
    1449             :          FORALL (ci=1:neq(cn))
    1450             :          expIGR(cg, cn, ci) = EXP(-r2Pi*img*DOT_PRODUCT(taual(:, natdPtr(cn) + ci), gPts(:, cg)))
    1451             :          muffintin(cg, :, :, cn, ci) = r4Pi_sVol*expIGR(cg, cn, ci)
    1452             :          END FORALL
    1453             : 
    1454             :          ! Multiplication with i^-L
    1455             :          FORALL (cl=0:lcutm(cn))
    1456             :          muffintin(cg, :, cl*cl + 1:(cl + 1)*(cl + 1), cn, :) = muffintin(cg, :, cl*cl + 1:(cl + 1)*(cl + 1), cn, :)/imgl(cl)
    1457             :          END FORALL
    1458             : 
    1459             :          END FORALL ntypesA
    1460             : 
    1461             :          ! Calculate the k+G, abs(k+G), and abs(k+G)^2
    1462             :          k_G(:, cg) = MATMUL(gPts(:, cg) + bk(:), bmat(:, :))
    1463             :          AbsK_G2(cg) = SUM(k_G(:, cg)**2)
    1464             :          AbsK_G(cg) = SQRT(AbsK_G2(cg))
    1465             : 
    1466             :          ! Spherical harmonics are calculated for all lm's and k+G's
    1467             :          Ylm(cg, :) = calcYlm(k_G(:, cg), maxlcutm)
    1468             : 
    1469             :          ! Perform the integration in eq.[2] for all muffin-tins
    1470             :          ntypesB:FORALL (cn=1:ntype)
    1471             : 
    1472             :          ! Multiplication with the spherical harmonic
    1473             :          FORALL (cl=1:(lcutm(cn) + 1)**2)
    1474             :          muffintin(cg, :, cl, cn, :) = muffintin(cg, :, cl, cn, :)*Ylm(cg, cl)
    1475             :          END FORALL
    1476             : 
    1477             :          ! Calculate the spherical bessel function
    1478             :          FORALL (cr=1:jri(cn))
    1479             :          sphbesK_Gr(cg, cr, :, cn) = calcSphBes(AbsK_G(cg)*rmsh(cr, cn), maxlcutm)
    1480             :          END FORALL
    1481             :          ! integrate the function and multiplicate to the result
    1482             :          FORALL (cl=0:lcutm(cn))
    1483             :          FORALL (ci=1:nindxm(cl, cn))
    1484             :          intgrMT(cg, ci, cl, cn) = pure_intgrf(rmsh(:, cn)*basm(:, ci, cl, cn)* &
    1485             :                                                sphbesK_Gr(cg, :, cl, cn), jri, jmtd, rmsh, dx, ntype, cn, gridf)
    1486             :          muffintin(cg, ci, cl*cl + 1:(cl + 1)*(cl + 1), cn, :) = &
    1487             :             muffintin(cg, ci, cl*cl + 1:(cl + 1)*(cl + 1), cn, :)*intgrMT(cg, ci, cl, cn)%value
    1488             :          END FORALL
    1489             :          END FORALL
    1490             : 
    1491             :          END FORALL ntypesB
    1492             : 
    1493             :          END FORALL gpoints
    1494             : 
    1495             :          ! Check if any of the integrations failed and abort if one did
    1496           0 :          IF (ANY(intgrMT%ierror == NEGATIVE_EXPONENT_ERROR)) THEN
    1497           0 :             IF (irank == 0) WRITE (6, *) 'intgrf: Warning! Negative exponent x in extrapolation a+c*r**x'
    1498           0 :          ELSEIF (ANY(intgrMT%ierror == NEGATIVE_EXPONENT_WARNING)) THEN
    1499           0 :             IF (irank == 0) WRITE (6, *) 'intgrf: Negative exponent x in extrapolation a+c*r**x'
    1500           0 :             STOP 'intgrf: Negative exponent x in extrapolation a+c*r**x'
    1501             :          END IF
    1502             : 
    1503             : ! Calculate the Fourier transformed potential
    1504             : !                                                 /          \
    1505             : !          /      | erf |       \      4 Pi      |   |k+G|^2  |
    1506             : ! V (k) = ( k + G | --- | k + G' ) = -------  Exp| - -------  | kronecker(G,G')             [1]
    1507             : !  G       \      |  r  |       /    |k+G|^2     |    4 w^2   |
    1508             : !                                                 \          /
    1509           0 :          WHERE (absK_G /= 0)
    1510             :          ! Calculate (k+G)^2 / 4*w^2
    1511             :          arg = -r1_4omega2*AbsK_G2
    1512             : 
    1513             :          ! The potential is calculated using eq.[1]
    1514             :          potential = r4Pi/AbsK_G2*EXP(arg)
    1515             : 
    1516             : ! Calculate the Fourier transformed potential for the full(!) potential
    1517             : !
    1518             : !          /          | erfc |          \      Pi
    1519             : ! V (k) = ( k + G = 0 | ---- | k + G = 0 ) = -----
    1520             : !  G       \          |   r  |          /     w^2
    1521             : !
    1522             :          ELSEWHERE
    1523             :          ! the negative sign is added to compensate the sign when the
    1524             :          ! contribution is added to the Coulomb matrix
    1525             :          potential = -pi_omega2
    1526             :          endwhere
    1527             : 
    1528           0 :          DEALLOCATE (gridf)
    1529             : 
    1530             :          !
    1531             :          ! Create pointer which correlate the position in the array with the
    1532             :          ! appropriate indices of the MT mixed basis function
    1533             :          !
    1534           0 :          cg = 0
    1535           0 :          DO cn = 1, ntype
    1536           0 :             DO ci = 1, neq(cn)
    1537           0 :                DO cl = 0, lcutm(cn)
    1538           0 :                   DO cm = -cl, cl
    1539           0 :                      DO cr = 1, nindxm(cl, cn)
    1540           0 :                         cg = cg + 1
    1541           0 :                         ptrType(cg) = cn
    1542           0 :                         ptrEq(cg) = ci
    1543           0 :                         ptrLM(cg) = (cl + 1)*cl + cm + 1
    1544           0 :                         ptrN(cg) = cr
    1545             :                      END DO
    1546             :                   END DO
    1547             :                END DO
    1548             :             END DO
    1549             :          END DO
    1550           0 :          IF (nbasp /= cg) STOP 'hsefunctional: wrong array size: nbasp'
    1551             : 
    1552             : #ifdef CPP_INVERSION
    1553             :          ! Symmetrize muffin tin fourier transform
    1554             :          DO ci = 1, nbasp
    1555             :             sym_muffintin(ci, :noGPts) = muffintin(:, ptrN(ci), ptrLM(ci), ptrType(ci), ptrEq(ci))
    1556             :          END DO
    1557             :          DO cg = 1, noGPts
    1558             :             CALL symmetrize(sym_muffintin(:, cg), 1, nbasp, 2, .FALSE., &
    1559             :                             ntype, ntype, neq, lcutm, maxlcutm, &
    1560             :                             nindxm, natd, invsat, invsatnr)
    1561             :          END DO
    1562             :          ! store the fourier transform of the muffin tin basis
    1563             :          known_fourier_trafo(:, :, ikpt) = REAL(sym_muffintin)
    1564             : #else
    1565             :          ! store the fourier transform of the muffin tin basis
    1566           0 :          DO ci = 1, nbasp
    1567           0 :             known_fourier_trafo(ci, :noGPts, ikpt) = CONJG(muffintin(:, ptrN(ci), ptrLM(ci), ptrType(ci), ptrEq(ci)))
    1568             :          END DO
    1569             : #endif
    1570             :          ! store the fourier transform of the potential
    1571           0 :          known_potential(:noGPts, ikpt) = potential
    1572             :          ! set the flag so that the fourier transform is not calculated again
    1573           0 :          already_known(ikpt) = .TRUE.
    1574           0 :          IF (MINVAL(ABS(potential)) > 1e-12) THEN
    1575           0 :             WRITE (*, *) 'hsefunctional: Warning! Smallest potential bigger than numerical precision!', MINVAL(ABS(potential))
    1576           0 :             WRITE (*, *) 'Perhaps you should increase the number of g-Points used and recompile!'
    1577             :          END IF
    1578             : 
    1579             :       END IF
    1580             : 
    1581             :       ! return the fourier transform
    1582           0 :       potential = known_potential(:noGPts, ikpt)
    1583           0 :       fourier_trafo = known_fourier_trafo(:, :noGPts, ikpt)
    1584             : 
    1585           0 :    END SUBROUTINE calculate_fourier_transform_once
    1586             : 
    1587             :    ! Correct the pure Coulomb Matrix with by subtracting the long-range component
    1588             :    !
    1589             :    !  /     |      |      \     /     |       |      \     /     |     |      \
    1590             :    ! ( M    | V    | M     ) = ( M    | V     | M     ) - ( M    | V   | M     )
    1591             :    !  \ k,I |  HSE |  k,J /     \ k,I |  coul |  k,J /     \ k,I |  LR |  k,J /
    1592             :    !
    1593             :    ! The long-range component is given py the potential
    1594             :    !
    1595             :    !         erf( w r )
    1596             :    ! V (r) = ----------
    1597             :    !  LR         r
    1598             :    !
    1599             :    ! Input:
    1600             :    ! rmsh     - array of radial mesh points
    1601             :    ! rmt      - radius of the muffin tin
    1602             :    ! dx       - logarithmic increment of radial mesh
    1603             :    ! jri      - number of radial mesh points
    1604             :    ! jmtd     - dimension of radial mesh points
    1605             :    ! nkptf    - number of total k-points
    1606             :    ! nkptd    - dimension of k-points
    1607             :    ! nkpti    - number of irreducible k-points in window (in KPTS)
    1608             :    ! bk       - k-vector for all irreduble k-points
    1609             :    ! bmat     - reciprocal lattice vector for all directions
    1610             :    ! vol      - volume of unit cell
    1611             :    ! ntype    - number of atom types
    1612             :    ! neq      - number of symmetry-equivalent atoms of atom type i
    1613             :    ! natd     - dimension of atoms in unit cell
    1614             :    ! taual    - vector of atom positions (internal coordinates)
    1615             :    ! lcutm    - l cutoff for mixed basis
    1616             :    ! maxlcutm - maximum of all these l cutoffs
    1617             :    ! nindxm   - number of radial functions of mixed basis
    1618             :    ! maxindxm - maximum of these numbers
    1619             :    ! gptm     - reciprocal lattice vectors of the mixed basis (internal coord.)
    1620             :    ! ngptm    - number of vectors
    1621             :    ! pgptm    - pointer to the appropriate g-vector
    1622             :    ! gptmd    - dimension of gptm
    1623             :    ! basm     - radial mixed basis functions (mt + inter)
    1624             :    ! lexp     - cutoff of spherical harmonics expansion of plane wave
    1625             :    ! maxbasm  - maximum number of mixed basis functions
    1626             :    ! nbasm    - number of mixed basis function
    1627             :    ! invsat   - number of inversion-symmetric atom
    1628             :    ! invsatnr - number of inversion-symmetric atom
    1629             :    ! Inout:
    1630             :    ! coulomb  - Coulomb matrix which is changed
    1631           0 :    SUBROUTINE change_coulombmatrix( &
    1632             :       ! Input
    1633           0 :       rmsh, rmt, dx, jri, jmtd, nkptf, nkptd, nkpti, bk, &
    1634           0 :       bmat, vol, ntype, neq, natd, taual, lcutm, maxlcutm, &
    1635           0 :       nindxm, maxindxm, gptm, ngptm, pgptm, gptmd, &
    1636           0 :       basm, lexp, maxbasm, nbasm, invsat, invsatnr, irank, &
    1637             :       ! Input & output
    1638           0 :       coulomb)
    1639             : 
    1640             :       USE m_trafo, ONLY: symmetrize
    1641             :       USE m_wrapper, ONLY: packmat, unpackmat, diagonalize, inverse
    1642             :       USE m_olap, ONLY: olap_pw
    1643             : 
    1644             :       IMPLICIT NONE
    1645             : 
    1646             :       ! scalar input
    1647             :       INTEGER, INTENT(IN)    :: natd, ntype, maxlcutm, lexp
    1648             :       INTEGER, INTENT(IN)    :: jmtd, nkptf, nkptd, nkpti
    1649             :       INTEGER, INTENT(IN)    :: maxindxm
    1650             :       INTEGER, INTENT(IN)    :: gptmd, irank
    1651             :       INTEGER, INTENT(IN)    :: maxbasm
    1652             :       REAL, INTENT(IN)       :: vol
    1653             : 
    1654             :       ! array input
    1655             :       INTEGER, INTENT(IN)    :: lcutm(ntype)
    1656             :       INTEGER, INTENT(IN)    :: nindxm(0:maxlcutm, ntype), neq(ntype)
    1657             :       INTEGER, INTENT(IN)    :: jri(ntype)
    1658             :       INTEGER, INTENT(IN)    :: gptm(3, gptmd)
    1659             :       INTEGER, INTENT(IN)    :: ngptm(nkptf)
    1660             :       INTEGER, INTENT(IN)    :: pgptm(MAXVAL(ngptm), nkptf)
    1661             :       INTEGER, INTENT(IN)    :: nbasm(nkptf)
    1662             :       INTEGER, INTENT(IN)    :: invsat(natd), invsatnr(natd)
    1663             : 
    1664             :       REAL, INTENT(IN)       :: bk(3, nkptd)
    1665             :       REAL, INTENT(IN)       :: rmsh(jmtd, ntype), rmt(ntype), dx(ntype)
    1666             :       REAL, INTENT(IN)       :: basm(jmtd, maxindxm, 0:maxlcutm, ntype)
    1667             :       REAL, INTENT(IN)       :: bmat(3, 3)
    1668             :       REAL, INTENT(IN)       :: taual(3, natd)
    1669             : 
    1670             :       ! array inout
    1671             : #ifdef CPP_INVERSION
    1672             :       REAL, INTENT(INOUT)    :: coulomb(maxbasm*(maxbasm + 1)/2, nkpti)
    1673             : #else
    1674             :       COMPLEX, INTENT(INOUT) :: coulomb(maxbasm*(maxbasm + 1)/2, nkpti)
    1675             : #endif
    1676             : 
    1677             :       ! private scalars
    1678             :       INTEGER                :: ikpt, itype, ieq, il, im, iindxm, idum, n1, n2, ok    ! counters and other helper variables
    1679             :       INTEGER                :: noGPts, nbasp                                 ! no used g-vectors, no MT functions
    1680             : 
    1681             :       ! private arrays
    1682           0 :       INTEGER, ALLOCATABLE   :: ptrType(:), ptrEq(:), ptrL(:), ptrM(:), ptrN(:)  ! Pointer
    1683             :       REAL                   :: potential(maxNoGPts)                         ! Fourier transformed potential
    1684             :       COMPLEX                :: muffintin(maxNoGPts, maxindxm, &               ! muffin-tin overlap integral
    1685           0 :                                           (maxlcutm + 1)**2, ntype, MAXVAL(neq))
    1686           0 :       COMPLEX                :: interstitial(maxNoGPts, gptmd)                ! interstistial overlap intergral
    1687           0 :       COMPLEX, ALLOCATABLE   :: coulmat(:, :)                                 ! helper array to symmetrize coulomb
    1688             : 
    1689             :       ! Check size of arrays
    1690           0 :       IF (nkpti > nkptd) STOP 'hsefunctional: missmatch in dimension of arrays'
    1691           0 :       nbasp = maxbasm - MAXVAL(ngptm)
    1692           0 :       IF (ANY(nbasm - ngptm /= nbasp)) STOP 'hsefunctional: wrong assignment of nbasp'
    1693             : 
    1694             :       !
    1695             :       ! Create pointer which correlate the position in the array with the
    1696             :       ! appropriate indices of the MT mixed basis function
    1697             :       !
    1698           0 :       ALLOCATE (ptrType(nbasp), ptrEq(nbasp), ptrL(nbasp), ptrM(nbasp), ptrN(nbasp))
    1699             :       nbasp = 0
    1700           0 :       DO itype = 1, ntype
    1701           0 :          DO ieq = 1, neq(itype)
    1702           0 :             DO il = 0, lcutm(itype)
    1703           0 :                DO im = -il, il
    1704           0 :                   DO iindxm = 1, nindxm(il, itype)
    1705           0 :                      nbasp = nbasp + 1
    1706           0 :                      ptrType(nbasp) = itype
    1707           0 :                      ptrEq(nbasp) = ieq
    1708           0 :                      ptrL(nbasp) = il
    1709           0 :                      ptrM(nbasp) = im
    1710           0 :                      ptrN(nbasp) = iindxm
    1711             :                   END DO
    1712             :                END DO
    1713             :             END DO
    1714             :          END DO
    1715             :       END DO
    1716             : 
    1717             :       !
    1718             :       ! Change the Coulomb matrix for all k-points
    1719             :       !
    1720           0 :       DO ikpt = 1, nkpti
    1721             :          ! use the same g-vectors as in the mixed basis
    1722             :          ! adjust the limit of the array if necessary
    1723           0 :          IF (ngptm(ikpt) < maxNoGPts) THEN
    1724           0 :             noGPts = ngptm(ikpt)
    1725             :          ELSE
    1726           0 :             noGPts = maxNoGPts
    1727             :          END IF
    1728             : 
    1729             :          !
    1730             :          ! Calculate the Fourier transform of the mixed basis and the potential
    1731             :          !
    1732             :          CALL calculate_fourier_transform( &
    1733             :             ! Input
    1734             :             rmsh, rmt, dx, jri, jmtd, bk(:, ikpt), &
    1735             :             bmat, vol, ntype, neq, natd, taual, lcutm, maxlcutm, &
    1736             :             nindxm, maxindxm, gptm, ngptm(ikpt), pgptm(:, ikpt), gptmd, &
    1737             :             basm, noGPts, irank, &
    1738             :             ! Output
    1739           0 :             potential, muffintin, interstitial)
    1740           0 :          interstitial = CONJG(interstitial)
    1741             : 
    1742             :          ! Helper matrix for temporary storage of the attenuated Coulomb matrix
    1743           0 :          ALLOCATE (coulmat(nbasm(ikpt), nbasm(ikpt)), stat=ok)
    1744           0 :          IF (ok /= 0) STOP 'hsefunctional: failure at matrix allocation'
    1745           0 :          coulmat = 0
    1746             :          !
    1747             :          ! Calculate the difference of the Coulomb matrix by the attenuation
    1748             :          !
    1749           0 :          DO n1 = 1, nbasp
    1750           0 :             DO n2 = 1, n1
    1751             :                ! muffin tin - muffin tin contribution
    1752             :                coulmat(n2, n1) = -gPtsSummation(noGpts, &
    1753             :                                                 muffintin(:, ptrN(n2), (ptrL(n2) + 1)*ptrL(n2) + ptrM(n2) + 1, ptrType(n2), ptrEq(n2)), potential, &
    1754           0 :                                                 CONJG(muffintin(:, ptrN(n1), (ptrL(n1) + 1)*ptrL(n1) + ptrM(n1) + 1, ptrType(n1), ptrEq(n1))))
    1755           0 :                coulmat(n1, n2) = CONJG(coulmat(n2, n1))
    1756             :             END DO
    1757             :          END DO
    1758           0 :          DO n1 = nbasp + 1, nbasm(ikpt)
    1759           0 :             DO n2 = 1, n1
    1760           0 :                IF (n2 <= nbasp) THEN
    1761             :                   ! muffin tin - interstitial contribution
    1762             :                   coulmat(n2, n1) = -gPtsSummation(noGPts, &
    1763             :                                                    muffintin(:, ptrN(n2), (ptrL(n2) + 1)*ptrL(n2) + ptrM(n2) + 1, ptrType(n2), ptrEq(n2)), &
    1764           0 :                                                    potential, CONJG(interstitial(:, pgptm(n1 - nbasp, ikpt))))
    1765           0 :                   coulmat(n1, n2) = CONJG(coulmat(n2, n1))
    1766             :                ELSE
    1767             :                   ! interstitial - interstitial contribution
    1768             :                   coulmat(n2, n1) = -gPtsSummation(noGPts, &
    1769             :                                                    interstitial(:, pgptm(n2 - nbasp, ikpt)), potential, &
    1770           0 :                                                    CONJG(interstitial(:, pgptm(n1 - nbasp, ikpt))))
    1771           0 :                   coulmat(n1, n2) = CONJG(coulmat(n2, n1))
    1772             :                END IF
    1773             :             END DO
    1774             :          END DO
    1775             : 
    1776             : #ifdef CPP_INVERSION
    1777             :          ! symmetrize matrix if system has inversion symmetry
    1778             :          CALL symmetrize(coulmat, nbasm(ikpt), nbasm(ikpt), 3, .FALSE., &
    1779             :                          ntype, ntype, neq, lcutm, maxlcutm, &
    1780             :                          nindxm, natd, invsat, invsatnr)
    1781             : #endif
    1782             :          ! add the changes to the Coulomb matrix
    1783           0 :          coulomb(:nbasm(ikpt)*(nbasm(ikpt) + 1)/2, ikpt) = packmat(coulmat) + coulomb(:nbasm(ikpt)*(nbasm(ikpt) + 1)/2, ikpt)
    1784           0 :          DEALLOCATE (coulmat)
    1785             : 
    1786             :       END DO
    1787             : 
    1788           0 :       DEALLOCATE (ptrType, ptrEq, ptrL, ptrM, ptrN)
    1789             : 
    1790           0 :    END SUBROUTINE change_coulombmatrix
    1791             : 
    1792             :    ! Correct the pure Coulomb Matrix with by subtracting the long-range component
    1793             :    ! during execution time
    1794             :    !
    1795             :    !  /     |      |      \     /     |       |      \     /     |     |      \
    1796             :    ! ( M    | V    | M     ) = ( M    | V     | M     ) - ( M    | V   | M     )
    1797             :    !  \ k,I |  HSE |  k,J /     \ k,I |  coul |  k,J /     \ k,I |  LR |  k,J /
    1798             :    !
    1799             :    ! The long-range component is given py the potential
    1800             :    !
    1801             :    !         erf( w r )
    1802             :    ! V (r) = ----------
    1803             :    !  LR         r
    1804             :    !
    1805             :    ! Input:
    1806             :    ! rmsh     - array of radial mesh points
    1807             :    ! rmt      - radius of the muffin tin
    1808             :    ! dx       - logarithmic increment of radial mesh
    1809             :    ! jri      - number of radial mesh points
    1810             :    ! jmtd     - dimension of radial mesh points
    1811             :    ! bk       - k-vector for this k-point
    1812             :    ! ikpt     - number of this k-point
    1813             :    ! nkptf    - number of total k-points
    1814             :    ! bmat     - reciprocal lattice vector for all directions
    1815             :    ! vol      - volume of unit cell
    1816             :    ! ntype    - number of atom types
    1817             :    ! neq      - number of symmetry-equivalent atoms of atom type i
    1818             :    ! natd     - dimension of atoms in unit cell
    1819             :    ! taual    - vector of atom positions (internal coordinates)
    1820             :    ! lcutm    - l cutoff for mixed basis
    1821             :    ! maxlcutm - maximum of all these l cutoffs
    1822             :    ! nindxm   - number of radial functions of mixed basis
    1823             :    ! maxindxm - maximum of these numbers
    1824             :    ! gptm     - reciprocal lattice vectors of the mixed basis (internal coord.)
    1825             :    ! ngptm    - number of vectors
    1826             :    ! pgptm    - pointer to the appropriate g-vector
    1827             :    ! gptmd    - dimension of gptm
    1828             :    ! basm     - radial mixed basis functions (mt + inter)
    1829             :    ! nbasm    - number of mixed basis function
    1830             :    ! nobd     - dimension of occupied bands
    1831             :    ! nbands   - number of bands
    1832             :    ! nsst     - size of indx
    1833             :    ! indx     - pointer to bands
    1834             :    ! invsat   - number of inversion-symmetric atom
    1835             :    ! invsatnr - number of inversion-symmetric atom
    1836             :    ! cprod    - scalar product of mixed basis and wavefunction product basis
    1837             :    ! wl_iks   -
    1838             :    ! n_q      -
    1839             :    ! Return:
    1840             :    ! Change of the Coulomb matrix
    1841           0 :    FUNCTION dynamic_hse_adjustment( &
    1842           0 :       rmsh, rmt, dx, jri, jmtd, bk, ikpt, nkptf, bmat, vol, &
    1843           0 :       ntype, neq, natd, taual, lcutm, maxlcutm, nindxm, maxindxm, &
    1844           0 :       gptm, ngptm, pgptm, gptmd, basm, nbasm, &
    1845           0 :       nobd, nbands, nsst, ibando, psize, indx, invsat, invsatnr, irank, &
    1846           0 :       cprod_r, cprod_c, l_real, wl_iks, n_q)
    1847             : 
    1848             :       USE m_trafo, ONLY: symmetrize
    1849             :       USE m_olap, ONLY: olap_pw, olap_pwp
    1850             :       USE m_wrapper, ONLY: diagonalize, dotprod, matvec, packmat, inverse
    1851             : 
    1852             :       IMPLICIT NONE
    1853             : 
    1854             :       ! scalar input
    1855             :       INTEGER, INTENT(IN)  :: natd, ntype, maxlcutm
    1856             :       INTEGER, INTENT(IN)  :: jmtd, ikpt, nkptf
    1857             :       INTEGER, INTENT(IN)  :: maxindxm
    1858             :       INTEGER, INTENT(IN)  :: gptmd, irank
    1859             :       INTEGER, INTENT(IN)  :: nbasm, nobd, nbands, ibando, psize
    1860             :       INTEGER, INTENT(IN)  :: n_q
    1861             :       REAL, INTENT(IN)     :: vol
    1862             : 
    1863             :       ! array input
    1864             :       INTEGER, INTENT(IN)  :: lcutm(ntype)
    1865             :       INTEGER, INTENT(IN)  :: nindxm(0:maxlcutm, ntype), neq(ntype)
    1866             :       INTEGER, INTENT(IN)  :: jri(ntype)
    1867             :       INTEGER, INTENT(IN)  :: gptm(3, gptmd)
    1868             :       INTEGER, INTENT(IN)  :: ngptm
    1869             :       INTEGER, INTENT(IN)  :: pgptm(ngptm)
    1870             :       INTEGER, INTENT(IN)  :: nsst(nbands), indx(nbands, nbands)
    1871             :       INTEGER, INTENT(IN)  :: invsat(natd), invsatnr(natd)
    1872             :       REAL, INTENT(IN)     :: bk(3)
    1873             :       REAL, INTENT(IN)     :: rmsh(jmtd, ntype), rmt(ntype), dx(ntype)
    1874             :       REAL, INTENT(IN)     :: basm(jmtd, maxindxm, 0:maxlcutm, ntype)
    1875             :       REAL, INTENT(IN)     :: bmat(3, 3)
    1876             :       REAL, INTENT(IN)     :: taual(3, natd)
    1877             :       REAL, INTENT(IN)     :: wl_iks(nobd)
    1878             :       REAL, INTENT(IN)     :: cprod_r(nbasm, psize, nbands)
    1879             :       COMPLEX, INTENT(IN)  :: cprod_c(nbasm, psize, nbands)
    1880             :       LOGICAL, INTENT(IN)   :: l_real
    1881             : 
    1882             :       ! return type definition
    1883             :       COMPLEX              :: dynamic_hse_adjustment(nbands, nbands)
    1884             : 
    1885             :       ! private scalars
    1886             :       INTEGER              :: noGpts, nbasp
    1887             :       INTEGER              :: itype, ieq, il, im, iindxm, ibasp
    1888             :       INTEGER              :: igpt, iobd, iobd0, isst, iband1, iband2
    1889             :       COMPLEX              :: cdum, gPtsSum
    1890             : 
    1891             :       ! private arrays
    1892             :       INTEGER, ALLOCATABLE :: ptrType(:), ptrEq(:), ptrL(:), ptrM(:), ptrN(:) ! pointer to reference muffintin-array
    1893             :       REAL                 :: potential(maxNoGPts)                        ! Fourier transformed potential
    1894             :       COMPLEX              :: muffintin(maxNoGPts, maxindxm, &              ! muffin-tin overlap integral
    1895             :                                         (maxlcutm + 1)**2, ntype, MAXVAL(neq))
    1896             :       COMPLEX              :: interstitial(maxNoGPts, gptmd)               ! interstistial overlap intergral
    1897             : #ifdef CPP_INVERSION
    1898             :       REAL                 :: fourier_trafo(nbasm - ngptm, maxNoGPts)        ! Fourier trafo of all mixed basis functions
    1899             : #else
    1900           0 :       COMPLEX              :: fourier_trafo(nbasm - ngptm, maxNoGPts)        ! Fourier trafo of all mixed basis functions
    1901             : #endif
    1902             : 
    1903           0 :       REAL                 :: cprod_fourier_trafo_r(maxNoGpts, psize, nbands)  ! Product of cprod and Fourier trafo
    1904           0 :       COMPLEX              :: cprod_fourier_trafo_c(maxNoGpts, psize, nbands)  ! Product of cprod and Fourier trafo
    1905             : 
    1906             :       ! Initialisation
    1907           0 :       dynamic_hse_adjustment = 0.0
    1908           0 :       potential = 0.0
    1909           0 :       fourier_trafo = 0.0
    1910           0 :       cprod_fourier_trafo_r = 0.0
    1911           0 :       cprod_fourier_trafo_c = 0.0
    1912           0 :       noGPts = MIN(ngptm, maxNoGPts)
    1913           0 :       nbasp = nbasm - ngptm
    1914             : 
    1915             :       !
    1916             :       ! Calculate the fourier transform of the mixed basis once
    1917             :       ! If it was already calculated load the old results
    1918             :       !
    1919             :       CALL calculate_fourier_transform_once( &
    1920             :          rmsh, rmt, dx, jri, jmtd, bk, ikpt, nkptf, &
    1921             :          bmat, vol, ntype, neq, natd, taual, lcutm, maxlcutm, &
    1922             :          nindxm, maxindxm, gptm, ngptm, pgptm, gptmd, &
    1923             :          nbasp, basm, noGPts, invsat, invsatnr, irank, &
    1924           0 :          potential, fourier_trafo)
    1925             : 
    1926             :       ! Calculate the Fourier transform of the 'normal' basis
    1927             :       ! by summing over the mixed basis
    1928           0 :       DO igpt = 1, noGPts
    1929           0 :          DO iobd0 = 1, psize!ibando,min(ibando+psize,nobd)
    1930           0 :             iobd = iobd0 + ibando - 1
    1931           0 :             IF (iobd > nobd) CYCLE
    1932           0 :             DO iband1 = 1, nbands
    1933           0 :                if (l_real) THEN
    1934             :                   cprod_fourier_trafo_r(igpt, iobd0, iband1) = &
    1935             :                      ! muffin tin contribution
    1936             :                      dotprod(fourier_trafo(:nbasp, igpt), cprod_r(:nbasp, iobd0, iband1)) &
    1937             :                      ! interstitial contribution (interstitial is kronecker_G,G')
    1938           0 :                      + cprod_r(nbasp + igpt, iobd0, iband1)
    1939             :                else
    1940             :                   cprod_fourier_trafo_c(igpt, iobd0, iband1) = &
    1941             :                      ! muffin tin contribution
    1942             :                      dotprod(cprod_c(:nbasp, iobd0, iband1), fourier_trafo(:nbasp, igpt)) &
    1943             :                      ! interstitial contribution (interstitial is kronecker_G,G')
    1944           0 :                      + CONJG(cprod_c(nbasp + igpt, iobd0, iband1))
    1945             :                endif
    1946             :             END DO
    1947             :          END DO
    1948             :       END DO
    1949             : 
    1950             :       ! Summation over all G-vectors
    1951           0 :       DO iband1 = 1, nbands
    1952           0 :          DO iobd0 = 1, psize!ibando,min(ibando+psize,nobd)
    1953           0 :             iobd = iobd0 + ibando - 1
    1954           0 :             IF (iobd > nobd) CYCLE
    1955             :             ! prefactor for all elements
    1956           0 :             cdum = wl_iks(iobd)/n_q
    1957           0 :             DO isst = 1, nsst(iband1)
    1958           0 :                iband2 = indx(isst, iband1)
    1959             :                ! do summation over G-vectors
    1960           0 :                if (l_real) THEN
    1961             :                   gPtsSum = cdum*gPtsSummation(noGPts, &
    1962             :                                                cprod_fourier_trafo_r(:, iobd0, iband1), &
    1963           0 :                                                potential, cprod_fourier_trafo_r(:, iobd0, iband2))
    1964             :                ELSE
    1965             :                   gPtsSum = cdum*gPtsSummation(noGPts, &
    1966             :                                                conjg(cprod_fourier_trafo_c(:, iobd0, iband1)), &
    1967           0 :                                                potential, cprod_fourier_trafo_c(:, iobd0, iband2))
    1968             :                endif
    1969             :                dynamic_hse_adjustment(iband2, iband1) = &
    1970           0 :                   dynamic_hse_adjustment(iband2, iband1) - gPtsSum
    1971             :             END DO
    1972             :          END DO
    1973             :       END DO
    1974             : 
    1975             :    END FUNCTION dynamic_hse_adjustment
    1976             : 
    1977             :    ! Helper function needed to use forall statements instead of do loop's
    1978             :    ! calls the harmonicsr subroutine from 'util.F'
    1979           0 :    PURE FUNCTION calcYlm(rvec, ll)
    1980             :       USE m_util, ONLY: harmonicsr
    1981             :       IMPLICIT NONE
    1982             :       REAL, INTENT(IN)    :: rvec(3)
    1983             :       INTEGER, INTENT(IN) :: ll
    1984             :       COMPLEX             :: calcYlm((ll + 1)**2)
    1985           0 :       CALL harmonicsr(calcYlm, rvec, ll)
    1986             :    END FUNCTION calcYlm
    1987             : 
    1988             :    ! Helper function needed to use forall statements instead of do loop's
    1989             :    ! calls the sphbessel subroutine from 'util.F'
    1990           0 :    PURE FUNCTION calcSphBes(x, l)
    1991             :       USE m_util, ONLY: sphbessel
    1992             :       IMPLICIT NONE
    1993             :       REAL, INTENT(IN)    :: x
    1994             :       INTEGER, INTENT(IN) :: l
    1995             :       REAL                :: calcSphBes(0:l)
    1996           0 :       CALL sphbessel(calcSphBes, x, l)
    1997             :    END FUNCTION calcSphBes
    1998             : 
    1999             :    ! Dot_product definition for two arrays (3d x 1d) with common size of first dimension
    2000           0 :    PURE FUNCTION my_dot_product3x1(first, second)
    2001             :       IMPLICIT NONE
    2002             :       INTEGER, INTENT(IN) :: first(:, :, :)
    2003             :       REAL, INTENT(IN)    :: second(SIZE(first, 1))
    2004             :       REAL                :: my_dot_product3x1(SIZE(first, 2), SIZE(first, 3))
    2005             :       INTEGER             :: ci, cj
    2006             : 
    2007           0 :       FORALL (ci=1:SIZE(first, 2), cj=1:SIZE(first, 3))
    2008             :       my_dot_product3x1(ci, cj) = DOT_PRODUCT(first(:, ci, cj), second)
    2009             :       END FORALL
    2010             :    END FUNCTION my_dot_product3x1
    2011             : 
    2012             :    ! Dot_product definition for two arrays (4d x 1d) with common size of first dimension
    2013           0 :    PURE FUNCTION my_dot_product4x1(first, second)
    2014             :       IMPLICIT NONE
    2015             :       INTEGER, INTENT(IN) :: first(:, :, :, :)
    2016             :       REAL, INTENT(IN)    :: second(SIZE(first, 1))
    2017             :       REAL                :: my_dot_product4x1(SIZE(first, 2), SIZE(first, 3), SIZE(first, 4))
    2018             :       INTEGER             :: ci, cj, ck
    2019             : 
    2020           0 :       FORALL (ci=1:SIZE(first, 2), cj=1:SIZE(first, 3), ck=1:SIZE(first, 4))
    2021             :       my_dot_product4x1(ci, cj, ck) = DOT_PRODUCT(first(:, ci, cj, ck), second)
    2022             :       END FORALL
    2023             :    END FUNCTION my_dot_product4x1
    2024             : 
    2025             :    ! Dot_product definition for two arrays (3d x 3d) with common size of all dimensions,
    2026             :    ! The scalar product is evaluated over the last index
    2027           0 :    PURE FUNCTION my_dot_product3x3(first, second)
    2028             :       IMPLICIT NONE
    2029             :       COMPLEX, INTENT(IN) :: first(:, :, :)
    2030             :       REAL, INTENT(IN) :: second(SIZE(first, 1), SIZE(first, 2), SIZE(first, 3))
    2031             :       COMPLEX             :: my_dot_product3x3(SIZE(first, 1), SIZE(first, 2))
    2032             :       INTEGER             :: ci, cj
    2033             : 
    2034           0 :       FORALL (ci=1:SIZE(first, 1), cj=1:SIZE(first, 2))
    2035             :       my_dot_product3x3(ci, cj) = DOT_PRODUCT(first(ci, cj, :), second(ci, cj, :))
    2036             :       END FORALL
    2037             :    END FUNCTION my_dot_product3x3
    2038             : 
    2039             :    ! Dot_product definition for two arrays (4d x 4d) with common size of all dimensions,
    2040             :    ! The scalar product is evaluated over the last index
    2041           0 :    PURE FUNCTION my_dot_product4x4(first, second)
    2042             :       IMPLICIT NONE
    2043             :       COMPLEX, INTENT(IN) :: first(:, :, :, :)
    2044             :       REAL, INTENT(IN) :: second(SIZE(first, 1), SIZE(first, 2), SIZE(first, 3), SIZE(first, 4))
    2045             :       COMPLEX             :: my_dot_product4x4(SIZE(first, 1), SIZE(first, 2), SIZE(first, 3))
    2046             :       INTEGER             :: ci, cj, ck
    2047             : 
    2048           0 :       FORALL (ci=1:SIZE(first, 1), cj=1:SIZE(first, 2), ck=1:SIZE(first, 3))
    2049             :       my_dot_product4x4(ci, cj, ck) = DOT_PRODUCT(first(ci, cj, ck, :), second(ci, cj, ck, :))
    2050             :       END FORALL
    2051             :    END FUNCTION my_dot_product4x4
    2052             : 
    2053             :    ! Sum over last index of 4d array
    2054           0 :    PURE FUNCTION my_sum3d(array)
    2055             :       IMPLICIT NONE
    2056             :       COMPLEX, INTENT(IN) :: array(:, :, :)
    2057             :       COMPLEX             :: my_sum3d(SIZE(array, 1), SIZE(array, 2))
    2058             :       INTEGER             :: ci, cj
    2059             : 
    2060           0 :       FORALL (ci=1:SIZE(array, 1), cj=1:SIZE(array, 2))
    2061             :       my_sum3d(ci, cj) = SUM(array(ci, cj, :))
    2062             :       END FORALL
    2063             :    END FUNCTION my_sum3d
    2064             : 
    2065             :    ! Sum over last index of 4d array
    2066           0 :    PURE FUNCTION my_sum4d(array)
    2067             :       IMPLICIT NONE
    2068             :       COMPLEX, INTENT(IN) :: array(:, :, :, :)
    2069             :       COMPLEX             :: my_sum4d(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3))
    2070             :       INTEGER             :: ci, cj, ck
    2071             : 
    2072           0 :       FORALL (ci=1:SIZE(array, 1), cj=1:SIZE(array, 2), ck=1:SIZE(array, 3))
    2073             :       my_sum4d(ci, cj, ck) = SUM(array(ci, cj, ck, :))
    2074             :       END FORALL
    2075             :    END FUNCTION my_sum4d
    2076             : 
    2077             :    ! Sum product of three vectors
    2078           0 :    COMPLEX PURE FUNCTION crc_gPtsSummation(noGpts, vec1, vec2, vec3)
    2079             :       IMPLICIT NONE
    2080             :       INTEGER, INTENT(IN) :: noGPts
    2081             :       COMPLEX, INTENT(IN) :: vec1(noGPts)
    2082             :       REAL, INTENT(IN)    :: vec2(noGPts)
    2083             :       COMPLEX, INTENT(IN) :: vec3(noGPts)
    2084           0 :       COMPLEX             :: temp(noGPts)
    2085           0 :       temp = vec1*vec3
    2086           0 :       crc_gPtsSummation = DOT_PRODUCT(temp, vec2)
    2087           0 :    END FUNCTION crc_gPtsSummation
    2088           0 :    REAL PURE FUNCTION rrr_gPtsSummation(noGpts, vec1, vec2, vec3)
    2089             :       IMPLICIT NONE
    2090             :       INTEGER, INTENT(IN) :: noGPts
    2091             :       REAL, INTENT(IN)    :: vec1(noGPts)
    2092             :       REAL, INTENT(IN)    :: vec2(noGPts)
    2093             :       REAL, INTENT(IN)    :: vec3(noGPts)
    2094           0 :       REAL                :: temp(noGPts)
    2095           0 :       temp = vec1*vec3
    2096           0 :       rrr_gPtsSummation = DOT_PRODUCT(temp, vec2)
    2097           0 :    END FUNCTION rrr_gPtsSummation
    2098             : 
    2099             :    ! Include the core valence exchange for the HSE functional
    2100             :    ! The HSE potential is the Coulomb one attenuated with the complementary
    2101             :    ! error function which can be expanded in Legendre polynomials (r' > r)
    2102             :    !
    2103             :    !                      -----                l+2n
    2104             :    ! erfc( w |r - r'| )    \                  r           /         \
    2105             :    ! ------------------ =   )    d  (w r')  ---------  P ( cos(r,r') )
    2106             :    !    | r - r' |         /      ln          l+2n+1    l \         /
    2107             :    !                      -----              r'
    2108             :    !                       l,n
    2109             :    !
    2110             :    ! Then an analytical integration of the angular part of the exchange
    2111             :    ! potential is possible and only the radial part remains
    2112             :    !
    2113             :    !          R                        r                              R
    2114             :    !           /   -----  /  d (w r)    /                              /           d  (w r')  \
    2115             :    !  4 Pi    |     \    |    ln       |            l+2n+2      l+2n  |             ln         |
    2116             :    ! -------  | dr   )   |  ---------  | dr' f(r') r'       +  r      | dr' f(r') -----------  |
    2117             :    ! 2 l + 1  |     /    |    l+2n+1   |                              |              l+2n-1    |
    2118             :    !         /     -----  \  r        /                              /              r'        /
    2119             :    !          0      n                 0                              r
    2120             :    !
    2121             :    ! The function f is a product of core and valence function
    2122             :    ! (note: in the code f is defined with an additional 1/r factor)
    2123             :    !
    2124           0 :    SUBROUTINE exchange_vccvHSE(nk, bkpt, nkptd, nkpt, nkpti, ntype, neq, natd, &
    2125           0 :                                lmax, lmaxd, nindx, maxindx, lmaxc, nindxc, &
    2126           0 :                                maxindxc, core1, core2, lcutm, maxlmindx, &
    2127           0 :                                bas1, bas2, jmtd, rmsh, dx, jri, jspd, jsp, &
    2128           0 :                                maxfac, fac, sfac, nv, neigd, nbasfcn, nbands, &
    2129           0 :                                gridf, nsymop, nsest, indx_sest, irank, &
    2130           0 :                                a_ex, nobd, w_iks, &
    2131           0 :                                mat_ex, te_hfex_core)
    2132             : 
    2133             :       USE m_constants
    2134             :       USE m_util
    2135             :       USE m_wrapper
    2136             : 
    2137             :       IMPLICIT NONE
    2138             : 
    2139             : !   -scalars -
    2140             :       INTEGER, INTENT(IN)      ::  ntype, jmtd, lmaxd, jspd, jsp, neigd
    2141             :       INTEGER, INTENT(IN)      ::  maxfac, nbands, maxlmindx, natd, &
    2142             :                                   maxindxc, nk, maxindx, nbasfcn
    2143             :       INTEGER, INTENT(IN)      ::  nkptd, nkpt, nkpti, irank
    2144             :       INTEGER, INTENT(IN)      ::  nsymop
    2145             :       REAL, INTENT(IN)      ::  a_ex
    2146             :       REAL, INTENT(INOUT)   ::  te_hfex_core
    2147             : 
    2148             : !   - arrays -
    2149             :       INTEGER, INTENT(IN)      ::  neq(ntype), lcutm(ntype), lmax(ntype), &
    2150             :                                   lmaxc(ntype), jri(ntype), nv(jspd), &
    2151             :                                   nindxc(0:MAXVAL(lmaxc), ntype), &
    2152             :                                   nindx(0:lmaxd, ntype)
    2153             :       INTEGER, INTENT(IN)      ::  nsest(nbands), indx_sest(nbands, nbands)
    2154             :       INTEGER, INTENT(IN)      ::  nobd(nkpt)
    2155             :       REAL, INTENT(IN)         ::  rmsh(jmtd, ntype), dx(ntype)
    2156             :       REAL, INTENT(IN)         ::  bas1(jmtd, maxindx, 0:lmaxd, ntype), &
    2157             :                                   bas2(jmtd, maxindx, 0:lmaxd, ntype)
    2158             :       REAL, INTENT(IN)         ::  core1(jmtd, maxindxc, 0:MAXVAL(lmaxc), ntype), &
    2159             :                                   core2(jmtd, maxindxc, 0:MAXVAL(lmaxc), ntype)
    2160             :       REAL, INTENT(IN)         ::  fac(0:maxfac), sfac(0:maxfac)
    2161             :       REAL, INTENT(IN)         ::  bkpt(3)
    2162             :       REAL, INTENT(IN)         ::  gridf(jmtd, ntype)
    2163             :       REAL, INTENT(IN)         ::  w_iks(neigd, nkptd, jspd)
    2164             : 
    2165             : #ifdef CPP_INVERSION
    2166             :       REAL, INTENT(INOUT)  ::  mat_ex(nbasfcn*(nbasfcn + 1)/2)
    2167             : #else
    2168             :       COMPLEX, INTENT(INOUT)  ::  mat_ex(nbasfcn*(nbasfcn + 1)/2)
    2169             : #endif
    2170             :       INTEGER, PARAMETER      ::  ncut = 5                  ! cut-off value of n-summation
    2171             :       INTEGER                 ::  cn                        ! counter for n-summation
    2172           0 :       REAL                    ::  d_ln(jmtd, 0:lmaxd, 0:ncut) ! expansion coefficients of erfc(wr)/r
    2173             :       ! in Legendre polynomials
    2174             : 
    2175             : !   - local scalars -
    2176             :       INTEGER                 ::  iatom, ieq, itype, ic, l, l1, l2, &
    2177             :                                  ll, lm, m, m1, m2, p1, p2, n, n1, n2, nn2, i, j
    2178             :       INTEGER                 ::  iband1, iband2, ndb1, ndb2, ic1, ic2
    2179             :       INTEGER                 ::  irecl_cmt
    2180             : 
    2181             :       REAL                    ::  rdum
    2182             :       REAL                    ::  sum_offdia
    2183             : 
    2184             :       COMPLEX                 ::  cdum
    2185             : !   - local arrays -
    2186           0 :       INTEGER, ALLOCATABLE     ::  larr(:), larr2(:)
    2187           0 :       INTEGER, ALLOCATABLE     ::  parr(:), parr2(:)
    2188             : 
    2189           0 :       REAL                    ::  sum_primf(jmtd), integrand(jmtd)
    2190           0 :       REAL                    ::  primf1(jmtd), primf2(jmtd)
    2191           0 :       REAL, ALLOCATABLE        ::  fprod(:, :), fprod2(:, :)
    2192           0 :       REAL, ALLOCATABLE        ::  integral(:, :)
    2193             : 
    2194           0 :       COMPLEX                 ::  cmt(neigd, maxlmindx, natd)
    2195           0 :       COMPLEX                 ::  exchange(nbands, nbands)
    2196           0 :       COMPLEX, ALLOCATABLE     ::  carr(:, :), carr2(:, :), carr3(:, :)
    2197             : 
    2198             :       LOGICAL                 ::  ldum(nbands, nbands)
    2199             : 
    2200             :       ! check if a_ex is consistent
    2201             : !     IF ( a_ex /= aMix_HSE ) STOP 'hsefunctional: inconsistent mixing!'
    2202             : 
    2203             :       ! read in mt wavefunction coefficients from file cmt
    2204           0 :       irecl_cmt = neigd*maxlmindx*natd*16
    2205           0 :       OPEN (unit=777, file='cmt', form='unformatted', access='direct', recl=irecl_cmt)
    2206           0 :       READ (777, rec=nk) cmt(:, :, :)
    2207           0 :       CLOSE (777)
    2208             : 
    2209           0 :       ALLOCATE (fprod(jmtd, 5), larr(5), parr(5))
    2210             : 
    2211           0 :       exchange = 0
    2212           0 :       iatom = 0
    2213           0 :       rdum = 0
    2214           0 :       DO itype = 1, ntype
    2215             :          ! Calculate the expansion coefficients of the potential in Legendre polynomials
    2216           0 :          d_ln(:, :lcutm(itype), :) = calculate_coefficients(rmsh(:, itype), lcutm(itype), ncut, fac)
    2217           0 :          DO ieq = 1, neq(itype)
    2218           0 :             iatom = iatom + 1
    2219           0 :             DO l1 = 0, lmaxc(itype)
    2220           0 :                DO p1 = 1, nindxc(l1, itype)
    2221             : 
    2222           0 :                   DO l = 0, lcutm(itype)
    2223             : 
    2224             :                      ! Define core-valence product functions
    2225             : 
    2226           0 :                      n = 0
    2227           0 :                      DO l2 = 0, lmax(itype)
    2228           0 :                         IF (l < ABS(l1 - l2) .OR. l > l1 + l2) CYCLE
    2229             : 
    2230           0 :                         DO p2 = 1, nindx(l2, itype)
    2231           0 :                            n = n + 1
    2232           0 :                            m = SIZE(fprod, 2)
    2233           0 :                            IF (n > m) THEN
    2234           0 :                               ALLOCATE (fprod2(jmtd, m), larr2(m), parr2(m))
    2235           0 :                               fprod2 = fprod; larr2 = larr; parr2 = parr
    2236           0 :                               DEALLOCATE (fprod, larr, parr)
    2237           0 :                               ALLOCATE (fprod(jmtd, m + 5), larr(m + 5), parr(m + 5))
    2238           0 :                               fprod(:, :m) = fprod2
    2239           0 :                               larr(:m) = larr2
    2240           0 :                               parr(:m) = parr2
    2241           0 :                               DEALLOCATE (fprod2, larr2, parr2)
    2242             :                            END IF
    2243             :                            fprod(:, n) = (core1(:, p1, l1, itype) &
    2244             :                                           *bas1(:, p2, l2, itype) &
    2245             :                                           + core2(:, p1, l1, itype) &
    2246           0 :                                           *bas2(:, p2, l2, itype))/rmsh(:, itype)
    2247           0 :                            larr(n) = l2
    2248           0 :                            parr(n) = p2
    2249             :                         END DO
    2250             :                      END DO
    2251             : 
    2252             :                      ! Evaluate radial integrals (special part of Coulomb matrix : contribution from single MT)
    2253             : 
    2254             :                      ALLOCATE (integral(n, n), carr(n, nbands), &
    2255           0 :                                carr2(n, nv(jsp)), carr3(n, nv(jsp)))
    2256             : 
    2257           0 :                      DO i = 1, n
    2258             :                         ! Initialization of n Summation
    2259           0 :                         sum_primf = 0.0
    2260             :                         ! Integration over r'
    2261           0 :                         DO cn = 0, ncut ! Summation over Legendre polynomials
    2262           0 :                            primf1 = 0.0
    2263           0 :                            primf2 = 0.0
    2264             :                            ! Calculate integral for 0 < r' < r
    2265             :                            CALL primitivef(primf1, fprod(:, i)*rmsh(:, itype)**(l + 2*cn + 1), &
    2266           0 :                                            rmsh, dx, jri, jmtd, itype, ntype)
    2267             :                            ! Calculate integral for r < r' < R
    2268             :                            CALL primitivef(primf2, d_ln(:, l, cn)*fprod(:, i)/rmsh(:, itype)**(l + 2*cn), &
    2269           0 :                                            rmsh, dx, jri, jmtd, -itype, ntype)  ! -itype is to enforce inward integration
    2270             :                            ! Multiplication with appropriate prefactors
    2271           0 :                            primf1 = primf1/rmsh(:, itype)**(l + 2*cn)*d_ln(:, l, cn)
    2272           0 :                            primf2 = primf2*rmsh(:, itype)**(l + 2*cn + 1)
    2273             :                            ! Summation over n
    2274           0 :                            sum_primf = sum_primf + primf1 + primf2
    2275             :                         END DO
    2276           0 :                         DO j = 1, n
    2277             :                            ! Integration over r
    2278           0 :                            integrand = fprod(:, j)*sum_primf
    2279             :                            integral(i, j) = fpi_const/(2*l + 1)* &
    2280             :                                             intgrf(integrand, jri, jmtd, rmsh, &
    2281           0 :                                                    dx, ntype, itype, gridf)
    2282             :                         END DO
    2283             : 
    2284             :                      END DO
    2285             : 
    2286             :                      ! Add everything up
    2287             : 
    2288           0 :                      DO m1 = -l1, l1
    2289           0 :                         DO m = -l, l
    2290           0 :                            m2 = m1 + m
    2291             : 
    2292           0 :                            carr = 0
    2293             :                            ic = 0
    2294           0 :                            DO n1 = 1, nbands
    2295             : 
    2296           0 :                               DO i = 1, n
    2297           0 :                                  ll = larr(i)
    2298           0 :                                  IF (ABS(m2) > ll) CYCLE
    2299             : 
    2300             :                                  lm = SUM((/((2*l2 + 1)*nindx(l2, itype), l2=0, ll - 1)/)) &
    2301           0 :                                       + (m2 + ll)*nindx(ll, itype) + parr(i)
    2302             : 
    2303             :                                  carr(i, n1) = cmt(n1, lm, iatom) &
    2304           0 :                                                *gaunt(l1, ll, l, m1, m2, m, maxfac, fac, sfac)
    2305             : 
    2306             :                               END DO
    2307           0 :                               DO n2 = 1, nsest(n1)!n1
    2308           0 :                                  nn2 = indx_sest(n2, n1)
    2309             :                                  exchange(nn2, n1) = exchange(nn2, n1) &
    2310           0 :                                                      + DOT_PRODUCT(carr(:, n1), MATMUL(integral, carr(:, nn2)))
    2311             :                               END DO
    2312             :                            END DO
    2313             :                         END DO
    2314             :                      END DO
    2315             : 
    2316           0 :                      DEALLOCATE (integral, carr, carr2, carr3)
    2317             : 
    2318             :                   END DO
    2319             :                END DO
    2320             :             END DO
    2321             :          END DO
    2322             :       END DO
    2323             : 
    2324             : #ifdef CPP_INVERSION
    2325             :       IF (ANY(ABS(aimag(exchange)) > 10.0**-10)) THEN
    2326             :          IF (irank == 0) WRITE (6, '(A)') 'exchangeCore: Warning! Unusually large imaginary component.'
    2327             :          WRITE (*, *) MAXVAL(ABS(aimag(exchange)))
    2328             :          STOP 'exchangeCore: Unusually large imaginary component.'
    2329             :       END IF
    2330             : #endif
    2331             : 
    2332           0 :       DO n1 = 1, nobd(nk)
    2333           0 :          te_hfex_core = te_hfex_core - a_ex*w_iks(n1, nk, jsp)*exchange(n1, n1)
    2334             :       END DO
    2335             : 
    2336             :       ! add the core-valence contribution to the exchange matrix mat_ex
    2337             : 
    2338             :       ic = 0
    2339             :       sum_offdia = 0
    2340           0 :       DO n1 = 1, nbands
    2341           0 :          DO n2 = 1, n1
    2342           0 :             ic = ic + 1
    2343           0 :             mat_ex(ic) = mat_ex(ic) + CONJG(exchange(n2, n1))/nsymop
    2344             :          END DO
    2345             :       END DO
    2346             : 
    2347           0 :    END SUBROUTINE exchange_vccvHSE
    2348             : 
    2349             :    ! Include the core core exchange for the HSE functional
    2350             :    ! The HSE potential is the Coulomb one attenuated with the complementary
    2351             :    ! error function which can be expanded in Legendre polynomials (r' > r)
    2352             :    !
    2353             :    !                      -----                l+2n
    2354             :    ! erfc( w |r - r'| )    \                  r           /         \
    2355             :    ! ------------------ =   )    d  (w r')  ---------  P ( cos(r,r') )
    2356             :    !    | r - r' |         /      ln          l+2n+1    l \         /
    2357             :    !                      -----              r'
    2358             :    !                       l,n
    2359             :    !
    2360             :    ! Then an analytical integration of the angular part of the exchange
    2361             :    ! potential is possible and only the radial part remains
    2362             :    !
    2363             :    !          R                        r                              R
    2364             :    !           /   -----  /  d (w r)    /                              /           d  (w r')  \
    2365             :    !  4 Pi    |     \    |    ln       |            l+2n+2      l+2n  |             ln         |
    2366             :    ! -------  | dr   )   |  ---------  | dr' f(r') r'       +  r      | dr' f(r') -----------  |
    2367             :    ! 2 l + 1  |     /    |    l+2n+1   |                              |              l+2n-1    |
    2368             :    !         /     -----  \  r        /                              /              r'        /
    2369             :    !          0      n                 0                              r
    2370             :    !
    2371             :    ! The function f is a product of core and valence function
    2372             :    ! (note: in the code f is defined with an additional 1/r factor)
    2373             :    !
    2374           0 :    SUBROUTINE exchange_ccccHSE( &
    2375             :       ! Input
    2376           0 :       nk, nkpti, nw, nwd, ntype, neq, natd, lmaxcd, lmaxc, &
    2377           0 :       nindxc, maxindxc, ncst, ncstd, jmtd, jri, &
    2378           0 :       rmsh, dx, lmaxd, core1, core2, bkpt, gridf, &
    2379           0 :       invsat, invsatnr, wtkpt, maxfac, fac, a_ex, irank, &
    2380             :       ! Output
    2381             :       te_hfex)
    2382             : 
    2383             :       USE m_constants
    2384             :       USE m_util
    2385             :       USE m_wrapper
    2386             :       USE m_gaunt
    2387             :       USE m_trafo
    2388             : 
    2389             :       IMPLICIT NONE
    2390             : 
    2391             :       ! - scalars -
    2392             :       INTEGER, INTENT(IN)    ::  nk, nkpti, nw, nwd, ntype, natd, ncstd
    2393             :       INTEGER, INTENT(IN)    ::  lmaxcd, lmaxd, irank
    2394             :       INTEGER, INTENT(IN)    ::  jmtd, maxindxc, maxfac
    2395             : 
    2396             :       REAL, INTENT(IN)    ::  a_ex
    2397             :       REAL, INTENT(INOUT) ::  te_hfex
    2398             : 
    2399             :       ! - arays -
    2400             :       INTEGER, INTENT(IN)    ::  neq(ntype), ncst(ntype), lmaxc(ntype)
    2401             :       INTEGER, INTENT(IN)    ::  nindxc(0:lmaxcd, ntype)
    2402             :       INTEGER, INTENT(IN)    ::  jri(ntype)
    2403             :       INTEGER, INTENT(IN)    ::  invsat(natd), invsatnr(natd)
    2404             : 
    2405             :       REAL, INTENT(IN)    ::  rmsh(jmtd, ntype), dx(ntype)
    2406             :       REAL, INTENT(IN)    ::  core1(jmtd, maxindxc, 0:lmaxcd, ntype), &
    2407             :                              core2(jmtd, maxindxc, 0:lmaxcd, ntype)
    2408             :       REAL, INTENT(IN)    ::  bkpt(3)
    2409             :       REAL, INTENT(IN)    ::  gridf(jmtd, ntype)
    2410             :       REAL, INTENT(IN)    ::  wtkpt(nkpti, nwd)
    2411             :       REAL, INTENT(IN)    ::  fac(0:maxfac)
    2412             : 
    2413             :       ! - local scalars -
    2414             :       INTEGER               ::  itype, ieq, icst, icst1, icst2, iatom, iatom0
    2415             :       INTEGER               ::  l1, l2, l, ll, llmax
    2416             :       INTEGER               ::  m1, m2, m, mm
    2417             :       INTEGER               ::  n1, n2, n
    2418             : 
    2419             :       REAL                  ::  rdum, rdum1
    2420             :       ! - local arrays -
    2421           0 :       INTEGER               ::  point(maxindxc, -lmaxcd:lmaxcd, 0:lmaxcd, natd)
    2422           0 :       REAL                  ::  rprod(jmtd), primf1(jmtd), primf2(jmtd), sum_primf(jmtd), integrand(jmtd)
    2423           0 :       COMPLEX               ::  exch(ncstd, ncstd)
    2424             : 
    2425             :       INTEGER, PARAMETER    ::  ncut = 5                     ! cut-off value of n-summation
    2426             :       INTEGER               ::  cn                           ! counter for n-summation
    2427           0 :       REAL                  ::  d_ln(jmtd, 0:2*lmaxcd, 0:ncut) ! expansion coefficients of erfc(wr)/r
    2428             :       ! in Legendre polynomials
    2429             :       CHARACTER*100         :: outtext
    2430             : 
    2431             : !   IF ( a_ex /= aMix_HSE ) STOP 'hsefunctional: mixing parameter inconsistent'
    2432             : 
    2433             : !   IF ( irank == 0 ) THEN
    2434             : !     WRITE(outtext,'(A)') new_line('n') // new_line('n') // '### core-core-core-core exchange ###'
    2435             : !     CALL writeout(outtext,irank)
    2436             : !     WRITE(outtext,'(A)') new_line('n') // '        k-point       band    exchange'
    2437             : !     CALL writeout(outtext,irank)
    2438             : !   END IF
    2439             : 
    2440             :       ! set up point
    2441           0 :       icst = 0
    2442           0 :       iatom = 0
    2443           0 :       DO itype = 1, ntype
    2444           0 :          DO ieq = 1, neq(itype)
    2445           0 :             iatom = iatom + 1
    2446           0 :             DO l = 0, lmaxc(itype)
    2447           0 :                DO m = -l, l
    2448           0 :                   DO n = 1, nindxc(l, itype)
    2449           0 :                      icst = icst + 1
    2450           0 :                      point(n, m, l, iatom) = icst
    2451             :                   END DO
    2452             :                END DO
    2453             :             END DO
    2454             :          END DO
    2455             :       END DO
    2456             : 
    2457           0 :       llmax = 2*lmaxcd
    2458           0 :       exch = 0
    2459           0 :       iatom0 = 0
    2460           0 :       DO itype = 1, ntype
    2461             :          ! Calculate the expansion coefficients of the potential in Legendre polynomials
    2462           0 :          d_ln(:, 0:2*lmaxc(itype), :) = calculate_coefficients(rmsh(:, itype), 2*lmaxc(itype), ncut, fac)
    2463             : 
    2464           0 :          DO l1 = 0, lmaxc(itype)  ! left core state
    2465             : 
    2466           0 :             DO l2 = 0, lmaxc(itype)  ! right core state
    2467           0 :                DO l = 0, lmaxc(itype)   ! occupied core state
    2468             : 
    2469           0 :                   DO ll = ABS(l1 - l), l1 + l
    2470           0 :                      IF (ll < ABS(l - l2) .OR. ll > l + l2) CYCLE
    2471           0 :                      IF (MOD(l + l1 + ll, 2) /= 0) CYCLE
    2472           0 :                      IF (MOD(l + l2 + ll, 2) /= 0) CYCLE
    2473             : 
    2474           0 :                      DO m1 = -l1, l1
    2475           0 :                         m2 = m1
    2476           0 :                         IF (ABS(m2) > l2) CYCLE
    2477           0 :                         DO m = -l, l
    2478           0 :                            mm = m - m1
    2479           0 :                            IF (ABS(mm) > ll) CYCLE
    2480           0 :                            rdum = fpi_const/(2*ll + 1)*gaunt1(l, ll, l1, m, mm, m1, llmax)*gaunt1(l, ll, l2, m, mm, m2, llmax)
    2481             : 
    2482           0 :                            DO n = 1, nindxc(l, itype)
    2483           0 :                               DO n2 = 1, nindxc(l2, itype)
    2484             :                                  rprod(:) = (core1(:, n, l, itype)*core1(:, n2, l2, itype) &
    2485           0 :                                              + core2(:, n, l, itype)*core2(:, n2, l2, itype))/rmsh(:, itype)
    2486             : 
    2487             :                                  ! Initialization of n Summation
    2488           0 :                                  sum_primf = 0.0
    2489             :                                  ! Integration over r'
    2490           0 :                                  DO cn = 0, ncut ! Summation over Legendre polynomials
    2491           0 :                                     primf1 = 0.0
    2492           0 :                                     primf2 = 0.0
    2493             :                                     ! Calculate integral for 0 < r' < r
    2494           0 :                                     CALL primitivef(primf1, rprod(:)*rmsh(:, itype)**(ll + 2*cn + 1), rmsh, dx, jri, jmtd, itype, ntype)
    2495             :                                     ! Calculate integral for r < r' < R
    2496             :                                     !-itype is to enforce inward integration
    2497           0 :                                     CALL primitivef(primf2, d_ln(:, ll, cn)*rprod(:)/rmsh(:, itype)**(ll + 2*cn), rmsh, dx, jri, jmtd, -itype, ntype)
    2498             :                                     ! Multiplication with appropriate prefactors
    2499           0 :                                     primf1 = primf1/rmsh(:, itype)**(ll + 2*cn)*d_ln(:, ll, cn)
    2500           0 :                                     primf2 = primf2*rmsh(:, itype)**(ll + 2*cn + 1)
    2501             :                                     ! Summation over n
    2502           0 :                                     sum_primf = sum_primf + primf1 + primf2
    2503             :                                  END DO
    2504             : 
    2505           0 :                                  DO n1 = 1, nindxc(l1, itype)
    2506             : 
    2507             :                                     rprod(:) = (core1(:, n, l, itype)*core1(:, n1, l1, itype) &
    2508           0 :                                                 + core2(:, n, l, itype)*core2(:, n1, l1, itype))/rmsh(:, itype)
    2509             : 
    2510           0 :                                     integrand = rprod*sum_primf
    2511             : 
    2512           0 :                                     rdum1 = rdum*intgrf(integrand, jri, jmtd, rmsh, dx, ntype, itype, gridf)
    2513             : 
    2514           0 :                                     iatom = iatom0
    2515           0 :                                     DO ieq = 1, neq(itype)
    2516           0 :                                        iatom = iatom + 1
    2517           0 :                                        icst1 = point(n1, m1, l1, iatom)
    2518           0 :                                        icst2 = point(n2, m2, l2, iatom)
    2519           0 :                                        exch(icst1, icst2) = exch(icst1, icst2) + rdum1
    2520             :                                     END DO
    2521             :                                  END DO  !n1
    2522             : 
    2523             :                               END DO  !n2
    2524             :                            END DO  !n
    2525             : 
    2526             :                         END DO  !m
    2527             :                      END DO  !m1
    2528             : 
    2529             :                   END DO  !ll
    2530             : 
    2531             :                END DO  !l
    2532             :             END DO  !l2
    2533             :          END DO  !l1
    2534           0 :          iatom0 = iatom0 + neq(itype)
    2535             :       END DO  !itype
    2536             : 
    2537             : # ifdef CPP_INVERSION
    2538             :       CALL symmetrize(exch, ncstd, ncstd, 3, .FALSE., &
    2539             :                       ntype, ntype, neq, lmaxc, lmaxcd, &
    2540             :                       nindxc, natd, invsat, invsatnr)
    2541             :       IF (ANY(ABS(aimag(exch)) > 1E-6)) STOP 'exchange_cccc: exch possesses significant imaginary part'
    2542             : # endif
    2543             : !   DO icst = 1,ncstd
    2544             : !     IF ( irank == 0 ) &
    2545             : !       WRITE(6,'(    ''  ('',F5.3,'','',F5.3,'','',F5.3,'')'',I4,1X,F12.5)')bkpt,icst,REAL(exch(icst,icst))*(-27.211608)
    2546             : !   END DO
    2547             : 
    2548             :       ! add core exchange contributions to the te_hfex
    2549             : 
    2550           0 :       DO icst1 = 1, ncstd
    2551           0 :          te_hfex = te_hfex - a_ex*wtkpt(nk, nw)*exch(icst1, icst1)
    2552             :       END DO
    2553             : 
    2554           0 :    END SUBROUTINE exchange_ccccHSE
    2555             : 
    2556             :    ! The expansion coefficients of the attenuated Coulomb interaction in Legendre polynomials are calculated.
    2557             :    ! The attenuation is expressed in terms of a complementary error function
    2558             :    !
    2559             :    !                    -----              l+2n
    2560             :    !  erfc(w|r-r'|)      \                r           /           \
    2561             :    ! --------------- =    )    d  (wr') ---------  P (  cos(r,r')  )
    2562             :    !     |r-r'|          /      ln        l+2n+1    l \           /
    2563             :    !                    -----            r'
    2564             :    !                     l,n
    2565             :    ! where w is omega of the HSE functional, P_l are the Legendre polynomials and r' > r
    2566             :    !
    2567             :    ! The coefficients are given as
    2568             :    !                                    l-1
    2569             :    !                                2  -----   i+1  2i+1
    2570             :    !                       1      -x    \     2    x
    2571             :    !  d (x) = erfc(x) + -------  e       )   ------------
    2572             :    !   l,0                ____          /     (2i + 1)!!
    2573             :    !                     V Pi          -----
    2574             :    !                                    i=0
    2575             :    !
    2576             :    !                                       n-1
    2577             :    !               n     2                -----      i+1  l+i+1  2l+2n+2i+1
    2578             :    !           (-1)    -x        2l+1      \     (-1)    2      x
    2579             :    ! d (x)  = ------- e     -------------   )   ----------------------------
    2580             :    !  l,n>0     ____         n (2l+2n+1)   /      i! (n-i-1)! (2l+2i+1)!!
    2581             :    !           V Pi                       -----
    2582             :    !                                       i=0
    2583             :    !
    2584             :    ! Input:  rmsh - grid of radial position
    2585             :    !         lmax - maximum l value
    2586             :    !         ncut - maximum value for which coefficients are calculated
    2587             :    !         fac  - (optional) known factorials
    2588             :    ! Return: d_ln's for all r in rmsh, all 0 <= l <= lmax and all n <= ncut
    2589           0 :    FUNCTION calculate_coefficients(rmsh, lmax, ncut, fac) RESULT(d_ln)
    2590             : 
    2591             :       USE m_constants
    2592             : 
    2593             :       IMPLICIT NONE
    2594             : 
    2595             :       REAL, INTENT(IN)           :: rmsh(:)
    2596             :       INTEGER, INTENT(IN)        :: lmax, ncut
    2597             :       REAL, INTENT(IN), OPTIONAL :: fac(0:)
    2598             : 
    2599             :       INTEGER :: ci, cn, cl, fac_size
    2600             :       REAL    :: r1_spi, r2l_1, ri, r2l_2i_1, rn, r2l_2n_1
    2601           0 :       REAL    :: fn(0:MAX(0, ncut - 1))
    2602           0 :       REAL    :: rx(SIZE(rmsh)), x2(SIZE(rmsh)), x2n(SIZE(rmsh))
    2603           0 :       REAL    :: init(SIZE(rmsh)), addend(SIZE(rmsh))
    2604           0 :       REAL    :: erfc_x(SIZE(rmsh)), exp_x2(SIZE(rmsh))
    2605             :       REAL    :: d_ln(SIZE(rmsh), 0:lmax, 0:ncut)
    2606             : 
    2607             :       ! Initialize factorials
    2608           0 :       IF (PRESENT(fac)) THEN
    2609           0 :          fac_size = MIN(ncut, SIZE(fac))
    2610           0 :          fn(0:fac_size - 1) = fac(0:fac_size - 1)
    2611             :       ELSE
    2612           0 :          fac_size = 1
    2613           0 :          fn(0) = 1.0
    2614             :       END IF
    2615           0 :       DO cn = fac_size, ncut - 1
    2616           0 :          fn(cn) = fn(cn - 1)*cn
    2617             :       END DO
    2618             : 
    2619             :       ! Initialize 1 / sqrt(pi)
    2620             :       r1_spi = 1.0/SQRT(pi_const)
    2621             : 
    2622             :       ! Calculate x and x^2 and the functions erfc(x) and exp(x^2)
    2623           0 :       rx = omega_HSe*rmsh
    2624           0 :       x2 = rx**2
    2625           0 :       erfc_x = erfc(rx)
    2626           0 :       exp_x2 = EXP(-x2)
    2627             : 
    2628             :       ! Calculate the first coefficient
    2629             :       !                                    l-1
    2630             :       !                                2  -----   i+1  2i+1
    2631             :       !                       1      -x    \     2    x
    2632             :       !  d (x) = erfc(x) + -------  e       )   ------------
    2633             :       !   l,0                ____          /     (2i + 1)!!
    2634             :       !                     V Pi          -----
    2635             :       !                                    i=0
    2636             :       ! Initialize the first coefficient
    2637             :       r2l_1 = 1.0
    2638           0 :       init = 2.0*rx*r1_spi*exp_x2
    2639           0 :       addend = init
    2640           0 :       IF (lmax >= 0) d_ln(:, 0, 0) = erfc_x
    2641           0 :       IF (lmax >= 1) d_ln(:, 1, 0) = addend
    2642             :       ! Iterative calculation of the other coefficients
    2643           0 :       DO cl = 2, lmax
    2644           0 :          r2l_1 = r2l_1 + 2.0
    2645           0 :          addend = 2.0*addend*x2/r2l_1
    2646           0 :          d_ln(:, cl, 0) = d_ln(:, cl - 1, 0) + addend
    2647             :       END DO
    2648             :       ! Adding the erfc(x) term which is present in all coefficients
    2649           0 :       FORALL (cl=1:lmax)
    2650             :       d_ln(:, cl, 0) = erfc_x + d_ln(:, cl, 0)
    2651             :       END FORALL
    2652             : 
    2653             :       ! Calculate all other coefficients
    2654             :       !                                       n-1
    2655             :       !               n     2                -----      i+1  l+i+1  2l+2n+2i+1
    2656             :       !           (-1)    -x        2l+1      \     (-1)    2      x
    2657             :       ! d (x)  = ------- e     -------------   )   ----------------------------
    2658             :       !  l,n>0     ____         n (2l+2n+1)   /      i! (n-i-1)! (2l+2i+1)!!
    2659             :       !           V Pi                       -----
    2660             :       !                                       i=0
    2661             :       ! Initialize helper variables
    2662             :       r2l_1 = -1.0
    2663           0 :       DO cl = 0, lmax
    2664             :          ! Calculation of the l-dependent part
    2665           0 :          r2l_1 = r2l_1 + 2.0
    2666           0 :          addend = init*x2               ! addend ~ 2^(l+1) x^(2l+1) / (2l-1)!!
    2667           0 :          init = 2.0*addend/r2l_1
    2668             :          ! Summation over all i
    2669             :          ! i = 0 term
    2670           0 :          DO cn = 1, ncut
    2671           0 :             d_ln(:, cl, cn) = addend/fn(cn - 1)
    2672             :          END DO
    2673             :          ! higher values of i
    2674             :          ri = 0.0
    2675             :          r2l_2i_1 = r2l_1
    2676           0 :          DO ci = 1, ncut - 1
    2677           0 :             ri = ri + 1.0
    2678           0 :             r2l_2i_1 = r2l_2i_1 + 2.0
    2679           0 :             addend = -2.0*addend*x2/(ri*r2l_2i_1)       ! addend ~ (-1)^i+1 2^(l+i+1) x^(2l+2i+1) (2l+1) / i!(2l+2i+1)!!
    2680           0 :             DO cn = ci + 1, ncut
    2681           0 :                d_ln(:, cl, cn) = d_ln(:, cl, cn) + addend/fn(cn - ci - 1)
    2682             :             END DO
    2683             :          END DO
    2684             :          r2l_2n_1 = r2l_1
    2685             :          ! Divide by n and l dependent part
    2686           0 :          DO cn = 1, ncut
    2687           0 :             r2l_2n_1 = r2l_2n_1 + 2.0
    2688           0 :             d_ln(:, cl, cn) = d_ln(:, cl, cn)/r2l_2n_1
    2689             :          END DO
    2690             :       END DO
    2691             :       ! Resolve only-n dependent part
    2692             :       rn = 1.0
    2693           0 :       x2n = x2
    2694           0 :       IF (ncut >= 1) THEN
    2695           0 :          FORALL (cl=0:lmax)
    2696             :          d_ln(:, cl, 1) = d_ln(:, cl, 1)*x2
    2697             :          END FORALL
    2698             :       END IF
    2699           0 :       DO cn = 2, ncut
    2700           0 :          rn = rn + 1.0
    2701           0 :          x2n = -x2n*x2
    2702           0 :          FORALL (cl=0:lmax)
    2703             :          d_ln(:, cl, cn) = d_ln(:, cl, cn)*x2n/rn
    2704             :          END FORALL
    2705             :       END DO
    2706             : 
    2707           0 :    END FUNCTION calculate_coefficients
    2708             : 
    2709             : END MODULE m_hsefunctional

Generated by: LCOV version 1.13