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