Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculates short range exchange part for wPBE functional and averaged
10 : !> PBE exchange-hole functional (omega = 0.0 )
11 : !> \par History
12 : !> Manuel Guidon (05.2007) : initial version
13 : !> \author Manuel Guidon (05.2007)
14 : ! **************************************************************************************************
15 : MODULE xc_xwpbe
16 : USE bibliography, ONLY: Heyd2004, &
17 : cite_reference
18 : USE input_section_types, ONLY: section_vals_type, &
19 : section_vals_val_get
20 : USE kinds, ONLY: dp
21 : USE mathconstants, ONLY: pi, &
22 : rootpi
23 : USE xc_derivative_desc, ONLY: deriv_norm_drho, &
24 : deriv_norm_drhoa, &
25 : deriv_norm_drhob, &
26 : deriv_rho, &
27 : deriv_rhoa, &
28 : deriv_rhob
29 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
30 : xc_dset_get_derivative
31 : USE xc_derivative_types, ONLY: xc_derivative_get, &
32 : xc_derivative_type
33 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
34 : USE xc_rho_set_types, ONLY: xc_rho_set_get, &
35 : xc_rho_set_type
36 : #include "../base/base_uses.f90"
37 : #:include '../common/expint.fypp'
38 :
39 : IMPLICIT NONE
40 :
41 : PRIVATE
42 :
43 : ! *** Global parameters ***
44 :
45 : PUBLIC :: xwpbe_lda_info, xwpbe_lda_eval, xwpbe_lsd_info, &
46 : xwpbe_lsd_eval
47 :
48 : REAL(KIND=dp), PARAMETER :: alpha1 = -1.128223946706117_dp, &
49 : alpha2 = 1.452736265762971_dp, &
50 : alpha3 = -1.243162299390327_dp, &
51 : alpha4 = 0.971824836115601_dp, &
52 : alpha5 = -0.568861079687373_dp, &
53 : alpha6 = 0.246880514820192_dp, &
54 : alpha7 = -0.065032363850763_dp, &
55 : alpha8 = 0.008401793031216_dp
56 : REAL(KIND=dp), PARAMETER :: beta = 1.455915450052607_dp, &
57 : beta2 = 2.0_dp
58 : REAL(KIND=dp), PARAMETER :: a1 = 0.00979681_dp, &
59 : a2 = 0.04108340_dp, &
60 : a3 = 0.18744000_dp, &
61 : a4 = 0.00120824_dp, &
62 : a5 = 0.0347188_dp
63 : REAL(KIND=dp), PARAMETER :: A = 1.0161144_dp, &
64 : B = -0.37170836_dp, &
65 : C = -0.077215461_dp, &
66 : DD = 0.57786348_dp, &
67 : E = -0.051955731_dp, &
68 : F1 = 0.47965830_dp, &
69 : F2 = 6.4753871_dp, &
70 : clda = -0.73855876638202240588423_dp
71 : REAL(KIND=dp), PARAMETER :: expcutoff = 700.0_dp, &
72 : exei1 = 4.0364_dp, &
73 : exei2 = 1.15198_dp, &
74 : exei3 = 5.03627_dp, &
75 : exei4 = 4.19160_dp
76 : REAL(KIND=dp), PARAMETER :: smax = 8.572844_dp, &
77 : sconst = 18.79622316_dp, &
78 : scutoff = 8.3_dp
79 : REAL(KIND=dp), PARAMETER :: gcutoff = 0.08_dp, &
80 : g1 = -0.02628417880_dp/E, &
81 : g2 = -0.07117647788_dp/E, &
82 : g3 = 0.08534541323_dp/E, &
83 : g4 = 0.0_dp
84 : REAL(KIND=dp), PARAMETER :: wcutoff = 14.0_dp
85 : REAL(KIND=dp), PARAMETER :: f12 = 0.5_dp, f14 = 0.25_dp, f158 = 15.0_dp/8.0_dp, &
86 : f1516 = 15.0_dp/16.0_dp, f24364 = 243.0_dp/64.0_dp, &
87 : f2716 = 27.0_dp/16.0_dp, f2732 = 27.0_dp/32.0_dp, &
88 : f34 = 0.75_dp, f32 = 1.5_dp, f38 = 0.375_dp, f68 = 0.75_dp, &
89 : f6561512 = 6561.0_dp/512.0_dp, f8132 = 81.0_dp/32.0_dp, &
90 : f8164 = 81.0_dp/64.0_dp, f729128 = 729.0_dp/128.0_dp, &
91 : f52 = 2.5_dp, f94 = 9.0_dp/4.0_dp, f916 = 9.0_dp/16.0_dp, &
92 : f89 = 8.0_dp/9.0_dp, f2187256 = 2187.0_dp/256.0_dp, &
93 : r1 = 1.0_dp, f98 = 9.0_dp/8.0_dp, r15 = 15.0_dp, &
94 : r3 = 3.0_dp, r4 = 4.0_dp, r16 = 16.0_dp, r8 = 8.0_dp, &
95 : r6 = 6.0_dp, r2 = 2.0_dp
96 :
97 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xwpbe'
98 :
99 : CONTAINS
100 :
101 : ! **************************************************************************************************
102 : !> \brief return various information on the functional
103 : !> \param reference string with the reference of the actual functional
104 : !> \param shortform string with the shortform of the functional name
105 : !> \param needs the components needed by this functional are set to
106 : !> true (does not set the unneeded components to false)
107 : !> \param max_deriv ...
108 : !> \par History
109 : !> 05.2007 created [Manuel Guidon]
110 : !> \author Manuel Guidon
111 : ! **************************************************************************************************
112 2963 : SUBROUTINE xwpbe_lda_info(reference, shortform, needs, max_deriv)
113 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
114 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
115 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
116 :
117 2963 : IF (PRESENT(reference)) THEN
118 9 : reference = "Jochen Heyd and Gustavo E. Scuseria, J. Chem. Phys., 120, 7274 {LDA version}"
119 : END IF
120 2963 : IF (PRESENT(shortform)) THEN
121 9 : shortform = "shortrange part of PBE exchange {LDA}"
122 : END IF
123 2963 : IF (PRESENT(needs)) THEN
124 2954 : needs%rho = .TRUE.
125 2954 : needs%norm_drho = .TRUE.
126 : END IF
127 : ! deriv > 1 are not correct
128 2963 : IF (PRESENT(max_deriv)) max_deriv = 1
129 2963 : END SUBROUTINE xwpbe_lda_info
130 :
131 : ! **************************************************************************************************
132 : !> \brief evaluates the screened hole averaged PBE exchange functional for lda
133 : !> \param rho_set the density where you want to evaluate the functional
134 : !> \param deriv_set place where to store the functional derivatives (they are
135 : !> added to the derivatives)
136 : !> \param order maximum degree of the derivative that should be evaluated
137 : !> \param xwpbe_params input parameters (scaling,omega)
138 : !> \par History
139 : !> 05.2007 created [Manuel Guidon]
140 : !> \author Manuel Guidon
141 : !> \note
142 : !> The current version provides code for derivatives up to second order.
143 : !> Using the maple sheet in cp2k/doc it is straightforward to produce routines
144 : !> for higher derivatives.
145 : ! **************************************************************************************************
146 18840 : SUBROUTINE xwpbe_lda_eval(rho_set, deriv_set, order, xwpbe_params)
147 :
148 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
149 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
150 : INTEGER, INTENT(IN) :: order
151 : TYPE(section_vals_type), POINTER :: xwpbe_params
152 :
153 : CHARACTER(len=*), PARAMETER :: routineN = 'xwpbe_lda_eval'
154 :
155 : INTEGER :: handle, npoints
156 : INTEGER, DIMENSION(2, 3) :: bo
157 : REAL(kind=dp) :: epsilon_norm_drho, epsilon_rho, omega, &
158 : sx, sx0
159 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
160 3768 : POINTER :: dummy, e_0, e_ndrho, e_rho, norm_drho, &
161 3768 : rho
162 : TYPE(xc_derivative_type), POINTER :: deriv
163 :
164 3768 : CALL timeset(routineN, handle)
165 :
166 3768 : CALL cite_reference(Heyd2004)
167 :
168 3768 : CALL section_vals_val_get(xwpbe_params, "SCALE_X", r_val=sx)
169 3768 : CALL section_vals_val_get(xwpbe_params, "SCALE_X0", r_val=sx0)
170 3768 : CALL section_vals_val_get(xwpbe_params, "OMEGA", r_val=omega)
171 :
172 : CALL xc_rho_set_get(rho_set, rho=rho, &
173 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
174 3768 : drho_cutoff=epsilon_norm_drho)
175 3768 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
176 :
177 3768 : dummy => rho
178 :
179 3768 : e_0 => dummy
180 3768 : e_rho => dummy
181 3768 : e_ndrho => dummy
182 :
183 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
184 3768 : allocate_deriv=.TRUE.)
185 3768 : CALL xc_derivative_get(deriv, deriv_data=e_0)
186 3768 : IF (order >= 1) THEN
187 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
188 3350 : allocate_deriv=.TRUE.)
189 3350 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
190 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
191 3350 : allocate_deriv=.TRUE.)
192 3350 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
193 : END IF
194 3768 : IF (order > 1) THEN
195 0 : CPABORT("derivatives bigger than 1 not implemented")
196 : END IF
197 :
198 : !$OMP PARALLEL DEFAULT(NONE) &
199 : !$OMP SHARED(npoints, order, rho, norm_drho, e_0, e_rho, e_ndrho) &
200 3768 : !$OMP SHARED(epsilon_rho, sx, sx0, omega)
201 :
202 : CALL xwpbe_lda_calc(npoints, order, rho=rho, norm_drho=norm_drho, &
203 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
204 : epsilon_rho=epsilon_rho, sx=sx, sx0=sx0, omega=omega)
205 :
206 : !$OMP END PARALLEL
207 :
208 3768 : CALL timestop(handle)
209 :
210 3768 : END SUBROUTINE xwpbe_lda_eval
211 :
212 : ! **************************************************************************************************
213 : !> \brief evaluates the screened hole averaged PBE exchange functional for lda
214 : !> \param npoints ...
215 : !> \param order maximum degree of the derivative that should be evaluated
216 : !> \param rho , ndrho: density and norm of the density gradient
217 : !> \param norm_drho ...
218 : !> \param e_0 ...
219 : !> \param e_rho ...
220 : !> \param e_ndrho ...
221 : !> \param epsilon_rho ...
222 : !> \param sx , sx0: scaling factor for omega!=0 and omega=0
223 : !> \param sx0 ...
224 : !> \param omega screening parameter
225 : !> \par History
226 : !> 05.2007 created [Manuel Guidon]
227 : !> \author Manuel Guidon
228 : !> \note
229 : !> In order to avoid numerical instabilities, this routine calls different
230 : !> subroutines. There are 4 routines for the case omega!=0 and 2 routines
231 : !> for omega=0.
232 : ! **************************************************************************************************
233 3768 : SUBROUTINE xwpbe_lda_calc(npoints, order, rho, norm_drho, e_0, e_rho, e_ndrho, &
234 : epsilon_rho, sx, sx0, omega)
235 :
236 : INTEGER, INTENT(in) :: npoints, order
237 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, norm_drho, e_0, e_rho, e_ndrho
238 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, sx0, omega
239 :
240 : INTEGER :: ip
241 : REAL(dp) :: gradient_before, gradient_increment, &
242 : my_ndrho, my_rho, pullback
243 : REAL(KIND=dp) :: ss, ss2, sscale, t1, t2, t3, t4, t5, t6, &
244 : t7, t8, ww
245 :
246 : !$OMP DO
247 :
248 : DO ip = 1, npoints
249 83977303 : my_rho = MAX(rho(ip), 0.0_dp)
250 83977303 : IF (my_rho > epsilon_rho) THEN
251 69821467 : gradient_before = e_ndrho(ip)
252 69821467 : my_ndrho = MAX(norm_drho(ip), 0.0_dp)
253 :
254 : !Do some precalculation in order to catch the correct branch afterwards
255 69821467 : sscale = 1.0_dp
256 69821467 : t1 = pi**2
257 69821467 : t2 = t1*my_rho
258 69821467 : t3 = t2**(0.1e1_dp/0.3e1_dp)
259 69821467 : t4 = 0.1e1_dp/t3
260 69821467 : t5 = omega*t4
261 69821467 : ww = 0.6933612743506347048433524e0_dp*t5
262 69821467 : t6 = my_ndrho*t4
263 69821467 : t7 = 0.1e1_dp/my_rho
264 69821467 : t8 = t7*sscale
265 69821467 : ss = 0.3466806371753173524216762e0_dp*t6*t8
266 69821467 : IF (ss > scutoff) THEN
267 46741279 : ss2 = ss*ss
268 46741279 : sscale = (smax*ss2 - sconst)/(ss2*ss)
269 : END IF
270 :
271 69821467 : IF (sx0 /= 0.0_dp) THEN
272 : !original PBE hole
273 59543031 : IF (ss*sscale > gcutoff) THEN
274 : CALL xwpbe_lda_calc_0(e_0(ip), e_rho(ip), e_ndrho(ip), my_rho, &
275 59263608 : my_ndrho, sscale, sx0, order)
276 : ELSE
277 : CALL xwpbe_lda_calc_01(e_0(ip), e_rho(ip), e_ndrho(ip), my_rho, &
278 279423 : my_ndrho, sscale, sx0, order)
279 : END IF
280 : END IF
281 :
282 69821467 : IF (sx /= 0.0_dp) THEN
283 24513811 : IF (omega == 0.0_dp) THEN
284 28 : IF (ss*sscale > gcutoff) THEN
285 : CALL xwpbe_lda_calc_0(e_0(ip), e_rho(ip), e_ndrho(ip), my_rho, &
286 20 : my_ndrho, sscale, sx, order)
287 : ELSE
288 : CALL xwpbe_lda_calc_01(e_0(ip), e_rho(ip), e_ndrho(ip), my_rho, &
289 8 : my_ndrho, sscale, sx, order)
290 : END IF
291 24513783 : ELSE IF (ww < wcutoff .AND. ss*sscale > gcutoff) THEN
292 : CALL xwpbe_lda_calc_1(e_0(ip), e_rho(ip), e_ndrho(ip), my_rho, &
293 17846301 : my_ndrho, omega, sscale, sx, order)
294 233512 : ELSE IF (ww < wcutoff .AND. ss*sscale <= gcutoff) THEN
295 : CALL xwpbe_lda_calc_2(e_0(ip), e_rho(ip), e_ndrho(ip), my_rho, &
296 233512 : my_ndrho, omega, sscale, sx, order)
297 6433970 : ELSE IF (ww >= wcutoff .AND. ss*sscale > gcutoff) THEN
298 : CALL xwpbe_lda_calc_3(e_0(ip), e_rho(ip), e_ndrho(ip), my_rho, &
299 6433848 : my_ndrho, omega, sscale, sx, order)
300 : ELSE
301 : CALL xwpbe_lda_calc_4(e_0(ip), e_rho(ip), e_ndrho(ip), my_rho, &
302 122 : my_ndrho, omega, sscale, sx, order)
303 : END IF
304 : END IF
305 69821467 : IF (order >= 1 .AND. ss > scutoff) THEN
306 : ! Pull back g_eff = g*f(s), s = g/(2*k_F*rho), through the large-s scaling.
307 : ! The generated expressions differentiate at fixed f(s).
308 31462880 : pullback = (-smax*ss2 + 3.0_dp*sconst)/(smax*ss2 - sconst)
309 31462880 : gradient_increment = e_ndrho(ip) - gradient_before
310 31462880 : e_rho(ip) = e_rho(ip) - gradient_increment*(4.0_dp*my_ndrho/(3.0_dp*my_rho))*pullback
311 31462880 : e_ndrho(ip) = gradient_before + gradient_increment*(1.0_dp + pullback)
312 : END IF
313 : END IF
314 : END DO
315 :
316 : !$OMP END DO
317 :
318 3768 : END SUBROUTINE xwpbe_lda_calc
319 :
320 : ! **************************************************************************************************
321 : !> \brief Evaluates the screened hole averaged PBE exchange functional for lda
322 : !> \param e_0 ...
323 : !> \param e_rho ...
324 : !> \param e_ndrho ...
325 : !> \param rho , ndrho: density and norm of the density gradient
326 : !> \param ndrho ...
327 : !> \param sscale scaling factor to enforce Lieb-Oxford bound
328 : !> \param sx0 scaling factor
329 : !> \param order maximum degree of the derivative that should be evaluated
330 : !> \par History
331 : !> 05.2007 created [Manuel Guidon]
332 : !> \author Manuel Guidon
333 : !> \note
334 : !> This routine evaluates the exact functional for omega=0.
335 : ! **************************************************************************************************
336 87518835 : SUBROUTINE xwpbe_lda_calc_0(e_0, e_rho, e_ndrho, rho, ndrho, sscale, sx0, order)
337 : REAL(KIND=dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho
338 : REAL(KIND=dp), INTENT(IN) :: rho, ndrho, sscale, sx0
339 : INTEGER, INTENT(IN) :: order
340 :
341 : REAL(KIND=dp) :: dQndrho, dQrho, Q, t1, t10, t100, t103, t105, t106, t109, t11, t110, t111, &
342 : t112, t114, t116, t117, t118, t119, t12, t120, t121, t122, t123, t125, t126, t127, t13, &
343 : t130, t134, t135, t136, t138, t139, t14, t142, t143, t145, t146, t147, t148, t149, t15, &
344 : t150, t152, t156, t158, t159, t16, t160, t163, t164, t165, t166, t169, t170, t172, t173, &
345 : t174, t175, t176, t18, t180, t183, t184, t185, t187, t188, t19, t190, t191, t192, t193, &
346 : t199, t2, t20, t200, t201, t202, t203, t207, t209, t21, t215, t218, t219, t22, t220, &
347 : t222, t223, t224, t225, t227, t228, t231, t233, t236, t237, t24
348 : REAL(KIND=dp) :: t240, t241, t242, t245, t246, t247, t249, t25, t250, t253, t254, t258, t26, &
349 : t261, t262, t263, t265, t266, t269, t272, t274, t275, t278, t28, t281, t282, t284, t285, &
350 : t287, t288, t29, t290, t291, t294, t297, t3, t300, t301, t303, t304, t306, t307, t31, &
351 : t314, t32, t321, t323, t326, t327, t33, t330, t334, t335, t336, t337, t338, t339, t34, &
352 : t340, t342, t343, t344, t345, t347, t355, t356, t357, t358, t359, t36, t361, t362, t363, &
353 : t367, t368, t37, t372, t374, t376, t377, t381, t383, t384, t387, t388, t39, t390, t391, &
354 : t394, t397, t398, t4, t40, t400, t401, t402, t404, t405, t406
355 : REAL(KIND=dp) :: t408, t409, t410, t411, t412, t413, t414, t415, t416, t417, t418, t42, &
356 : t423, t425, t426, t429, t430, t432, t433, t435, t44, t440, t443, t444, t446, t447, t449, &
357 : t46, t463, t465, t47, t471, t472, t475, t476, t478, t479, t480, t484, t485, t489, t491, &
358 : t495, t497, t5, t500, t501, t504, t505, t507, t508, t510, t511, t513, t514, t515, t52, &
359 : t520, t53, t54, t56, t58, t59, t6, t61, t63, t65, t66, t67, t68, t69, t7, t70, t71, t72, &
360 : t73, t74, t76, t77, t78, t79, t8, t80, t81, t83, t84, t86, t88, t89, t9, t91, t92, t93, &
361 : t94, t95, t97, t98
362 :
363 87518835 : t1 = ndrho**2
364 87518835 : t2 = a1*t1
365 87518835 : t3 = r2**2
366 87518835 : t4 = 0.1e1_dp/t3
367 87518835 : t5 = t2*t4
368 87518835 : t6 = pi**2
369 87518835 : t7 = r3*t6
370 87518835 : t8 = t7*rho
371 87518835 : t9 = t8**(0.1e1_dp/0.3e1_dp)
372 87518835 : t10 = t9**2
373 87518835 : t11 = 0.1e1_dp/t10
374 87518835 : t12 = rho**2
375 87518835 : t13 = 0.1e1_dp/t12
376 87518835 : t14 = t11*t13
377 87518835 : t15 = sscale**2
378 87518835 : t16 = t14*t15
379 87518835 : t18 = t1**2
380 87518835 : t19 = a2*t18
381 87518835 : t20 = t3**2
382 87518835 : t21 = 0.1e1_dp/t20
383 87518835 : t22 = t19*t21
384 87518835 : t24 = 0.1e1_dp/t9/t8
385 87518835 : t25 = t12**2
386 87518835 : t26 = 0.1e1_dp/t25
387 87518835 : t28 = t15**2
388 87518835 : t29 = t24*t26*t28
389 87518835 : t31 = t5*t16 + t22*t29
390 87518835 : t32 = f94*t31
391 87518835 : t33 = a3*t18
392 87518835 : t34 = t33*t21
393 87518835 : t36 = t18*ndrho
394 87518835 : t37 = a4*t36
395 87518835 : t39 = 0.1e1_dp/t20/r2
396 87518835 : t40 = t37*t39
397 87518835 : t42 = 0.1e1_dp/t10/t8
398 87518835 : t44 = 0.1e1_dp/t25/rho
399 87518835 : t46 = t28*sscale
400 87518835 : t47 = t42*t44*t46
401 87518835 : t52 = 0.1e1_dp/t20/t3
402 87518835 : t53 = a5*t18*t1*t52
403 87518835 : t54 = r3**2
404 87518835 : t56 = t6**2
405 87518835 : t58 = 0.1e1_dp/t54/t56
406 87518835 : t59 = t25**2
407 87518835 : t61 = t28*t15
408 87518835 : t63 = t58/t59*t61
409 87518835 : t65 = r1 + t34*t29 + t40*t47 + t53*t63
410 87518835 : t66 = 0.1e1_dp/t65
411 87518835 : t67 = t66*t1
412 87518835 : t68 = t32*t67
413 87518835 : t69 = t4*t11
414 87518835 : t70 = t13*t15
415 87518835 : t71 = 0.1e1_dp/A
416 87518835 : t72 = t70*t71
417 87518835 : t73 = t69*t72
418 87518835 : Q = t68*t73
419 87518835 : t74 = rho**(0.1e1_dp/0.3e1_dp)
420 87518835 : t76 = t74*rho*f89
421 87518835 : t77 = B*f12
422 87518835 : t78 = t1*t4
423 87518835 : t79 = t78*t11
424 87518835 : t80 = t31*t66
425 87518835 : t81 = t70*t80
426 87518835 : t83 = t79*t81 + DD
427 87518835 : t84 = 0.1e1_dp/t83
428 87518835 : t86 = F2*t31
429 87518835 : t88 = F1 + t86*t66
430 87518835 : t89 = t70*t88
431 87518835 : t91 = t79*t89 + r1
432 87518835 : t92 = f12*t91
433 87518835 : t93 = t83**2
434 87518835 : t94 = 0.1e1_dp/t93
435 87518835 : t95 = C*t94
436 87518835 : t97 = f34*pi
437 87518835 : t98 = rootpi
438 87518835 : t100 = r6*C
439 87518835 : t103 = r4*B
440 87518835 : t105 = r8*A
441 87518835 : t106 = t93*t83
442 87518835 : t109 = t98*(r15*E + t100*t91*t83 + t103*t93 + t105*t106)
443 87518835 : t110 = 0.1e1_dp/r16
444 87518835 : t111 = SQRT(t83)
445 87518835 : t112 = t111*t106
446 87518835 : t114 = t110/t112
447 87518835 : t116 = SQRT(A)
448 87518835 : t117 = EXP(Q)
449 87518835 : t118 = t116*t117
450 87518835 : t119 = f32*ndrho
451 87518835 : t120 = 0.1e1_dp/r2
452 87518835 : t121 = t119*t120
453 87518835 : t122 = 0.1e1_dp/t9
454 87518835 : t123 = 0.1e1_dp/rho
455 87518835 : t125 = t80*t71
456 87518835 : t126 = SQRT(t125)
457 87518835 : t127 = sscale*t126
458 87518835 : t130 = erfc(t121*t122*t123*t127)
459 87518835 : t134 = 0.1e1_dp/f1516
460 87518835 : t135 = (t97 + t109*t114 - t97*t118*t130)*t134
461 87518835 : t136 = 0.1e1_dp/t98
462 87518835 : t138 = 0.1e1_dp/E
463 87518835 : t139 = t136*t112*t138
464 87518835 : t142 = (-t135*t139 + r1)*E
465 87518835 : t143 = 0.1e1_dp/t106
466 87518835 : t145 = f12*A
467 87518835 : t146 = exei(Q, t117)
468 87518835 : t147 = t78*t14
469 87518835 : t148 = t15*t31
470 87518835 : t149 = t66*t84
471 87518835 : t150 = t148*t149
472 87518835 : t152 = LOG(t147*t150)
473 : t156 = (t77*t84 + t92*t95 + t142*t143 + t145*(t146 + t152)) &
474 87518835 : *Clda
475 87518835 : e_0 = e_0 + (-t76*t156)*sx0
476 87518835 : IF (order >= 1) THEN
477 64894530 : t158 = t4*t42
478 64894530 : t159 = t2*t158
479 64894530 : t160 = t70*t7
480 64894530 : t163 = t12*rho
481 64894530 : t164 = 0.1e1_dp/t163
482 64894530 : t165 = t11*t164
483 64894530 : t166 = t165*t15
484 64894530 : t169 = t54*t56
485 64894530 : t170 = t169*t12
486 64894530 : t172 = 0.1e1_dp/t9/t170
487 64894530 : t173 = t21*t172
488 64894530 : t174 = t19*t173
489 64894530 : t175 = t26*t28
490 64894530 : t176 = t175*t7
491 64894530 : t180 = t24*t44*t28
492 : t183 = -0.2e1_dp/0.3e1_dp*t159*t160 - (2._dp*t5*t166) - 0.4e1_dp/ &
493 64894530 : 0.3e1_dp*t174*t176 - (4._dp*t22*t180)
494 64894530 : t184 = f94*t183
495 64894530 : t185 = t184*t67
496 64894530 : t187 = t65**2
497 64894530 : t188 = 0.1e1_dp/t187
498 64894530 : t190 = t188*t1*t4
499 64894530 : t191 = t32*t190
500 64894530 : t192 = t15*t71
501 64894530 : t193 = t33*t173
502 64894530 : t199 = 0.1e1_dp/t10/t170
503 64894530 : t200 = t39*t199
504 64894530 : t201 = t37*t200
505 64894530 : t202 = t44*t46
506 64894530 : t203 = t202*t7
507 64894530 : t207 = 0.1e1_dp/t25/t12
508 64894530 : t209 = t42*t207*t46
509 64894530 : t215 = t58/t59/rho*t61
510 : t218 = -0.4e1_dp/0.3e1_dp*t193*t176 - (4._dp*t34*t180) - 0.5e1_dp &
511 64894530 : /0.3e1_dp*t201*t203 - (5._dp*t40*t209) - (8._dp*t53*t215)
512 64894530 : t219 = t192*t218
513 64894530 : t220 = t14*t219
514 64894530 : t222 = t67*t4
515 64894530 : t223 = t32*t222
516 64894530 : t224 = t42*t13
517 64894530 : t225 = t224*t15
518 64894530 : t227 = t71*r3*t6
519 64894530 : t228 = t225*t227
520 64894530 : t231 = t164*t15
521 64894530 : t233 = t69*t231*t71
522 : dQrho = t185*t73 - t191*t220 - 0.2e1_dp/0.3e1_dp*t223*t228 - (2._dp &
523 64894530 : *t68*t233)
524 64894530 : t236 = a1*ndrho
525 64894530 : t237 = t236*t4
526 64894530 : t240 = t1*ndrho
527 64894530 : t241 = a2*t240
528 64894530 : t242 = t241*t21
529 64894530 : t245 = 2._dp*t237*t16 + 4._dp*t242*t29
530 64894530 : t246 = f94*t245
531 64894530 : t247 = t246*t67
532 64894530 : t249 = a3*t240
533 64894530 : t250 = t249*t21
534 64894530 : t253 = a4*t18
535 64894530 : t254 = t253*t39
536 64894530 : t258 = a5*t36*t52
537 64894530 : t261 = 4._dp*t250*t29 + 5._dp*t254*t47 + 6._dp*t258*t63
538 64894530 : t262 = t192*t261
539 64894530 : t263 = t14*t262
540 64894530 : t265 = t66*ndrho
541 64894530 : t266 = t32*t265
542 64894530 : dQndrho = t247*t73 - t191*t263 + 2._dp*t266*t73
543 64894530 : t269 = t74*f89
544 64894530 : t272 = t78*t224
545 64894530 : t274 = t66*r3*t6
546 64894530 : t275 = t148*t274
547 64894530 : t278 = t231*t80
548 64894530 : t281 = t183*t66
549 64894530 : t282 = t70*t281
550 64894530 : t284 = t188*t218
551 64894530 : t285 = t148*t284
552 : t287 = -0.2e1_dp/0.3e1_dp*t272*t275 - (2._dp*t79*t278) + (t79 &
553 64894530 : *t282) - t147*t285
554 64894530 : t288 = t94*t287
555 64894530 : t290 = t15*t88
556 64894530 : t291 = t290*t7
557 64894530 : t294 = t231*t88
558 64894530 : t297 = F2*t183
559 64894530 : t300 = t297*t66 - t86*t284
560 64894530 : t301 = t70*t300
561 : t303 = -0.2e1_dp/0.3e1_dp*t272*t291 - (2._dp*t79*t294) + (t79 &
562 64894530 : *t301)
563 64894530 : t304 = f12*t303
564 64894530 : t306 = C*t143
565 64894530 : t307 = t306*t287
566 64894530 : t314 = t83*t287
567 : t321 = t98*(t100*t303*t83 + t100*t91*t287 + 2._dp*t103*t314 &
568 64894530 : + 3._dp*t105*t93*t287)
569 64894530 : t323 = t93**2
570 64894530 : t326 = t110/t111/t323
571 64894530 : t327 = t326*t287
572 64894530 : t330 = t97*t116
573 64894530 : t334 = rootpi
574 64894530 : t335 = 0.1e1_dp/t334
575 64894530 : t336 = t117*t335
576 64894530 : t337 = f32**2
577 64894530 : t338 = t337*t1
578 64894530 : t339 = t338*t69
579 64894530 : t340 = t70*t125
580 64894530 : t342 = 1.0_dp/t117
581 64894530 : t343 = t120*t24
582 64894530 : t344 = t119*t343
583 64894530 : t345 = t123*sscale
584 64894530 : t347 = t126*r3*t6
585 64894530 : t355 = t119*t120*t122
586 64894530 : t356 = 0.1e1_dp/t126
587 64894530 : t357 = t281*t71
588 64894530 : t358 = t31*t188
589 64894530 : t359 = t71*t218
590 64894530 : t361 = t357 - t358*t359
591 64894530 : t362 = t356*t361
592 64894530 : t363 = t345*t362
593 : t367 = t342*(-t344*t345*t347/0.3e1_dp - t121*t122*t13*t127 &
594 64894530 : + t355*t363/0.2e1_dp)
595 64894530 : t368 = t336*t367
596 : t372 = (t321*t114 - 0.7e1_dp/0.2e1_dp*t109*t327 - (t330*dQrho &
597 64894530 : *t117*t130) + (2._dp*t330*t368))*t134
598 64894530 : t374 = t135*t136
599 64894530 : t376 = t111*t93*t138
600 64894530 : t377 = t376*t287
601 64894530 : t381 = (-t372*t139 - 0.7e1_dp/0.2e1_dp*t374*t377)*E
602 64894530 : t383 = 0.1e1_dp/t323
603 64894530 : t384 = t383*t287
604 64894530 : t387 = dexeirho(Q, dQrho, t146)
605 64894530 : t388 = t78*t225
606 64894530 : t390 = t84*r3*t6
607 64894530 : t391 = t80*t390
608 64894530 : t394 = t78*t165
609 64894530 : t397 = t15*t183
610 64894530 : t398 = t397*t149
611 64894530 : t400 = t188*t84
612 64894530 : t401 = t400*t218
613 64894530 : t402 = t148*t401
614 64894530 : t404 = t66*t94
615 64894530 : t405 = t404*t287
616 64894530 : t406 = t148*t405
617 : t408 = -0.2e1_dp/0.3e1_dp*t388*t391 - (2._dp*t394*t150) + t147 &
618 64894530 : *t398 - t147*t402 - t147*t406
619 64894530 : t409 = 0.1e1_dp/t1
620 64894530 : t410 = t408*t409
621 64894530 : t411 = t3*t10
622 64894530 : t412 = t410*t411
623 64894530 : t413 = 0.1e1_dp/t15
624 64894530 : t414 = t12*t413
625 64894530 : t415 = 0.1e1_dp/t31
626 64894530 : t416 = t415*t65
627 64894530 : t417 = t416*t83
628 64894530 : t418 = t414*t417
629 : t423 = (-t77*t288 + t304*t95 - 2._dp*t92*t307 + t381*t143 - 3._dp &
630 64894530 : *t142*t384 + t145*(t387 + t412*t418))*Clda
631 64894530 : e_rho = e_rho + (-0.4e1_dp/0.3e1_dp*t269*t156 - t76*t423)*sx0
632 64894530 : t425 = ndrho*t4
633 64894530 : t426 = t425*t11
634 64894530 : t429 = t245*t66
635 64894530 : t430 = t70*t429
636 64894530 : t432 = t188*t261
637 64894530 : t433 = t148*t432
638 64894530 : t435 = 2._dp*t426*t81 + t79*t430 - t147*t433
639 64894530 : t440 = F2*t245
640 64894530 : t443 = t440*t66 - t86*t432
641 64894530 : t444 = t70*t443
642 64894530 : t446 = 2._dp*t426*t89 + t79*t444
643 64894530 : t447 = f12*t446
644 64894530 : t449 = t306*t435
645 : t463 = t98*(t100*t446*t83 + t100*t91*t435 + 2._dp*t103*t83 &
646 64894530 : *t435 + 3._dp*t105*t93*t435)
647 64894530 : t465 = t326*t435
648 64894530 : t471 = f32*t120
649 64894530 : t472 = t471*t122
650 64894530 : t475 = t429*t71
651 64894530 : t476 = t71*t261
652 64894530 : t478 = t475 - t358*t476
653 64894530 : t479 = t356*t478
654 64894530 : t480 = t345*t479
655 64894530 : t484 = t342*(t472*t345*t126 + t355*t480/0.2e1_dp)
656 64894530 : t485 = t336*t484
657 : t489 = (t463*t114 - 0.7e1_dp/0.2e1_dp*t109*t465 - (t330*dQndrho &
658 64894530 : *t117*t130) + (2._dp*t330*t485))*t134
659 64894530 : t491 = t376*t435
660 64894530 : t495 = (-t489*t139 - 0.7e1_dp/0.2e1_dp*t374*t491)*E
661 64894530 : t497 = t383*t435
662 64894530 : t500 = dexeindrho(Q, dQndrho, t146)
663 64894530 : t501 = t425*t14
664 64894530 : t504 = t15*t245
665 64894530 : t505 = t504*t149
666 64894530 : t507 = t400*t261
667 64894530 : t508 = t148*t507
668 64894530 : t510 = t404*t435
669 64894530 : t511 = t148*t510
670 64894530 : t513 = 2._dp*t501*t150 + t147*t505 - t147*t508 - t147*t511
671 64894530 : t514 = t513*t409
672 64894530 : t515 = t514*t411
673 : t520 = (-t77*t94*t435 + t447*t95 - 2._dp*t92*t449 + t495*t143 &
674 64894530 : - 3._dp*t142*t497 + t145*(t500 + t515*t418))*Clda
675 64894530 : e_ndrho = e_ndrho + (-t76*t520)*sx0
676 : END IF
677 :
678 87518835 : END SUBROUTINE xwpbe_lda_calc_0
679 :
680 : ! **************************************************************************************************
681 : !> \brief Evaluates the screened hole averaged PBE exchange functional for lda
682 : !> \param e_0 ...
683 : !> \param e_rho ...
684 : !> \param e_ndrho ...
685 : !> \param rho , ndrho: density and norm of the density gradient
686 : !> \param ndrho ...
687 : !> \param sscale scaling factor to enforce Lieb-Oxford bound
688 : !> \param sx0 scaling factor
689 : !> \param order maximum degree of the derivative that should be evaluated
690 : !> \par History
691 : !> 05.2007 created [Manuel Guidon]
692 : !> \author Manuel Guidon
693 : !> \note
694 : !> This routine evaluates the functional for omega=0 using a taylor
695 : !> expansion for the parameter G.
696 : ! **************************************************************************************************
697 307195 : SUBROUTINE xwpbe_lda_calc_01(e_0, e_rho, e_ndrho, rho, ndrho, sscale, sx0, order)
698 : REAL(KIND=dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho
699 : REAL(KIND=dp), INTENT(IN) :: rho, ndrho, sscale, sx0
700 : INTEGER, INTENT(IN) :: order
701 :
702 : REAL(KIND=dp) :: dQndrho, dQrho, Q, t1, t10, t100, t101, t103, t104, t107, t109, t11, t111, &
703 : t112, t113, t114, t115, t116, t118, t12, t122, t124, t125, t126, t129, t13, t130, t131, &
704 : t132, t135, t136, t139, t14, t140, t141, t142, t146, t149, t15, t150, t151, t153, t154, &
705 : t156, t157, t158, t159, t16, t165, t166, t167, t168, t169, t173, t175, t18, t181, t184, &
706 : t185, t186, t188, t189, t19, t190, t191, t193, t194, t197, t199, t2, t20, t202, t203, &
707 : t206, t207, t208, t21, t211, t212, t213, t215, t216, t219, t22, t220, t224, t227, t228, &
708 : t229, t231, t232, t235, t238, t24, t240, t241, t244, t247, t248
709 : REAL(KIND=dp) :: t25, t250, t251, t253, t254, t256, t257, t26, t260, t263, t266, t267, t270, &
710 : t272, t273, t276, t277, t28, t280, t283, t288, t29, t293, t294, t297, t299, t3, t300, &
711 : t301, t304, t305, t307, t308, t31, t311, t314, t315, t317, t318, t319, t32, t321, t322, &
712 : t323, t325, t326, t327, t328, t329, t33, t330, t331, t332, t333, t334, t335, t34, t340, &
713 : t342, t343, t346, t347, t349, t350, t352, t357, t36, t360, t361, t364, t366, t37, t371, &
714 : t372, t375, t376, t379, t380, t383, t385, t388, t389, t39, t392, t393, t395, t396, t398, &
715 : t399, t4, t40, t401, t402, t403, t408, t42, t44, t46, t47, t5
716 : REAL(KIND=dp) :: t52, t53, t54, t56, t58, t59, t6, t61, t63, t65, t66, t67, t68, t69, t7, &
717 : t70, t71, t72, t73, t74, t76, t77, t78, t79, t8, t80, t81, t83, t84, t86, t88, t89, t9, &
718 : t92, t93, t94, t95, t97, t98
719 :
720 307195 : IF (ndrho == 0.0_dp) THEN
721 : ! The cancelling exp(Q)*E1(Q)+log(Q) terms tend to -Euler's constant.
722 : t1 = -f89*clda*(B/(2.0_dp*DD) + C/(2.0_dp*DD**2) + E/DD**3 + &
723 102 : 0.5_dp*A*(-0.577215664901532860606512_dp + LOG(A/(f94*DD))))
724 102 : e_0 = e_0 + sx0*t1*rho**(4.0_dp/3.0_dp)
725 102 : IF (order >= 1) e_rho = e_rho + sx0*(4.0_dp/3.0_dp)*t1*rho**(1.0_dp/3.0_dp)
726 102 : RETURN
727 : END IF
728 :
729 307093 : t1 = ndrho**2
730 307093 : t2 = a1*t1
731 307093 : t3 = r2**2
732 307093 : t4 = 0.1e1_dp/t3
733 307093 : t5 = t2*t4
734 307093 : t6 = pi**2
735 307093 : t7 = r3*t6
736 307093 : t8 = t7*rho
737 307093 : t9 = t8**(0.1e1_dp/0.3e1_dp)
738 307093 : t10 = t9**2
739 307093 : t11 = 0.1e1_dp/t10
740 307093 : t12 = rho**2
741 307093 : t13 = 0.1e1_dp/t12
742 307093 : t14 = t11*t13
743 307093 : t15 = sscale**2
744 307093 : t16 = t14*t15
745 307093 : t18 = t1**2
746 307093 : t19 = a2*t18
747 307093 : t20 = t3**2
748 307093 : t21 = 0.1e1_dp/t20
749 307093 : t22 = t19*t21
750 307093 : t24 = 0.1e1_dp/t9/t8
751 307093 : t25 = t12**2
752 307093 : t26 = 0.1e1_dp/t25
753 307093 : t28 = t15**2
754 307093 : t29 = t24*t26*t28
755 307093 : t31 = t5*t16 + t22*t29
756 307093 : t32 = f94*t31
757 307093 : t33 = a3*t18
758 307093 : t34 = t33*t21
759 307093 : t36 = t18*ndrho
760 307093 : t37 = a4*t36
761 307093 : t39 = 0.1e1_dp/t20/r2
762 307093 : t40 = t37*t39
763 307093 : t42 = 0.1e1_dp/t10/t8
764 307093 : t44 = 0.1e1_dp/t25/rho
765 307093 : t46 = t28*sscale
766 307093 : t47 = t42*t44*t46
767 307093 : t52 = 0.1e1_dp/t20/t3
768 307093 : t53 = a5*t18*t1*t52
769 307093 : t54 = r3**2
770 307093 : t56 = t6**2
771 307093 : t58 = 0.1e1_dp/t54/t56
772 307093 : t59 = t25**2
773 307093 : t61 = t28*t15
774 307093 : t63 = t58/t59*t61
775 307093 : t65 = r1 + t34*t29 + t40*t47 + t53*t63
776 307093 : t66 = 0.1e1_dp/t65
777 307093 : t67 = t66*t1
778 307093 : t68 = t32*t67
779 307093 : t69 = t4*t11
780 307093 : t70 = t13*t15
781 307093 : t71 = 0.1e1_dp/A
782 307093 : t72 = t70*t71
783 307093 : t73 = t69*t72
784 307093 : Q = t68*t73
785 307093 : t74 = rho**(0.1e1_dp/0.3e1_dp)
786 307093 : t76 = t74*rho*f89
787 307093 : t77 = B*f12
788 307093 : t78 = t1*t4
789 307093 : t79 = t78*t11
790 307093 : t80 = t31*t66
791 307093 : t81 = t70*t80
792 307093 : t83 = t79*t81 + DD
793 307093 : t84 = 0.1e1_dp/t83
794 307093 : t86 = F2*t31
795 307093 : t88 = F1 + t86*t66
796 307093 : t89 = t70*t88
797 307093 : t92 = f12*(t79*t89 + r1)
798 307093 : t93 = t83**2
799 307093 : t94 = 0.1e1_dp/t93
800 307093 : t95 = C*t94
801 307093 : t97 = g2*t1
802 307093 : t98 = t97*t4
803 307093 : t100 = g3*t18
804 307093 : t101 = t100*t21
805 307093 : t103 = g1 + t98*t16 + t101*t29
806 307093 : t104 = t70*t103
807 307093 : t107 = (t79*t104 + r1)*E
808 307093 : t109 = 0.1e1_dp/t93/t83
809 307093 : t111 = f12*A
810 307093 : t112 = exei(Q)
811 307093 : t113 = t78*t14
812 307093 : t114 = t15*t31
813 307093 : t115 = t66*t84
814 307093 : t116 = t114*t115
815 307093 : t118 = LOG(t113*t116)
816 : t122 = (t77*t84 + t92*t95 + t107*t109 + t111*(t112 + t118)) &
817 307093 : *Clda
818 307093 : e_0 = e_0 + (-t76*t122)*sx0
819 307093 : IF (order >= 1) THEN
820 300256 : t124 = t4*t42
821 300256 : t125 = t2*t124
822 300256 : t126 = t70*t7
823 300256 : t129 = t12*rho
824 300256 : t130 = 0.1e1_dp/t129
825 300256 : t131 = t11*t130
826 300256 : t132 = t131*t15
827 300256 : t135 = t54*t56
828 300256 : t136 = t135*t12
829 300256 : t139 = t21/t9/t136
830 300256 : t140 = t19*t139
831 300256 : t141 = t26*t28
832 300256 : t142 = t141*t7
833 300256 : t146 = t24*t44*t28
834 : t149 = -0.2e1_dp/0.3e1_dp*t125*t126 - (2._dp*t5*t132) - 0.4e1_dp/ &
835 300256 : 0.3e1_dp*t140*t142 - (4._dp*t22*t146)
836 300256 : t150 = f94*t149
837 300256 : t151 = t150*t67
838 300256 : t153 = t65**2
839 300256 : t154 = 0.1e1_dp/t153
840 300256 : t156 = t154*t1*t4
841 300256 : t157 = t32*t156
842 300256 : t158 = t15*t71
843 300256 : t159 = t33*t139
844 300256 : t165 = 0.1e1_dp/t10/t136
845 300256 : t166 = t39*t165
846 300256 : t167 = t37*t166
847 300256 : t168 = t44*t46
848 300256 : t169 = t168*t7
849 300256 : t173 = 0.1e1_dp/t25/t12
850 300256 : t175 = t42*t173*t46
851 300256 : t181 = t58/t59/rho*t61
852 : t184 = -0.4e1_dp/0.3e1_dp*t159*t142 - (4._dp*t34*t146) - 0.5e1_dp &
853 300256 : /0.3e1_dp*t167*t169 - (5._dp*t40*t175) - (8._dp*t53*t181)
854 300256 : t185 = t158*t184
855 300256 : t186 = t14*t185
856 300256 : t188 = t67*t4
857 300256 : t189 = t32*t188
858 300256 : t190 = t42*t13
859 300256 : t191 = t190*t15
860 300256 : t193 = t71*r3*t6
861 300256 : t194 = t191*t193
862 300256 : t197 = t130*t15
863 300256 : t199 = t69*t197*t71
864 : dQrho = t151*t73 - t157*t186 - 0.2e1_dp/0.3e1_dp*t189*t194 - (2._dp &
865 300256 : *t68*t199)
866 300256 : t202 = a1*ndrho
867 300256 : t203 = t202*t4
868 300256 : t206 = t1*ndrho
869 300256 : t207 = a2*t206
870 300256 : t208 = t207*t21
871 300256 : t211 = 2._dp*t203*t16 + 4._dp*t208*t29
872 300256 : t212 = f94*t211
873 300256 : t213 = t212*t67
874 300256 : t215 = a3*t206
875 300256 : t216 = t215*t21
876 300256 : t219 = a4*t18
877 300256 : t220 = t219*t39
878 300256 : t224 = a5*t36*t52
879 300256 : t227 = 4._dp*t216*t29 + 5._dp*t220*t47 + 6._dp*t224*t63
880 300256 : t228 = t158*t227
881 300256 : t229 = t14*t228
882 300256 : t231 = t66*ndrho
883 300256 : t232 = t32*t231
884 300256 : dQndrho = t213*t73 - t157*t229 + 2._dp*t232*t73
885 300256 : t235 = t74*f89
886 300256 : t238 = t78*t190
887 300256 : t240 = t66*r3*t6
888 300256 : t241 = t114*t240
889 300256 : t244 = t197*t80
890 300256 : t247 = t149*t66
891 300256 : t248 = t70*t247
892 300256 : t250 = t154*t184
893 300256 : t251 = t114*t250
894 : t253 = -0.2e1_dp/0.3e1_dp*t238*t241 - (2._dp*t79*t244) + (t79 &
895 300256 : *t248) - t113*t251
896 300256 : t254 = t94*t253
897 300256 : t256 = t15*t88
898 300256 : t257 = t256*t7
899 300256 : t260 = t197*t88
900 300256 : t263 = F2*t149
901 300256 : t266 = t263*t66 - t86*t250
902 300256 : t267 = t70*t266
903 : t270 = f12*(-0.2e1_dp/0.3e1_dp*t238*t257 - (2._dp*t79*t260) + &
904 300256 : (t79*t267))
905 300256 : t272 = C*t109
906 300256 : t273 = t272*t253
907 300256 : t276 = t15*t103
908 300256 : t277 = t276*t7
909 300256 : t280 = t197*t103
910 300256 : t283 = t97*t124
911 300256 : t288 = t100*t139
912 : t293 = -0.2e1_dp/0.3e1_dp*t283*t126 - (2._dp*t98*t132) - 0.4e1_dp &
913 300256 : /0.3e1_dp*t288*t142 - (4._dp*t101*t146)
914 300256 : t294 = t70*t293
915 : t297 = (-0.2e1_dp/0.3e1_dp*t238*t277 - (2._dp*t79*t280) + (t79 &
916 300256 : *t294))*E
917 300256 : t299 = t93**2
918 300256 : t300 = 0.1e1_dp/t299
919 300256 : t301 = t300*t253
920 300256 : t304 = dexeirho(Q, dQrho, t112)
921 300256 : t305 = t78*t191
922 300256 : t307 = t84*r3*t6
923 300256 : t308 = t80*t307
924 300256 : t311 = t78*t131
925 300256 : t314 = t15*t149
926 300256 : t315 = t314*t115
927 300256 : t317 = t154*t84
928 300256 : t318 = t317*t184
929 300256 : t319 = t114*t318
930 300256 : t321 = t66*t94
931 300256 : t322 = t321*t253
932 300256 : t323 = t114*t322
933 : t325 = -0.2e1_dp/0.3e1_dp*t305*t308 - (2._dp*t311*t116) + t113 &
934 300256 : *t315 - t113*t319 - t113*t323
935 300256 : t326 = 0.1e1_dp/t1
936 300256 : t327 = t325*t326
937 300256 : t328 = t3*t10
938 300256 : t329 = t327*t328
939 300256 : t330 = 0.1e1_dp/t15
940 300256 : t331 = t12*t330
941 300256 : t332 = 0.1e1_dp/t31
942 300256 : t333 = t332*t65
943 300256 : t334 = t333*t83
944 300256 : t335 = t331*t334
945 : t340 = (-t77*t254 + t270*t95 - 2._dp*t92*t273 + t297*t109 - 3._dp &
946 300256 : *t107*t301 + t111*(t304 + t329*t335))*Clda
947 300256 : e_rho = e_rho + (-0.4e1_dp/0.3e1_dp*t235*t122 - t76*t340)*sx0
948 300256 : t342 = ndrho*t4
949 300256 : t343 = t342*t11
950 300256 : t346 = t211*t66
951 300256 : t347 = t70*t346
952 300256 : t349 = t154*t227
953 300256 : t350 = t114*t349
954 300256 : t352 = 2._dp*t343*t81 + t79*t347 - t113*t350
955 300256 : t357 = F2*t211
956 300256 : t360 = t357*t66 - t86*t349
957 300256 : t361 = t70*t360
958 300256 : t364 = f12*(2._dp*t343*t89 + t79*t361)
959 300256 : t366 = t272*t352
960 300256 : t371 = g2*ndrho
961 300256 : t372 = t371*t4
962 300256 : t375 = g3*t206
963 300256 : t376 = t375*t21
964 300256 : t379 = 2._dp*t372*t16 + 4._dp*t376*t29
965 300256 : t380 = t70*t379
966 300256 : t383 = (2._dp*t343*t104 + t79*t380)*E
967 300256 : t385 = t300*t352
968 300256 : t388 = dexeindrho(Q, dQndrho, t112)
969 300256 : t389 = t342*t14
970 300256 : t392 = t15*t211
971 300256 : t393 = t392*t115
972 300256 : t395 = t317*t227
973 300256 : t396 = t114*t395
974 300256 : t398 = t321*t352
975 300256 : t399 = t114*t398
976 300256 : t401 = 2._dp*t389*t116 + t113*t393 - t113*t396 - t113*t399
977 300256 : t402 = t401*t326
978 300256 : t403 = t402*t328
979 : t408 = (-t77*t94*t352 + t364*t95 - 2._dp*t92*t366 + t383*t109 &
980 300256 : - 3._dp*t107*t385 + t111*(t388 + t403*t335))*Clda
981 300256 : e_ndrho = e_ndrho + (-t76*t408)*sx0
982 : END IF
983 :
984 : END SUBROUTINE xwpbe_lda_calc_01
985 :
986 : ! **************************************************************************************************
987 : !> \brief Evaluates the screened hole averaged PBE exchange functional for lda.
988 : !> \param e_0 ...
989 : !> \param e_rho ...
990 : !> \param e_ndrho ...
991 : !> \param rho , ndrho: density and norm of the density gradient
992 : !> \param ndrho ...
993 : !> \param omega scaling factor
994 : !> \param sscale scaling factor to enforce Lieb-Oxford bound
995 : !> \param sx scaling factor
996 : !> \param order maximum degree of the derivative that should be evaluated
997 : !> \par History
998 : !> 05.2007 created [Manuel Guidon]
999 : !> \author Manuel Guidon
1000 : !> \note
1001 : !> This routine evaluates the exact functional for omega!=0.
1002 : ! **************************************************************************************************
1003 24453786 : SUBROUTINE xwpbe_lda_calc_1(e_0, e_rho, e_ndrho, rho, ndrho, omega, sscale, sx, order)
1004 : REAL(KIND=dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho
1005 : REAL(KIND=dp), INTENT(IN) :: rho, ndrho, omega, sscale, sx
1006 : INTEGER, INTENT(IN) :: order
1007 :
1008 : REAL(KIND=dp) :: dQndrho, dQrho, Q, t1, t10, t1001, t1009, t1011, t102, t1024, t1035, t104, &
1009 : t105, t1052, t1066, t1068, t1069, t108, t109, t11, t110, t111, t113, t115, t116, t117, &
1010 : t118, t119, t12, t120, t121, t123, t124, t125, t126, t127, t128, t129, t13, t131, t132, &
1011 : t133, t136, t14, t140, t141, t142, t144, t145, t147, t148, t149, t15, t150, t151, t152, &
1012 : t153, t154, t155, t156, t157, t158, t159, t16, t162, t163, t166, t167, t168, t169, t17, &
1013 : t170, t174, t176, t177, t18, t180, t181, t185, t186, t189, t19, t190, t192, t193, t194, &
1014 : t195, t197, t198, t199, t2, t200, t202, t205, t209, t21, t210
1015 : REAL(KIND=dp) :: t213, t214, t219, t22, t221, t222, t223, t224, t227, t228, t23, t233, t237, &
1016 : t239, t24, t241, t243, t246, t247, t25, t250, t251, t252, t253, t255, t256, t257, t258, &
1017 : t263, t264, t265, t267, t27, t270, t271, t273, t275, t276, t28, t280, t284, t285, t286, &
1018 : t288, t289, t29, t290, t293, t294, t295, t298, t3, t303, t305, t307, t31, t310, t311, &
1019 : t312, t313, t315, t316, t318, t319, t32, t321, t326, t329, t330, t332, t333, t335, t336, &
1020 : t338, t339, t34, t340, t342, t343, t345, t346, t347, t348, t349, t35, t351, t352, t353, &
1021 : t354, t357, t358, t36, t361, t362, t363, t364, t368, t371, t372
1022 : REAL(KIND=dp) :: t373, t374, t375, t376, t377, t378, t38, t383, t384, t385, t386, t39, t390, &
1023 : t392, t398, t4, t401, t402, t403, t404, t406, t409, t41, t411, t412, t415, t416, t419, &
1024 : t42, t420, t421, t424, t425, t426, t428, t429, t432, t433, t437, t44, t440, t441, t442, &
1025 : t444, t446, t449, t452, t453, t454, t457, t46, t461, t463, t465, t468, t469, t472, t475, &
1026 : t478, t479, t48, t481, t486, t49, t493, t495, t498, t499, t5, t502, t503, t504, t507, &
1027 : t508, t509, t510, t511, t513, t514, t515, t517, t518, t522, t525, t529, t530, t531, t532, &
1028 : t533, t534, t535, t537, t538, t539, t54, t540, t542, t55, t550
1029 : REAL(KIND=dp) :: t551, t552, t553, t554, t556, t557, t558, t56, t562, t563, t567, t569, &
1030 : t571, t572, t575, t576, t577, t578, t58, t581, t582, t583, t587, t588, t589, t592, t596, &
1031 : t597, t6, t60, t602, t603, t608, t61, t613, t617, t620, t621, t624, t626, t629, t63, &
1032 : t631, t634, t638, t639, t644, t645, t65, t653, t655, t656, t657, t658, t659, t663, t666, &
1033 : t669, t67, t673, t679, t68, t681, t685, t689, t69, t690, t691, t697, t698, t7, t70, t701, &
1034 : t71, t711, t713, t717, t718, t72, t721, t728, t73, t735, t736, t739, t74, t740, t746, &
1035 : t747, t748, t75, t752, t753, t754, t755, t759, t761, t763, t765, t769
1036 : REAL(KIND=dp) :: t77, t771, t772, t779, t78, t780, t781, t783, t784, t787, t791, t792, t799, &
1037 : t8, t80, t800, t801, t803, t804, t808, t81, t810, t812, t815, t816, t817, t82, t820, &
1038 : t823, t826, t83, t832, t834, t836, t84, t842, t845, t846, t848, t849, t85, t851, t865, &
1039 : t867, t87, t870, t871, t873, t874, t876, t877, t88, t880, t884, t885, t888, t889, t891, &
1040 : t892, t893, t897, t898, t9, t90, t902, t904, t907, t908, t909, t91, t910, t916, t918, &
1041 : t919, t92, t93, t930, t932, t933, t937, t94, t942, t945, t95, t951, t954, t958, t96, &
1042 : t962, t965, t968, t97, t971, t972, t979, t982, t988, t99
1043 :
1044 24453786 : t1 = ndrho**2
1045 24453786 : t2 = r2**2
1046 24453786 : t3 = 0.1e1_dp/t2
1047 24453786 : t4 = t1*t3
1048 24453786 : t5 = pi**2
1049 24453786 : t6 = r3*t5
1050 24453786 : t7 = t6*rho
1051 24453786 : t8 = t7**(0.1e1_dp/0.3e1_dp)
1052 24453786 : t9 = t8**2
1053 24453786 : t10 = 0.1e1_dp/t9
1054 24453786 : t11 = t4*t10
1055 24453786 : t12 = rho**2
1056 24453786 : t13 = 0.1e1_dp/t12
1057 24453786 : t14 = sscale**2
1058 24453786 : t15 = t13*t14
1059 24453786 : t16 = a1*t1
1060 24453786 : t17 = t16*t3
1061 24453786 : t18 = t10*t13
1062 24453786 : t19 = t18*t14
1063 24453786 : t21 = t1**2
1064 24453786 : t22 = a2*t21
1065 24453786 : t23 = t2**2
1066 24453786 : t24 = 0.1e1_dp/t23
1067 24453786 : t25 = t22*t24
1068 24453786 : t27 = 0.1e1_dp/t8/t7
1069 24453786 : t28 = t12**2
1070 24453786 : t29 = 0.1e1_dp/t28
1071 24453786 : t31 = t14**2
1072 24453786 : t32 = t27*t29*t31
1073 24453786 : t34 = t17*t19 + t25*t32
1074 24453786 : t35 = a3*t21
1075 24453786 : t36 = t35*t24
1076 24453786 : t38 = t21*ndrho
1077 24453786 : t39 = a4*t38
1078 24453786 : t41 = 0.1e1_dp/t23/r2
1079 24453786 : t42 = t39*t41
1080 24453786 : t44 = 0.1e1_dp/t9/t7
1081 24453786 : t46 = 0.1e1_dp/t28/rho
1082 24453786 : t48 = t31*sscale
1083 24453786 : t49 = t44*t46*t48
1084 24453786 : t54 = 0.1e1_dp/t23/t2
1085 24453786 : t55 = a5*t21*t1*t54
1086 24453786 : t56 = r3**2
1087 24453786 : t58 = t5**2
1088 24453786 : t60 = 0.1e1_dp/t56/t58
1089 24453786 : t61 = t28**2
1090 24453786 : t63 = t31*t14
1091 24453786 : t65 = t60/t61*t63
1092 24453786 : t67 = r1 + t36*t32 + t42*t49 + t55*t65
1093 24453786 : t68 = 0.1e1_dp/t67
1094 24453786 : t69 = t34*t68
1095 24453786 : t70 = t15*t69
1096 24453786 : t71 = t11*t70
1097 24453786 : t72 = omega**2
1098 24453786 : t73 = beta*t72
1099 24453786 : t74 = t73*t10
1100 24453786 : t75 = t71 + t74
1101 24453786 : t77 = 0.1e1_dp/A
1102 24453786 : Q = f94*t75*t77
1103 24453786 : t78 = rho**(0.1e1_dp/0.3e1_dp)
1104 24453786 : t80 = t78*rho*f89
1105 24453786 : t81 = B*f12
1106 24453786 : t82 = t71 + DD
1107 24453786 : t83 = 0.1e1_dp/t82
1108 24453786 : t84 = t81*t83
1109 24453786 : t85 = F2*t34
1110 24453786 : t87 = F1 + t85*t68
1111 24453786 : t88 = t15*t87
1112 24453786 : t90 = t11*t88 + r1
1113 24453786 : t91 = f12*t90
1114 24453786 : t92 = t82**2
1115 24453786 : t93 = 0.1e1_dp/t92
1116 24453786 : t94 = C*t93
1117 24453786 : t95 = t91*t94
1118 24453786 : t96 = f34*pi
1119 24453786 : t97 = rootpi
1120 24453786 : t99 = r6*C
1121 24453786 : t102 = r4*B
1122 24453786 : t104 = r8*A
1123 24453786 : t105 = t92*t82
1124 24453786 : t108 = t97*(r15*E + t99*t90*t82 + t102*t92 + t104*t105)
1125 24453786 : t109 = 0.1e1_dp/r16
1126 24453786 : t110 = SQRT(t82)
1127 24453786 : t111 = t110*t105
1128 24453786 : t113 = t109/t111
1129 24453786 : t115 = SQRT(A)
1130 24453786 : t116 = f94*t34
1131 24453786 : t117 = t68*t1
1132 24453786 : t118 = t116*t117
1133 24453786 : t119 = t3*t10
1134 24453786 : t120 = t15*t77
1135 24453786 : t121 = t119*t120
1136 24453786 : t123 = EXP(t118*t121)
1137 24453786 : t124 = t115*t123
1138 24453786 : t125 = f32*ndrho
1139 24453786 : t126 = 0.1e1_dp/r2
1140 24453786 : t127 = t125*t126
1141 24453786 : t128 = 0.1e1_dp/t8
1142 24453786 : t129 = 0.1e1_dp/rho
1143 24453786 : t131 = t69*t77
1144 24453786 : t132 = SQRT(t131)
1145 24453786 : t133 = sscale*t132
1146 24453786 : t136 = erfc(t127*t128*t129*t133)
1147 24453786 : t140 = 0.1e1_dp/f1516
1148 24453786 : t141 = (t96 + t108*t113 - t96*t124*t136)*t140
1149 24453786 : t142 = 0.1e1_dp/t97
1150 24453786 : t144 = 0.1e1_dp/E
1151 24453786 : t145 = t142*t111*t144
1152 24453786 : t147 = -t141*t145 + r1
1153 24453786 : t148 = t147*E
1154 24453786 : t149 = 0.1e1_dp/t105
1155 24453786 : t150 = t148*t149
1156 24453786 : t151 = f158*E
1157 24453786 : t152 = t147*t83
1158 24453786 : t153 = t72*t10
1159 24453786 : t154 = t71 + DD + t153
1160 24453786 : t155 = t154**2
1161 24453786 : t156 = t155**2
1162 24453786 : t157 = t156*t154
1163 24453786 : t158 = SQRT(t157)
1164 24453786 : t159 = 0.1e1_dp/t158
1165 24453786 : t162 = SQRT(t154)
1166 24453786 : t163 = 0.1e1_dp/t162
1167 24453786 : t166 = f68*C
1168 24453786 : t167 = t90*t83
1169 24453786 : t168 = t155*t154
1170 24453786 : t169 = SQRT(t168)
1171 24453786 : t170 = 0.1e1_dp/t169
1172 : t174 = (-t151*t152*t159 - t81*t83*t163 - t166*t167*t170) &
1173 24453786 : *omega
1174 24453786 : t176 = f52*E
1175 24453786 : t177 = t147*t93
1176 24453786 : t180 = f12*C
1177 24453786 : t181 = t90*t93
1178 24453786 : t185 = t72*omega
1179 24453786 : t186 = (-t176*t177*t159 - t180*t181*t170)*t185
1180 24453786 : t189 = 0.1e1_dp/r3/t5
1181 24453786 : t190 = t189*t129
1182 24453786 : t192 = t72**2
1183 24453786 : t193 = t192*omega
1184 24453786 : t194 = t159*t193
1185 24453786 : t195 = t194*t44
1186 24453786 : t197 = f12*A
1187 24453786 : t198 = exei(Q)
1188 24453786 : t199 = t71 + DD + t74
1189 24453786 : t200 = 0.1e1_dp/t199
1190 24453786 : t202 = LOG(t75*t200)
1191 24453786 : t205 = SQRT(t199)
1192 24453786 : t209 = t115*f34
1193 24453786 : t210 = exer(Q)
1194 24453786 : t213 = (t197*t97/t205 - t209*t210)*alpha1
1195 24453786 : t214 = omega*t128
1196 24453786 : t219 = (t197*t200 - f98*t198)*alpha2
1197 24453786 : t221 = A*f14
1198 24453786 : t222 = t199**2
1199 24453786 : t223 = t222*t199
1200 24453786 : t224 = SQRT(t223)
1201 24453786 : t227 = SQRT(t75)
1202 24453786 : t228 = 0.1e1_dp/t227
1203 24453786 : t233 = 0.1e1_dp/t115
1204 : t237 = (t97*(t221/t224 - f98*t228) + f2716*t210*t233)*alpha3 &
1205 24453786 : *t185
1206 24453786 : t239 = 0.1e1_dp/t75
1207 24453786 : t241 = 0.1e1_dp/t222
1208 24453786 : t243 = f8132*t77
1209 24453786 : t246 = (-f98*t239 + t197*t241 + t243*t198)*alpha4
1210 24453786 : t247 = t192*t27
1211 24453786 : t250 = t75**2
1212 24453786 : t251 = t250*t75
1213 24453786 : t252 = SQRT(t251)
1214 24453786 : t253 = 0.1e1_dp/t252
1215 24453786 : t255 = f38*A
1216 24453786 : t256 = t222**2
1217 24453786 : t257 = t256*t199
1218 24453786 : t258 = SQRT(t257)
1219 24453786 : t263 = A**2
1220 24453786 : t264 = t263*A
1221 24453786 : t265 = SQRT(t264)
1222 24453786 : t267 = f24364/t265
1223 : t270 = (t97*(t243*t228 - f916*t253 + t255/t258) - t267*t210) &
1224 24453786 : *alpha5
1225 24453786 : t271 = t193*t44
1226 24453786 : t273 = 0.1e1_dp/t223
1227 24453786 : t275 = 0.1e1_dp/t250
1228 24453786 : t276 = f98*t275
1229 24453786 : t280 = f729128/t263
1230 24453786 : t284 = t192*t72
1231 : t285 = (A*t273 - t276 + t243*r1*t239 - t280*t198)*alpha6 &
1232 24453786 : *t284
1233 24453786 : t286 = t60*t13
1234 24453786 : t288 = f1516*A
1235 24453786 : t289 = t256*t223
1236 24453786 : t290 = SQRT(t289)
1237 24453786 : t293 = t250**2
1238 24453786 : t294 = t293*t75
1239 24453786 : t295 = SQRT(t294)
1240 24453786 : t298 = f8164*t77
1241 24453786 : t303 = t263**2
1242 24453786 : t305 = SQRT(t303*A)
1243 24453786 : t307 = f2187256/t305
1244 : t310 = (t97*(t288/t290 - f2732/t295 + t298*t253 - t280*t228) &
1245 24453786 : + t307*t210)*alpha7
1246 24453786 : t311 = t192*t185
1247 24453786 : t312 = t56*t58
1248 24453786 : t313 = t312*t12
1249 24453786 : t315 = 0.1e1_dp/t8/t313
1250 24453786 : t316 = t311*t315
1251 24453786 : t318 = r3*A
1252 24453786 : t319 = 0.1e1_dp/t256
1253 24453786 : t321 = 0.1e1_dp/t251
1254 24453786 : t326 = f6561512/t264
1255 : t329 = (t318*t319 - f94*t321 + t243*t275 - t280*t239 + t326 &
1256 24453786 : *t198)*alpha8
1257 24453786 : t330 = t192**2
1258 24453786 : t332 = 0.1e1_dp/t9/t313
1259 24453786 : t333 = t330*t332
1260 : t335 = t84 + t95 + t150 + t174*t128 + t186*t190 - t150*t195 + &
1261 : t197*(t198 + t202) + t213*t214 + t219*t153 + t237*t190 + &
1262 24453786 : t246*t247 + t270*t271 + t285*t286 + t310*t316 + t329*t333
1263 24453786 : t336 = t335*Clda
1264 24453786 : e_0 = e_0 + (-t80*t336)*sx
1265 24453786 : IF (order >= 1) THEN
1266 15526556 : t338 = t44*t13
1267 15526556 : t339 = t4*t338
1268 15526556 : t340 = t14*t34
1269 15526556 : t342 = t68*r3*t5
1270 15526556 : t343 = t340*t342
1271 15526556 : t345 = 0.2e1_dp/0.3e1_dp*t339*t343
1272 15526556 : t346 = t12*rho
1273 15526556 : t347 = 0.1e1_dp/t346
1274 15526556 : t348 = t347*t14
1275 15526556 : t349 = t348*t69
1276 15526556 : t351 = 2._dp*t11*t349
1277 15526556 : t352 = t3*t44
1278 15526556 : t353 = t16*t352
1279 15526556 : t354 = t15*t6
1280 15526556 : t357 = t10*t347
1281 15526556 : t358 = t357*t14
1282 15526556 : t361 = t24*t315
1283 15526556 : t362 = t22*t361
1284 15526556 : t363 = t29*t31
1285 15526556 : t364 = t363*t6
1286 15526556 : t368 = t27*t46*t31
1287 : t371 = -0.2e1_dp/0.3e1_dp*t353*t354 - (2._dp*t17*t358) - 0.4e1_dp &
1288 15526556 : /0.3e1_dp*t362*t364 - (4._dp*t25*t368)
1289 15526556 : t372 = t371*t68
1290 15526556 : t373 = t15*t372
1291 15526556 : t374 = t11*t373
1292 15526556 : t375 = t4*t18
1293 15526556 : t376 = t67**2
1294 15526556 : t377 = 0.1e1_dp/t376
1295 15526556 : t378 = t35*t361
1296 15526556 : t383 = t41*t332
1297 15526556 : t384 = t39*t383
1298 15526556 : t385 = t46*t48
1299 15526556 : t386 = t385*t6
1300 15526556 : t390 = 0.1e1_dp/t28/t12
1301 15526556 : t392 = t44*t390*t48
1302 15526556 : t398 = t60/t61/rho*t63
1303 : t401 = -0.4e1_dp/0.3e1_dp*t378*t364 - (4._dp*t36*t368) - 0.5e1_dp &
1304 15526556 : /0.3e1_dp*t384*t386 - (5._dp*t42*t392) - (8._dp*t55*t398)
1305 15526556 : t402 = t377*t401
1306 15526556 : t403 = t340*t402
1307 15526556 : t404 = t375*t403
1308 15526556 : t406 = t44*r3*t5
1309 15526556 : t409 = -t345 - t351 + t374 - t404 - 0.2e1_dp/0.3e1_dp*t73*t406
1310 15526556 : dQrho = f94*t409*t77
1311 15526556 : t411 = ndrho*t3
1312 15526556 : t412 = t411*t10
1313 15526556 : t415 = a1*ndrho
1314 15526556 : t416 = t415*t3
1315 15526556 : t419 = t1*ndrho
1316 15526556 : t420 = a2*t419
1317 15526556 : t421 = t420*t24
1318 15526556 : t424 = 2._dp*t416*t19 + 4._dp*t421*t32
1319 15526556 : t425 = t424*t68
1320 15526556 : t426 = t15*t425
1321 15526556 : t428 = a3*t419
1322 15526556 : t429 = t428*t24
1323 15526556 : t432 = a4*t21
1324 15526556 : t433 = t432*t41
1325 15526556 : t437 = a5*t38*t54
1326 15526556 : t440 = 4._dp*t429*t32 + 5._dp*t433*t49 + 6._dp*t437*t65
1327 15526556 : t441 = t377*t440
1328 15526556 : t442 = t340*t441
1329 15526556 : t444 = 2._dp*t412*t70 + t11*t426 - t375*t442
1330 15526556 : dQndrho = f94*t444*t77
1331 15526556 : t446 = t78*f89
1332 15526556 : t449 = t60*t347
1333 15526556 : t452 = C*t149
1334 15526556 : t453 = -t345 - t351 + t374 - t404
1335 15526556 : t454 = t452*t453
1336 15526556 : t457 = t329*t330
1337 15526556 : t461 = t56*r3*t58*t5*t346
1338 15526556 : t463 = 0.1e1_dp/t9/t461
1339 15526556 : t465 = t463*r3*t5
1340 15526556 : t468 = t14*t87
1341 15526556 : t469 = t468*t6
1342 15526556 : t472 = t348*t87
1343 15526556 : t475 = F2*t371
1344 15526556 : t478 = t475*t68 - t85*t402
1345 15526556 : t479 = t15*t478
1346 : t481 = -0.2e1_dp/0.3e1_dp*t339*t469 - (2._dp*t11*t472) + (t11 &
1347 15526556 : *t479)
1348 15526556 : t486 = t82*t453
1349 : t493 = t97*(t99*t481*t82 + t99*t90*t453 + 2._dp*t102*t486 &
1350 15526556 : + 3._dp*t104*t92*t453)
1351 15526556 : t495 = t92**2
1352 15526556 : t498 = t109/t110/t495
1353 15526556 : t499 = t498*t453
1354 15526556 : t502 = t96*t115
1355 15526556 : t503 = f94*t371
1356 15526556 : t504 = t503*t117
1357 15526556 : t507 = t377*t1*t3
1358 15526556 : t508 = t116*t507
1359 15526556 : t509 = t14*t77
1360 15526556 : t510 = t509*t401
1361 15526556 : t511 = t18*t510
1362 15526556 : t513 = t117*t3
1363 15526556 : t514 = t116*t513
1364 15526556 : t515 = t338*t14
1365 15526556 : t517 = t77*r3*t5
1366 15526556 : t518 = t515*t517
1367 15526556 : t522 = t119*t348*t77
1368 : t525 = t504*t121 - t508*t511 - 0.2e1_dp/0.3e1_dp*t514*t518 - (2._dp &
1369 15526556 : *t118*t522)
1370 15526556 : t529 = rootpi
1371 15526556 : t530 = 0.1e1_dp/t529
1372 15526556 : t531 = t123*t530
1373 15526556 : t532 = f32**2
1374 15526556 : t533 = t532*t1
1375 15526556 : t534 = t533*t119
1376 15526556 : t535 = t15*t131
1377 15526556 : t537 = EXP(-t534*t535)
1378 15526556 : t538 = t126*t27
1379 15526556 : t539 = t125*t538
1380 15526556 : t540 = t129*sscale
1381 15526556 : t542 = t132*r3*t5
1382 15526556 : t550 = t125*t126*t128
1383 15526556 : t551 = 0.1e1_dp/t132
1384 15526556 : t552 = t372*t77
1385 15526556 : t553 = t34*t377
1386 15526556 : t554 = t77*t401
1387 15526556 : t556 = t552 - t553*t554
1388 15526556 : t557 = t551*t556
1389 15526556 : t558 = t540*t557
1390 : t562 = t537*(-t539*t540*t542/0.3e1_dp - t127*t128*t13*t133 &
1391 15526556 : + t550*t558/0.2e1_dp)
1392 15526556 : t563 = t531*t562
1393 : t567 = (t493*t113 - 0.7e1_dp/0.2e1_dp*t108*t499 - (t502*t525 &
1394 15526556 : *t123*t136) + (2._dp*t502*t563))*t140
1395 15526556 : t569 = t141*t142
1396 15526556 : t571 = t110*t92*t144
1397 15526556 : t572 = t571*t453
1398 15526556 : t575 = -t567*t145 - 0.7e1_dp/0.2e1_dp*t569*t572
1399 15526556 : t576 = t575*E
1400 15526556 : t577 = t576*t149
1401 15526556 : t578 = t189*t13
1402 15526556 : t581 = 0.1e1_dp/t158/t157
1403 15526556 : t582 = t149*t581
1404 15526556 : t583 = t148*t582
1405 15526556 : t587 = -t345 - t351 + t374 - t404 - 0.2e1_dp/0.3e1_dp*t72*t44*t6
1406 15526556 : t588 = t156*t587
1407 15526556 : t589 = t271*t588
1408 15526556 : t592 = t219*t72
1409 15526556 : t596 = 0.1e1_dp/t224/t223
1410 15526556 : t597 = t596*t222
1411 15526556 : t602 = 0.1e1_dp/t227/t75
1412 15526556 : t603 = f98*t602
1413 15526556 : t608 = dexerrho(Q, dQrho)
1414 : t613 = (t97*(-0.3e1_dp/0.2e1_dp*t221*t597*t409 + t603*t409/ &
1415 15526556 : 0.2e1_dp) + f2716*t608*t233)*alpha3*t185
1416 15526556 : t617 = f12*t481
1417 15526556 : t620 = 0.1e1_dp/t495
1418 15526556 : t621 = t620*t453
1419 15526556 : t624 = t213*omega
1420 15526556 : t626 = t27*r3*t5
1421 15526556 : t629 = t246*t192
1422 15526556 : t631 = t315*r3*t5
1423 15526556 : t634 = t602*t409
1424 15526556 : t638 = 0.1e1_dp/t252/t251
1425 15526556 : t639 = f916*t638
1426 15526556 : t644 = 0.1e1_dp/t258/t257
1427 15526556 : t645 = t644*t256
1428 : t653 = (t97*(-t243*t634/0.2e1_dp + 0.3e1_dp/0.2e1_dp*t639*t250* &
1429 15526556 : t409 - 0.5e1_dp/0.2e1_dp*t255*t645*t409) - t267*t608)*alpha5
1430 : t655 = -(2._dp*t285*t449) - (2._dp*t91*t454) - 0.8e1_dp/0.3e1_dp &
1431 : *t457*t465 + t577 - t186*t578 + 0.5e1_dp/0.2e1_dp*t583*t589 &
1432 : - 0.2e1_dp/0.3e1_dp*t592*t406 + t613*t190 - t81*t93*t453 + &
1433 : t617*t94 - t237*t578 - (3._dp*t148*t621) - t624*t626/0.3e1_dp &
1434 15526556 : - 0.4e1_dp/0.3e1_dp*t629*t631 + t653*t271
1435 15526556 : t656 = t149*t159
1436 15526556 : t657 = t148*t656
1437 15526556 : t658 = t193*t332
1438 15526556 : t659 = t658*t6
1439 15526556 : t663 = t273*t409
1440 15526556 : t666 = dexeirho(Q, dQrho)
1441 15526556 : t669 = (t276*t409 - 2._dp*t197*t663 + t243*t666)*alpha4
1442 15526556 : t673 = t97/t205/t199
1443 15526556 : t679 = (-t197*t673*t409/0.2e1_dp - t209*t608)*alpha1
1444 15526556 : t681 = t241*t409
1445 15526556 : t685 = (-t197*t681 - f98*t666)*alpha2
1446 15526556 : t689 = 0.1e1_dp/t290/t289
1447 15526556 : t690 = t256*t222
1448 15526556 : t691 = t689*t690
1449 15526556 : t697 = f2732/t295/t294
1450 15526556 : t698 = t293*t409
1451 15526556 : t701 = t638*t250
1452 : t711 = (t97*(-0.7e1_dp/0.2e1_dp*t288*t691*t409 + 0.5e1_dp/0.2e1_dp &
1453 : *t697*t698 - 0.3e1_dp/0.2e1_dp*t298*t701*t409 + t280*t634/ &
1454 15526556 : 0.2e1_dp) + t307*t608)*alpha7
1455 15526556 : t713 = 0.1e1_dp/t257
1456 15526556 : t717 = 0.1e1_dp/t293
1457 15526556 : t718 = f94*t717
1458 15526556 : t721 = t321*t409
1459 : t728 = (-4._dp*t318*t713*t409 + 3._dp*t718*t409 - 2._dp*t243*t721 &
1460 15526556 : + t280*t275*t409 + t326*t666)*alpha8
1461 15526556 : t735 = t176*t147
1462 15526556 : t736 = t656*t453
1463 15526556 : t739 = t93*t581
1464 15526556 : t740 = t739*t588
1465 15526556 : t746 = t180*t90
1466 15526556 : t747 = t149*t170
1467 15526556 : t748 = t747*t453
1468 15526556 : t752 = 0.1e1_dp/t169/t168
1469 15526556 : t753 = t93*t752
1470 15526556 : t754 = t155*t587
1471 15526556 : t755 = t753*t754
1472 : t759 = (-t176*t575*t93*t159 + (2._dp*t735*t736) + 0.5e1_dp/ &
1473 : 0.2e1_dp*(t735)*(t740) - t180*t481*t93*t170 + (2._dp &
1474 15526556 : *t746*t748) + 0.3e1_dp/0.2e1_dp*(t746)*(t755))*t185
1475 15526556 : t761 = t310*t311
1476 15526556 : t763 = 0.1e1_dp/t8/t461
1477 15526556 : t765 = t763*r3*t5
1478 15526556 : t769 = t75*t241
1479 15526556 : t771 = t409*t200 - t769*t409
1480 15526556 : t772 = t771*t239
1481 15526556 : t779 = t151*t147
1482 15526556 : t780 = t93*t159
1483 15526556 : t781 = t780*t453
1484 15526556 : t783 = t83*t581
1485 15526556 : t784 = t783*t588
1486 15526556 : t787 = t93*t163
1487 15526556 : t791 = 0.1e1_dp/t162/t154
1488 15526556 : t792 = t83*t791
1489 15526556 : t799 = t166*t90
1490 15526556 : t800 = t93*t170
1491 15526556 : t801 = t800*t453
1492 15526556 : t803 = t83*t752
1493 15526556 : t804 = t803*t754
1494 : t808 = (-t151*t575*t83*t159 + t779*t781 + 0.5e1_dp/0.2e1_dp*t779 &
1495 : *t784 + t81*t787*t453 + t81*t792*t587/0.2e1_dp - t166 &
1496 : *t481*t83*t170 + t799*t801 + 0.3e1_dp/0.2e1_dp*t799*t804)* &
1497 15526556 : omega
1498 15526556 : t810 = t148*t620
1499 15526556 : t812 = t194*t44*t453
1500 15526556 : t815 = t270*t193
1501 15526556 : t816 = t332*r3
1502 15526556 : t817 = t816*t5
1503 15526556 : t820 = A*t319
1504 15526556 : t823 = f98*t321
1505 15526556 : t826 = r1*t275
1506 : t832 = (-3._dp*t820*t409 + 2._dp*t823*t409 - t243*t826*t409 - t280 &
1507 15526556 : *t666)*alpha6*t284
1508 : t834 = 0.5e1_dp/0.3e1_dp*t657*t659 + t669*t247 + t679*t214 + t685 &
1509 : *t153 - t577*t195 + t711*t316 + t728*t333 - t174*t626 &
1510 : /0.3e1_dp + t759*t190 - 0.7e1_dp/0.3e1_dp*t761*t765 + t197*(t666 &
1511 : + t772*t199) + t808*t128 + (3._dp*t810*t812) - 0.5e1_dp/0.3e1_dp &
1512 15526556 : *t815*t817 + t832*t286
1513 15526556 : t836 = (t655 + t834)*Clda
1514 15526556 : e_rho = e_rho + (-0.4e1_dp/0.3e1_dp*t446*t336 - t80*t836)*sx
1515 15526556 : t842 = F2*t424
1516 15526556 : t845 = t842*t68 - t85*t441
1517 15526556 : t846 = t15*t845
1518 15526556 : t848 = 2._dp*t412*t88 + t11*t846
1519 15526556 : t849 = f12*t848
1520 15526556 : t851 = t452*t444
1521 : t865 = t97*(t99*t848*t82 + t99*t90*t444 + 2._dp*t102*t82 &
1522 15526556 : *t444 + 3._dp*t104*t92*t444)
1523 15526556 : t867 = t498*t444
1524 15526556 : t870 = f94*t424
1525 15526556 : t871 = t870*t117
1526 15526556 : t873 = t509*t440
1527 15526556 : t874 = t18*t873
1528 15526556 : t876 = t68*ndrho
1529 15526556 : t877 = t116*t876
1530 15526556 : t880 = t871*t121 - t508*t874 + 2._dp*t877*t121
1531 15526556 : t884 = f32*t126
1532 15526556 : t885 = t884*t128
1533 15526556 : t888 = t425*t77
1534 15526556 : t889 = t77*t440
1535 15526556 : t891 = t888 - t553*t889
1536 15526556 : t892 = t551*t891
1537 15526556 : t893 = t540*t892
1538 15526556 : t897 = t537*(t885*t540*t132 + t550*t893/0.2e1_dp)
1539 15526556 : t898 = t531*t897
1540 : t902 = (t865*t113 - 0.7e1_dp/0.2e1_dp*t108*t867 - (t502*t880 &
1541 15526556 : *t123*t136) + (2._dp*t502*t898))*t140
1542 15526556 : t904 = t571*t444
1543 15526556 : t907 = -t902*t145 - 0.7e1_dp/0.2e1_dp*t569*t904
1544 15526556 : t908 = t907*E
1545 15526556 : t909 = t908*t149
1546 15526556 : t910 = t620*t444
1547 15526556 : t916 = t780*t444
1548 15526556 : t918 = t156*t444
1549 15526556 : t919 = t783*t918
1550 15526556 : t930 = t800*t444
1551 15526556 : t932 = t155*t444
1552 15526556 : t933 = t803*t932
1553 : t937 = (-t151*t907*t83*t159 + t779*t916 + 0.5e1_dp/0.2e1_dp*t779 &
1554 : *t919 + t81*t787*t444 + t81*t792*t444/0.2e1_dp - t166 &
1555 : *t848*t83*t170 + t799*t930 + 0.3e1_dp/0.2e1_dp*t799*t933)* &
1556 15526556 : omega
1557 15526556 : t942 = t656*t444
1558 15526556 : t945 = t739*t918
1559 15526556 : t951 = t747*t444
1560 15526556 : t954 = t753*t932
1561 : t958 = (-t176*t907*t93*t159 + (2._dp*t735*t942) + 0.5e1_dp/ &
1562 : 0.2e1_dp*(t735)*(t945) - t180*t848*t93*t170 + (2._dp &
1563 15526556 : *t746*t951) + 0.3e1_dp/0.2e1_dp*(t746)*(t954))*t185
1564 15526556 : t962 = t194*t44*t444
1565 15526556 : t965 = t271*t918
1566 15526556 : t968 = dexeindrho(Q, dQndrho)
1567 15526556 : t971 = t444*t200 - t769*t444
1568 15526556 : t972 = t971*t239
1569 15526556 : t979 = dexerndrho(Q, dQndrho)
1570 15526556 : t982 = (-t197*t673*t444/0.2e1_dp - t209*t979)*alpha1
1571 15526556 : t988 = (-t197*t241*t444 - f98*t968)*alpha2
1572 : t1001 = (t97*(-0.3e1_dp/0.2e1_dp*t221*t597*t444 + t603*t444/ &
1573 15526556 : 0.2e1_dp) + f2716*t979*t233)*alpha3*t185
1574 15526556 : t1009 = (t276*t444 - 2._dp*t197*t273*t444 + t243*t968)*alpha4
1575 15526556 : t1011 = t602*t444
1576 : t1024 = (t97*(-t243*t1011/0.2e1_dp + 0.3e1_dp/0.2e1_dp*t639*t250 &
1577 15526556 : *t444 - 0.5e1_dp/0.2e1_dp*t255*t645*t444) - t267*t979)*alpha5
1578 : t1035 = (-3._dp*t820*t444 + 2._dp*t823*t444 - t243*t826*t444 - &
1579 15526556 : t280*t968)*alpha6*t284
1580 : t1052 = (t97*(-0.7e1_dp/0.2e1_dp*t288*t691*t444 + 0.5e1_dp/0.2e1_dp &
1581 : *t697*t293*t444 - 0.3e1_dp/0.2e1_dp*t298*t701*t444 + t280 &
1582 15526556 : *t1011/0.2e1_dp) + t307*t979)*alpha7
1583 : t1066 = (-4._dp*t318*t713*t444 + 3._dp*t718*t444 - 2._dp*t243*t321 &
1584 15526556 : *t444 + t280*t275*t444 + t326*t968)*alpha8
1585 : t1068 = -t81*t93*t444 + t849*t94 - (2._dp*t91*t851) + t909 &
1586 : - (3._dp*t148*t910) + t937*t128 + t958*t190 - t909*t195 &
1587 : + (3._dp*t810*t962) + 0.5e1_dp/0.2e1_dp*t583*t965 + t197*(t968 &
1588 : + t972*t199) + t982*t214 + t988*t153 + t1001*t190 + t1009 &
1589 : *t247 + t1024*t271 + t1035*t286 + t1052*t316 + t1066* &
1590 15526556 : t333
1591 15526556 : t1069 = t1068*Clda
1592 15526556 : e_ndrho = e_ndrho + (-t80*t1069)*sx
1593 : END IF
1594 :
1595 24453786 : END SUBROUTINE xwpbe_lda_calc_1
1596 :
1597 : ! **************************************************************************************************
1598 : !> \brief Evaluates the screened hole averaged PBE exchange functional for lda.
1599 : !> \param e_0 ...
1600 : !> \param e_rho ...
1601 : !> \param e_ndrho ...
1602 : !> \param rho , ndrho: density and norm of the density gradient
1603 : !> \param ndrho ...
1604 : !> \param omega scaling factor
1605 : !> \param sscale scaling factor to enforce Lieb-Oxford bound
1606 : !> \param sx scaling factor
1607 : !> \param order maximum degree of the derivative that should be evaluated
1608 : !> \par History
1609 : !> 05.2007 created [Manuel Guidon]
1610 : !> \author Manuel Guidon
1611 : !> \note
1612 : !> This routine evaluates the functional for omega!=0 using a taylor
1613 : !> expansion for the parameter G.
1614 : ! **************************************************************************************************
1615 239869 : SUBROUTINE xwpbe_lda_calc_2(e_0, e_rho, e_ndrho, rho, ndrho, omega, sscale, sx, order)
1616 : REAL(KIND=dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho
1617 : REAL(KIND=dp), INTENT(IN) :: rho, ndrho, omega, sscale, sx
1618 : INTEGER, INTENT(IN) :: order
1619 :
1620 : REAL(KIND=dp) :: dQndrho, dQrho, Q, t1, t10, t100, t102, t103, t105, t106, t108, t109, t11, &
1621 : t110, t111, t112, t113, t114, t115, t116, t117, t118, t12, t121, t122, t125, t126, t127, &
1622 : t128, t129, t13, t133, t134, t136, t137, t14, t140, t141, t145, t146, t149, t15, t151, &
1623 : t153, t154, t155, t156, t158, t159, t16, t160, t161, t163, t166, t167, t17, t171, t172, &
1624 : t173, t176, t177, t18, t182, t184, t185, t186, t187, t19, t190, t191, t196, t2, t200, &
1625 : t202, t204, t206, t209, t21, t210, t213, t214, t215, t216, t218, t219, t22, t220, t221, &
1626 : t226, t227, t228, t23, t230, t233, t234, t236, t238, t239, t24
1627 : REAL(KIND=dp) :: t243, t247, t248, t249, t25, t251, t252, t253, t256, t257, t258, t261, &
1628 : t266, t268, t27, t270, t273, t274, t275, t276, t278, t279, t28, t281, t282, t284, t289, &
1629 : t29, t292, t293, t295, t296, t298, t299, t3, t301, t302, t303, t305, t306, t308, t309, &
1630 : t31, t310, t311, t312, t314, t315, t316, t317, t32, t320, t321, t324, t325, t326, t327, &
1631 : t331, t334, t335, t336, t337, t338, t339, t34, t340, t341, t346, t347, t348, t349, t35, &
1632 : t353, t355, t36, t361, t364, t365, t366, t367, t369, t372, t374, t375, t378, t379, t38, &
1633 : t382, t383, t384, t387, t388, t389, t39, t391, t392, t395, t396, t4
1634 : REAL(KIND=dp) :: t400, t403, t404, t405, t407, t409, t41, t413, t416, t417, t42, t420, t423, &
1635 : t428, t433, t434, t436, t44, t440, t441, t442, t443, t447, t448, t452, t453, t454, t457, &
1636 : t458, t46, t461, t464, t467, t468, t470, t474, t475, t476, t48, t480, t481, t482, t483, &
1637 : t487, t489, t49, t491, t496, t5, t500, t503, t505, t506, t507, t510, t512, t513, t514, &
1638 : t516, t519, t522, t523, t524, t530, t531, t535, t536, t54, t541, t542, t549, t55, t551, &
1639 : t552, t558, t559, t56, t561, t565, t566, t569, t574, t577, t579, t58, t583, t584, t585, &
1640 : t587, t588, t591, t595, t596, t6, t60, t603, t604, t605, t607
1641 : REAL(KIND=dp) :: t608, t61, t612, t615, t620, t622, t626, t628, t629, t63, t633, t634, t635, &
1642 : t638, t641, t644, t65, t650, t652, t656, t658, t659, t660, t665, t67, t670, t671, t679, &
1643 : t68, t681, t685, t687, t689, t69, t692, t695, t697, t699, t7, t70, t702, t705, t709, t71, &
1644 : t710, t714, t72, t723, t725, t727, t73, t733, t736, t737, t739, t74, t740, t742, t747, &
1645 : t748, t75, t751, t752, t755, t756, t758, t759, t760, t761, t767, t769, t77, t770, t78, &
1646 : t781, t783, t784, t788, t793, t796, t8, t80, t802, t805, t809, t81, t813, t816, t819, &
1647 : t82, t822, t823, t83, t830, t833, t839, t84, t85, t852, t860
1648 : REAL(KIND=dp) :: t862, t87, t875, t88, t886, t9, t90, t903, t91, t917, t919, t92, t920, t93, &
1649 : t94, t95, t96, t97, t99
1650 :
1651 239869 : t1 = ndrho**2
1652 239869 : t2 = r2**2
1653 239869 : t3 = 0.1e1_dp/t2
1654 239869 : t4 = t1*t3
1655 239869 : t5 = pi**2
1656 239869 : t6 = r3*t5
1657 239869 : t7 = t6*rho
1658 239869 : t8 = t7**(0.1e1_dp/0.3e1_dp)
1659 239869 : t9 = t8**2
1660 239869 : t10 = 0.1e1_dp/t9
1661 239869 : t11 = t4*t10
1662 239869 : t12 = rho**2
1663 239869 : t13 = 0.1e1_dp/t12
1664 239869 : t14 = sscale**2
1665 239869 : t15 = t13*t14
1666 239869 : t16 = a1*t1
1667 239869 : t17 = t16*t3
1668 239869 : t18 = t10*t13
1669 239869 : t19 = t18*t14
1670 239869 : t21 = t1**2
1671 239869 : t22 = a2*t21
1672 239869 : t23 = t2**2
1673 239869 : t24 = 0.1e1_dp/t23
1674 239869 : t25 = t22*t24
1675 239869 : t27 = 0.1e1_dp/t8/t7
1676 239869 : t28 = t12**2
1677 239869 : t29 = 0.1e1_dp/t28
1678 239869 : t31 = t14**2
1679 239869 : t32 = t27*t29*t31
1680 239869 : t34 = t17*t19 + t25*t32
1681 239869 : t35 = a3*t21
1682 239869 : t36 = t35*t24
1683 239869 : t38 = t21*ndrho
1684 239869 : t39 = a4*t38
1685 239869 : t41 = 0.1e1_dp/t23/r2
1686 239869 : t42 = t39*t41
1687 239869 : t44 = 0.1e1_dp/t9/t7
1688 239869 : t46 = 0.1e1_dp/t28/rho
1689 239869 : t48 = t31*sscale
1690 239869 : t49 = t44*t46*t48
1691 239869 : t54 = 0.1e1_dp/t23/t2
1692 239869 : t55 = a5*t21*t1*t54
1693 239869 : t56 = r3**2
1694 239869 : t58 = t5**2
1695 239869 : t60 = 0.1e1_dp/t56/t58
1696 239869 : t61 = t28**2
1697 239869 : t63 = t31*t14
1698 239869 : t65 = t60/t61*t63
1699 239869 : t67 = r1 + t36*t32 + t42*t49 + t55*t65
1700 239869 : t68 = 0.1e1_dp/t67
1701 239869 : t69 = t34*t68
1702 239869 : t70 = t15*t69
1703 239869 : t71 = t11*t70
1704 239869 : t72 = omega**2
1705 239869 : t73 = beta*t72
1706 239869 : t74 = t73*t10
1707 239869 : t75 = t71 + t74
1708 239869 : t77 = 0.1e1_dp/A
1709 239869 : Q = f94*t75*t77
1710 239869 : t78 = rho**(0.1e1_dp/0.3e1_dp)
1711 239869 : t80 = t78*rho*f89
1712 239869 : t81 = B*f12
1713 239869 : t82 = t71 + DD
1714 239869 : t83 = 0.1e1_dp/t82
1715 239869 : t84 = t81*t83
1716 239869 : t85 = F2*t34
1717 239869 : t87 = F1 + t85*t68
1718 239869 : t88 = t15*t87
1719 239869 : t90 = t11*t88 + r1
1720 239869 : t91 = f12*t90
1721 239869 : t92 = t82**2
1722 239869 : t93 = 0.1e1_dp/t92
1723 239869 : t94 = C*t93
1724 239869 : t95 = t91*t94
1725 239869 : t96 = g2*t1
1726 239869 : t97 = t96*t3
1727 239869 : t99 = g3*t21
1728 239869 : t100 = t99*t24
1729 239869 : t102 = g1 + t97*t19 + t100*t32
1730 239869 : t103 = t15*t102
1731 239869 : t105 = t11*t103 + r1
1732 239869 : t106 = t105*E
1733 239869 : t108 = 0.1e1_dp/t92/t82
1734 239869 : t109 = t106*t108
1735 239869 : t110 = f158*E
1736 239869 : t111 = t105*t83
1737 239869 : t112 = t72*t10
1738 239869 : t113 = t71 + DD + t112
1739 239869 : t114 = t113**2
1740 239869 : t115 = t114**2
1741 239869 : t116 = t115*t113
1742 239869 : t117 = SQRT(t116)
1743 239869 : t118 = 0.1e1_dp/t117
1744 239869 : t121 = SQRT(t113)
1745 239869 : t122 = 0.1e1_dp/t121
1746 239869 : t125 = f68*C
1747 239869 : t126 = t90*t83
1748 239869 : t127 = t114*t113
1749 239869 : t128 = SQRT(t127)
1750 239869 : t129 = 0.1e1_dp/t128
1751 : t133 = (-t110*t111*t118 - t81*t83*t122 - t125*t126*t129) &
1752 239869 : *omega
1753 239869 : t134 = 0.1e1_dp/t8
1754 239869 : t136 = f52*E
1755 239869 : t137 = t105*t93
1756 239869 : t140 = f12*C
1757 239869 : t141 = t90*t93
1758 239869 : t145 = t72*omega
1759 239869 : t146 = (-t136*t137*t118 - t140*t141*t129)*t145
1760 239869 : t149 = 0.1e1_dp/r3/t5
1761 239869 : t151 = t149/rho
1762 239869 : t153 = t72**2
1763 239869 : t154 = t153*omega
1764 239869 : t155 = t118*t154
1765 239869 : t156 = t155*t44
1766 239869 : t158 = f12*A
1767 239869 : t159 = exei(Q)
1768 239869 : t160 = t71 + DD + t74
1769 239869 : t161 = 0.1e1_dp/t160
1770 239869 : t163 = LOG(t75*t161)
1771 239869 : t166 = rootpi
1772 239869 : t167 = SQRT(t160)
1773 239869 : t171 = SQRT(A)
1774 239869 : t172 = t171*f34
1775 239869 : t173 = exer(Q)
1776 239869 : t176 = (t158*t166/t167 - t172*t173)*alpha1
1777 239869 : t177 = omega*t134
1778 239869 : t182 = (t158*t161 - f98*t159)*alpha2
1779 239869 : t184 = A*f14
1780 239869 : t185 = t160**2
1781 239869 : t186 = t185*t160
1782 239869 : t187 = SQRT(t186)
1783 239869 : t190 = SQRT(t75)
1784 239869 : t191 = 0.1e1_dp/t190
1785 239869 : t196 = 0.1e1_dp/t171
1786 : t200 = (t166*(t184/t187 - f98*t191) + f2716*t173*t196)* &
1787 239869 : alpha3*t145
1788 239869 : t202 = 0.1e1_dp/t75
1789 239869 : t204 = 0.1e1_dp/t185
1790 239869 : t206 = f8132*t77
1791 239869 : t209 = (-f98*t202 + t158*t204 + t206*t159)*alpha4
1792 239869 : t210 = t153*t27
1793 239869 : t213 = t75**2
1794 239869 : t214 = t213*t75
1795 239869 : t215 = SQRT(t214)
1796 239869 : t216 = 0.1e1_dp/t215
1797 239869 : t218 = f38*A
1798 239869 : t219 = t185**2
1799 239869 : t220 = t219*t160
1800 239869 : t221 = SQRT(t220)
1801 239869 : t226 = A**2
1802 239869 : t227 = t226*A
1803 239869 : t228 = SQRT(t227)
1804 239869 : t230 = f24364/t228
1805 : t233 = (t166*(t206*t191 - f916*t216 + t218/t221) - t230*t173) &
1806 239869 : *alpha5
1807 239869 : t234 = t154*t44
1808 239869 : t236 = 0.1e1_dp/t186
1809 239869 : t238 = 0.1e1_dp/t213
1810 239869 : t239 = f98*t238
1811 239869 : t243 = f729128/t226
1812 239869 : t247 = t153*t72
1813 : t248 = (A*t236 - t239 + t206*r1*t202 - t243*t159)*alpha6 &
1814 239869 : *t247
1815 239869 : t249 = t60*t13
1816 239869 : t251 = f1516*A
1817 239869 : t252 = t219*t186
1818 239869 : t253 = SQRT(t252)
1819 239869 : t256 = t213**2
1820 239869 : t257 = t256*t75
1821 239869 : t258 = SQRT(t257)
1822 239869 : t261 = f8164*t77
1823 239869 : t266 = t226**2
1824 239869 : t268 = SQRT(t266*A)
1825 239869 : t270 = f2187256/t268
1826 : t273 = (t166*(t251/t253 - f2732/t258 + t261*t216 - t243*t191) &
1827 239869 : + t270*t173)*alpha7
1828 239869 : t274 = t153*t145
1829 239869 : t275 = t56*t58
1830 239869 : t276 = t275*t12
1831 239869 : t278 = 0.1e1_dp/t8/t276
1832 239869 : t279 = t274*t278
1833 239869 : t281 = r3*A
1834 239869 : t282 = 0.1e1_dp/t219
1835 239869 : t284 = 0.1e1_dp/t214
1836 239869 : t289 = f6561512/t227
1837 : t292 = (t281*t282 - f94*t284 + t206*t238 - t243*t202 + t289 &
1838 239869 : *t159)*alpha8
1839 239869 : t293 = t153**2
1840 239869 : t295 = 0.1e1_dp/t9/t276
1841 239869 : t296 = t293*t295
1842 : t298 = t84 + t95 + t109 + t133*t134 + t146*t151 - t109*t156 + &
1843 : t158*(t159 + t163) + t176*t177 + t182*t112 + t200*t151 + &
1844 239869 : t209*t210 + t233*t234 + t248*t249 + t273*t279 + t292*t296
1845 239869 : t299 = t298*Clda
1846 239869 : e_0 = e_0 + (-t80*t299)*sx
1847 239869 : IF (order >= 1) THEN
1848 236157 : t301 = t44*t13
1849 236157 : t302 = t4*t301
1850 236157 : t303 = t14*t34
1851 236157 : t305 = t68*r3*t5
1852 236157 : t306 = t303*t305
1853 236157 : t308 = 0.2e1_dp/0.3e1_dp*t302*t306
1854 236157 : t309 = t12*rho
1855 236157 : t310 = 0.1e1_dp/t309
1856 236157 : t311 = t310*t14
1857 236157 : t312 = t311*t69
1858 236157 : t314 = 2._dp*t11*t312
1859 236157 : t315 = t3*t44
1860 236157 : t316 = t16*t315
1861 236157 : t317 = t15*t6
1862 236157 : t320 = t10*t310
1863 236157 : t321 = t320*t14
1864 236157 : t324 = t24*t278
1865 236157 : t325 = t22*t324
1866 236157 : t326 = t29*t31
1867 236157 : t327 = t326*t6
1868 236157 : t331 = t27*t46*t31
1869 : t334 = -0.2e1_dp/0.3e1_dp*t316*t317 - (2._dp*t17*t321) - 0.4e1_dp &
1870 236157 : /0.3e1_dp*t325*t327 - (4._dp*t25*t331)
1871 236157 : t335 = t334*t68
1872 236157 : t336 = t15*t335
1873 236157 : t337 = t11*t336
1874 236157 : t338 = t4*t18
1875 236157 : t339 = t67**2
1876 236157 : t340 = 0.1e1_dp/t339
1877 236157 : t341 = t35*t324
1878 236157 : t346 = t41*t295
1879 236157 : t347 = t39*t346
1880 236157 : t348 = t46*t48
1881 236157 : t349 = t348*t6
1882 236157 : t353 = 0.1e1_dp/t28/t12
1883 236157 : t355 = t44*t353*t48
1884 236157 : t361 = t60/t61/rho*t63
1885 : t364 = -0.4e1_dp/0.3e1_dp*t341*t327 - (4._dp*t36*t331) - 0.5e1_dp &
1886 236157 : /0.3e1_dp*t347*t349 - (5._dp*t42*t355) - (8._dp*t55*t361)
1887 236157 : t365 = t340*t364
1888 236157 : t366 = t303*t365
1889 236157 : t367 = t338*t366
1890 236157 : t369 = t44*r3*t5
1891 236157 : t372 = -t308 - t314 + t337 - t367 - 0.2e1_dp/0.3e1_dp*t73*t369
1892 236157 : dQrho = f94*t372*t77
1893 236157 : t374 = ndrho*t3
1894 236157 : t375 = t374*t10
1895 236157 : t378 = a1*ndrho
1896 236157 : t379 = t378*t3
1897 236157 : t382 = t1*ndrho
1898 236157 : t383 = a2*t382
1899 236157 : t384 = t383*t24
1900 236157 : t387 = 2._dp*t379*t19 + 4._dp*t384*t32
1901 236157 : t388 = t387*t68
1902 236157 : t389 = t15*t388
1903 236157 : t391 = a3*t382
1904 236157 : t392 = t391*t24
1905 236157 : t395 = a4*t21
1906 236157 : t396 = t395*t41
1907 236157 : t400 = a5*t38*t54
1908 236157 : t403 = 4._dp*t392*t32 + 5._dp*t396*t49 + 6._dp*t400*t65
1909 236157 : t404 = t340*t403
1910 236157 : t405 = t303*t404
1911 236157 : t407 = 2._dp*t375*t70 + t11*t389 - t338*t405
1912 236157 : dQndrho = f94*t407*t77
1913 236157 : t409 = t78*f89
1914 236157 : t413 = t27*r3*t5
1915 236157 : t416 = t14*t102
1916 236157 : t417 = t416*t6
1917 236157 : t420 = t311*t102
1918 236157 : t423 = t96*t315
1919 236157 : t428 = t99*t324
1920 : t433 = -0.2e1_dp/0.3e1_dp*t423*t317 - (2._dp*t97*t321) - 0.4e1_dp &
1921 236157 : /0.3e1_dp*t428*t327 - (4._dp*t100*t331)
1922 236157 : t434 = t15*t433
1923 : t436 = -0.2e1_dp/0.3e1_dp*t302*t417 - (2._dp*t11*t420) + (t11 &
1924 236157 : *t434)
1925 236157 : t440 = t136*t105
1926 236157 : t441 = t108*t118
1927 236157 : t442 = -t308 - t314 + t337 - t367
1928 236157 : t443 = t441*t442
1929 236157 : t447 = 0.1e1_dp/t117/t116
1930 236157 : t448 = t93*t447
1931 236157 : t452 = -t308 - t314 + t337 - t367 - 0.2e1_dp/0.3e1_dp*t72*t44*t6
1932 236157 : t453 = t115*t452
1933 236157 : t454 = t448*t453
1934 236157 : t457 = t14*t87
1935 236157 : t458 = t457*t6
1936 236157 : t461 = t311*t87
1937 236157 : t464 = F2*t334
1938 236157 : t467 = t464*t68 - t85*t365
1939 236157 : t468 = t15*t467
1940 : t470 = -0.2e1_dp/0.3e1_dp*t302*t458 - (2._dp*t11*t461) + (t11 &
1941 236157 : *t468)
1942 236157 : t474 = t140*t90
1943 236157 : t475 = t108*t129
1944 236157 : t476 = t475*t442
1945 236157 : t480 = 0.1e1_dp/t128/t127
1946 236157 : t481 = t93*t480
1947 236157 : t482 = t114*t452
1948 236157 : t483 = t481*t482
1949 : t487 = (-t136*t436*t93*t118 + (2._dp*t440*t443) + 0.5e1_dp/ &
1950 : 0.2e1_dp*(t440)*(t454) - t140*t470*t93*t129 + (2._dp &
1951 236157 : *t474*t476) + 0.3e1_dp/0.2e1_dp*(t474)*(t483))*t145
1952 236157 : t489 = t209*t153
1953 236157 : t491 = t278*r3*t5
1954 236157 : t496 = t166/t167/t160
1955 236157 : t500 = dexerrho(Q, dQrho)
1956 236157 : t503 = (-t158*t496*t372/0.2e1_dp - t172*t500)*alpha1
1957 236157 : t505 = t106*t441
1958 236157 : t506 = t154*t295
1959 236157 : t507 = t506*t6
1960 236157 : t510 = f12*t470
1961 236157 : t512 = t92**2
1962 236157 : t513 = 0.1e1_dp/t512
1963 236157 : t514 = t106*t513
1964 236157 : t516 = t155*t44*t442
1965 236157 : t519 = t149*t13
1966 236157 : t522 = 0.1e1_dp/t253/t252
1967 236157 : t523 = t219*t185
1968 236157 : t524 = t522*t523
1969 236157 : t530 = f2732/t258/t257
1970 236157 : t531 = t256*t372
1971 236157 : t535 = 0.1e1_dp/t215/t214
1972 236157 : t536 = t535*t213
1973 236157 : t541 = 0.1e1_dp/t190/t75
1974 236157 : t542 = t541*t372
1975 : t549 = (t166*(-0.7e1_dp/0.2e1_dp*t251*t524*t372 + 0.5e1_dp/0.2e1_dp &
1976 : *t530*t531 - 0.3e1_dp/0.2e1_dp*t261*t536*t372 + t243*t542 &
1977 236157 : /0.2e1_dp) + t270*t500)*alpha7
1978 236157 : t551 = C*t108
1979 236157 : t552 = t551*t442
1980 236157 : t558 = t436*E
1981 236157 : t559 = t558*t108
1982 236157 : t561 = 0.1e1_dp/t220
1983 236157 : t565 = 0.1e1_dp/t256
1984 236157 : t566 = f94*t565
1985 236157 : t569 = t284*t372
1986 236157 : t574 = dexeirho(Q, dQrho)
1987 : t577 = (-4._dp*t281*t561*t372 + 3._dp*t566*t372 - 2._dp*t206*t569 &
1988 236157 : + t243*t238*t372 + t289*t574)*alpha8
1989 : t579 = -t133*t413/0.3e1_dp + t487*t151 - 0.4e1_dp/0.3e1_dp*t489* &
1990 : t491 + t503*t177 + 0.5e1_dp/0.3e1_dp*t505*t507 + t510*t94 + (3._dp &
1991 : *t514*t516) - t146*t519 + t549*t279 - (2._dp*t91*t552) &
1992 : - t81*t93*t442 - t200*t519 - t559*t156 + t577*t296 &
1993 236157 : + t559
1994 236157 : t583 = t110*t105
1995 236157 : t584 = t93*t118
1996 236157 : t585 = t584*t442
1997 236157 : t587 = t83*t447
1998 236157 : t588 = t587*t453
1999 236157 : t591 = t93*t122
2000 236157 : t595 = 0.1e1_dp/t121/t113
2001 236157 : t596 = t83*t595
2002 236157 : t603 = t125*t90
2003 236157 : t604 = t93*t129
2004 236157 : t605 = t604*t442
2005 236157 : t607 = t83*t480
2006 236157 : t608 = t607*t482
2007 : t612 = (-t110*t436*t83*t118 + t583*t585 + 0.5e1_dp/0.2e1_dp*t583 &
2008 : *t588 + t81*t591*t442 + t81*t596*t452/0.2e1_dp - t125 &
2009 : *t470*t83*t129 + t603*t605 + 0.3e1_dp/0.2e1_dp*t603*t608)* &
2010 236157 : omega
2011 236157 : t615 = t236*t372
2012 236157 : t620 = (t239*t372 - 2._dp*t158*t615 + t206*t574)*alpha4
2013 236157 : t622 = t513*t442
2014 236157 : t626 = t75*t204
2015 236157 : t628 = t372*t161 - t626*t372
2016 236157 : t629 = t628*t202
2017 236157 : t633 = t233*t154
2018 236157 : t634 = t295*r3
2019 236157 : t635 = t634*t5
2020 236157 : t638 = A*t282
2021 236157 : t641 = f98*t284
2022 236157 : t644 = r1*t238
2023 : t650 = (-3._dp*t638*t372 + 2._dp*t641*t372 - t206*t644*t372 - t243 &
2024 236157 : *t574)*alpha6*t247
2025 236157 : t652 = t204*t372
2026 236157 : t656 = (-t158*t652 - f98*t574)*alpha2
2027 236157 : t658 = t108*t447
2028 236157 : t659 = t106*t658
2029 236157 : t660 = t234*t453
2030 236157 : t665 = f916*t535
2031 236157 : t670 = 0.1e1_dp/t221/t220
2032 236157 : t671 = t670*t219
2033 : t679 = (t166*(-t206*t542/0.2e1_dp + 0.3e1_dp/0.2e1_dp*t665*t213 &
2034 236157 : *t372 - 0.5e1_dp/0.2e1_dp*t218*t671*t372) - t230*t500)*alpha5
2035 236157 : t681 = t292*t293
2036 236157 : t685 = t56*r3*t58*t5*t309
2037 236157 : t687 = 0.1e1_dp/t9/t685
2038 236157 : t689 = t687*r3*t5
2039 236157 : t692 = t60*t310
2040 236157 : t695 = t273*t274
2041 236157 : t697 = 0.1e1_dp/t8/t685
2042 236157 : t699 = t697*r3*t5
2043 236157 : t702 = t176*omega
2044 236157 : t705 = t182*t72
2045 236157 : t709 = 0.1e1_dp/t187/t186
2046 236157 : t710 = t709*t185
2047 236157 : t714 = f98*t541
2048 : t723 = (t166*(-0.3e1_dp/0.2e1_dp*t184*t710*t372 + t714*t372/ &
2049 236157 : 0.2e1_dp) + f2716*t500*t196)*alpha3*t145
2050 : t725 = t612*t134 + t620*t210 - (3._dp*t106*t622) + t158*(t574 &
2051 : + t629*t160) - 0.5e1_dp/0.3e1_dp*t633*t635 + t650*t249 + &
2052 : t656*t112 + 0.5e1_dp/0.2e1_dp*t659*t660 + t679*t234 - 0.8e1_dp/ &
2053 : 0.3e1_dp*t681*t689 - (2._dp*t248*t692) - 0.7e1_dp/0.3e1_dp*t695 &
2054 : *t699 - t702*t413/0.3e1_dp - 0.2e1_dp/0.3e1_dp*t705*t369 + t723 &
2055 236157 : *t151
2056 236157 : t727 = (t579 + t725)*Clda
2057 236157 : e_rho = e_rho + (-0.4e1_dp/0.3e1_dp*t409*t299 - t80*t727)*sx
2058 236157 : t733 = F2*t387
2059 236157 : t736 = t733*t68 - t85*t404
2060 236157 : t737 = t15*t736
2061 236157 : t739 = 2._dp*t375*t88 + t11*t737
2062 236157 : t740 = f12*t739
2063 236157 : t742 = t551*t407
2064 236157 : t747 = g2*ndrho
2065 236157 : t748 = t747*t3
2066 236157 : t751 = g3*t382
2067 236157 : t752 = t751*t24
2068 236157 : t755 = 2._dp*t748*t19 + 4._dp*t752*t32
2069 236157 : t756 = t15*t755
2070 236157 : t758 = 2._dp*t375*t103 + t11*t756
2071 236157 : t759 = t758*E
2072 236157 : t760 = t759*t108
2073 236157 : t761 = t513*t407
2074 236157 : t767 = t584*t407
2075 236157 : t769 = t115*t407
2076 236157 : t770 = t587*t769
2077 236157 : t781 = t604*t407
2078 236157 : t783 = t114*t407
2079 236157 : t784 = t607*t783
2080 : t788 = (-t110*t758*t83*t118 + t583*t767 + 0.5e1_dp/0.2e1_dp*t583 &
2081 : *t770 + t81*t591*t407 + t81*t596*t407/0.2e1_dp - t125 &
2082 : *t739*t83*t129 + t603*t781 + 0.3e1_dp/0.2e1_dp*t603*t784)* &
2083 236157 : omega
2084 236157 : t793 = t441*t407
2085 236157 : t796 = t448*t769
2086 236157 : t802 = t475*t407
2087 236157 : t805 = t481*t783
2088 : t809 = (-t136*t758*t93*t118 + (2._dp*t440*t793) + 0.5e1_dp/ &
2089 : 0.2e1_dp*(t440)*(t796) - t140*t739*t93*t129 + (2._dp &
2090 236157 : *t474*t802) + 0.3e1_dp/0.2e1_dp*(t474)*(t805))*t145
2091 236157 : t813 = t155*t44*t407
2092 236157 : t816 = t234*t769
2093 236157 : t819 = dexeindrho(Q, dQndrho)
2094 236157 : t822 = t407*t161 - t626*t407
2095 236157 : t823 = t822*t202
2096 236157 : t830 = dexerndrho(Q, dQndrho)
2097 236157 : t833 = (-t158*t496*t407/0.2e1_dp - t172*t830)*alpha1
2098 236157 : t839 = (-t158*t204*t407 - f98*t819)*alpha2
2099 : t852 = (t166*(-0.3e1_dp/0.2e1_dp*t184*t710*t407 + t714*t407/ &
2100 236157 : 0.2e1_dp) + f2716*t830*t196)*alpha3*t145
2101 236157 : t860 = (t239*t407 - 2._dp*t158*t236*t407 + t206*t819)*alpha4
2102 236157 : t862 = t541*t407
2103 : t875 = (t166*(-t206*t862/0.2e1_dp + 0.3e1_dp/0.2e1_dp*t665*t213 &
2104 236157 : *t407 - 0.5e1_dp/0.2e1_dp*t218*t671*t407) - t230*t830)*alpha5
2105 : t886 = (-3._dp*t638*t407 + 2._dp*t641*t407 - t206*t644*t407 - t243 &
2106 236157 : *t819)*alpha6*t247
2107 : t903 = (t166*(-0.7e1_dp/0.2e1_dp*t251*t524*t407 + 0.5e1_dp/0.2e1_dp &
2108 : *t530*t256*t407 - 0.3e1_dp/0.2e1_dp*t261*t536*t407 + t243 &
2109 236157 : *t862/0.2e1_dp) + t270*t830)*alpha7
2110 : t917 = (-4._dp*t281*t561*t407 + 3._dp*t566*t407 - 2._dp*t206*t284 &
2111 236157 : *t407 + t243*t238*t407 + t289*t819)*alpha8
2112 : t919 = -t81*t93*t407 + t740*t94 - (2._dp*t91*t742) + t760 &
2113 : - (3._dp*t106*t761) + t788*t134 + t809*t151 - t760*t156 &
2114 : + (3._dp*t514*t813) + 0.5e1_dp/0.2e1_dp*t659*t816 + t158*(t819 &
2115 : + t823*t160) + t833*t177 + t839*t112 + t852*t151 + t860 &
2116 236157 : *t210 + t875*t234 + t886*t249 + t903*t279 + t917*t296
2117 236157 : t920 = t919*Clda
2118 236157 : e_ndrho = e_ndrho + (-t80*t920)*sx
2119 : END IF
2120 :
2121 239869 : END SUBROUTINE xwpbe_lda_calc_2
2122 :
2123 : ! **************************************************************************************************
2124 : !> \brief Evaluates the screened hole averaged PBE exchange functional for lda.
2125 : !> \param e_0 ...
2126 : !> \param e_rho ...
2127 : !> \param e_ndrho ...
2128 : !> \param rho , ndrho: density and norm of the density gradient
2129 : !> \param ndrho ...
2130 : !> \param omega screening parameter
2131 : !> \param sscale scaling factor to enforce Lieb-Oxford bound
2132 : !> \param sx scaling factor
2133 : !> \param order maximum degree of the derivative that should be evaluated
2134 : !> \par History
2135 : !> 05.2007 created [Manuel Guidon]
2136 : !> \author Manuel Guidon
2137 : !> \note
2138 : !> This routine evaluates the functional for omega!=0 using a simple
2139 : !> gaussian expansion for large ww.
2140 : ! **************************************************************************************************
2141 7448304 : SUBROUTINE xwpbe_lda_calc_3(e_0, e_rho, e_ndrho, rho, ndrho, omega, sscale, sx, order)
2142 : REAL(KIND=dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho
2143 : REAL(KIND=dp), INTENT(IN) :: rho, ndrho, omega, sscale, sx
2144 : INTEGER, INTENT(IN) :: order
2145 :
2146 : REAL(KIND=dp) :: dQndrho, dQrho, Q, t1, t10, t102, t104, t105, t108, t109, t11, t110, t111, &
2147 : t113, t115, t116, t117, t118, t119, t12, t120, t121, t123, t124, t125, t126, t127, t128, &
2148 : t129, t13, t131, t132, t133, t136, t14, t140, t141, t142, t144, t145, t147, t148, t149, &
2149 : t15, t150, t151, t152, t154, t155, t156, t157, t158, t159, t16, t162, t163, t166, t167, &
2150 : t168, t169, t17, t170, t174, t176, t177, t18, t180, t181, t185, t186, t189, t19, t190, &
2151 : t192, t193, t194, t195, t197, t198, t199, t2, t200, t202, t206, t208, t209, t21, t210, &
2152 : t212, t213, t215, t216, t217, t218, t219, t22, t221, t222, t223
2153 : REAL(KIND=dp) :: t224, t227, t228, t23, t231, t232, t234, t235, t236, t237, t238, t24, t242, &
2154 : t245, t246, t247, t248, t249, t25, t250, t251, t252, t258, t259, t260, t261, t262, t266, &
2155 : t268, t27, t274, t277, t278, t279, t28, t280, t285, t287, t288, t29, t291, t292, t295, &
2156 : t296, t297, t3, t300, t301, t302, t304, t305, t308, t309, t31, t313, t316, t317, t318, &
2157 : t32, t320, t322, t325, t328, t329, t332, t335, t338, t339, t34, t341, t342, t344, t345, &
2158 : t35, t352, t359, t36, t361, t364, t365, t368, t369, t370, t373, t374, t375, t376, t377, &
2159 : t379, t38, t380, t381, t383, t384, t388, t39, t391, t395, t396
2160 : REAL(KIND=dp) :: t397, t398, t399, t4, t400, t401, t403, t404, t405, t406, t408, t41, t416, &
2161 : t417, t418, t419, t42, t420, t422, t423, t424, t428, t429, t433, t435, t437, t438, t44, &
2162 : t441, t442, t443, t444, t445, t451, t452, t453, t456, t457, t46, t461, t462, t463, t466, &
2163 : t470, t471, t478, t479, t48, t480, t483, t484, t485, t486, t49, t490, t493, t499, t5, &
2164 : t500, t501, t504, t505, t511, t512, t513, t516, t517, t521, t523, t526, t528, t531, t532, &
2165 : t533, t534, t537, t538, t539, t54, t542, t544, t545, t546, t548, t549, t55, t550, t554, &
2166 : t555, t56, t561, t564, t565, t567, t568, t570, t58, t584, t586
2167 : REAL(KIND=dp) :: t589, t590, t592, t593, t595, t596, t599, t6, t60, t603, t604, t607, t608, &
2168 : t61, t610, t611, t612, t616, t617, t621, t623, t626, t627, t628, t629, t63, t635, t637, &
2169 : t638, t649, t65, t651, t652, t656, t661, t664, t67, t670, t673, t677, t68, t681, t684, &
2170 : t687, t69, t690, t691, t695, t696, t7, t70, t71, t72, t73, t74, t75, t77, t78, t8, t80, &
2171 : t81, t82, t83, t84, t85, t87, t88, t9, t90, t91, t92, t93, t94, t95, t96, t97, t99
2172 :
2173 7448304 : t1 = ndrho**2
2174 7448304 : t2 = r2**2
2175 7448304 : t3 = 0.1e1_dp/t2
2176 7448304 : t4 = t1*t3
2177 7448304 : t5 = pi**2
2178 7448304 : t6 = r3*t5
2179 7448304 : t7 = t6*rho
2180 7448304 : t8 = t7**(0.1e1_dp/0.3e1_dp)
2181 7448304 : t9 = t8**2
2182 7448304 : t10 = 0.1e1_dp/t9
2183 7448304 : t11 = t4*t10
2184 7448304 : t12 = rho**2
2185 7448304 : t13 = 0.1e1_dp/t12
2186 7448304 : t14 = sscale**2
2187 7448304 : t15 = t13*t14
2188 7448304 : t16 = a1*t1
2189 7448304 : t17 = t16*t3
2190 7448304 : t18 = t10*t13
2191 7448304 : t19 = t18*t14
2192 7448304 : t21 = t1**2
2193 7448304 : t22 = a2*t21
2194 7448304 : t23 = t2**2
2195 7448304 : t24 = 0.1e1_dp/t23
2196 7448304 : t25 = t22*t24
2197 7448304 : t27 = 0.1e1_dp/t8/t7
2198 7448304 : t28 = t12**2
2199 7448304 : t29 = 0.1e1_dp/t28
2200 7448304 : t31 = t14**2
2201 7448304 : t32 = t27*t29*t31
2202 7448304 : t34 = t17*t19 + t25*t32
2203 7448304 : t35 = a3*t21
2204 7448304 : t36 = t35*t24
2205 7448304 : t38 = t21*ndrho
2206 7448304 : t39 = a4*t38
2207 7448304 : t41 = 0.1e1_dp/t23/r2
2208 7448304 : t42 = t39*t41
2209 7448304 : t44 = 0.1e1_dp/t9/t7
2210 7448304 : t46 = 0.1e1_dp/t28/rho
2211 7448304 : t48 = t31*sscale
2212 7448304 : t49 = t44*t46*t48
2213 7448304 : t54 = 0.1e1_dp/t23/t2
2214 7448304 : t55 = a5*t21*t1*t54
2215 7448304 : t56 = r3**2
2216 7448304 : t58 = t5**2
2217 7448304 : t60 = 0.1e1_dp/t56/t58
2218 7448304 : t61 = t28**2
2219 7448304 : t63 = t31*t14
2220 7448304 : t65 = t60/t61*t63
2221 7448304 : t67 = r1 + t36*t32 + t42*t49 + t55*t65
2222 7448304 : t68 = 0.1e1_dp/t67
2223 7448304 : t69 = t34*t68
2224 7448304 : t70 = t15*t69
2225 7448304 : t71 = t11*t70
2226 7448304 : t72 = omega**2
2227 7448304 : t73 = beta2*t72
2228 7448304 : t74 = t73*t10
2229 7448304 : t75 = t71 + t74
2230 7448304 : t77 = 0.1e1_dp/A
2231 7448304 : Q = f94*t75*t77
2232 7448304 : t78 = rho**(0.1e1_dp/0.3e1_dp)
2233 7448304 : t80 = t78*rho*f89
2234 7448304 : t81 = B*f12
2235 7448304 : t82 = t71 + DD
2236 7448304 : t83 = 0.1e1_dp/t82
2237 7448304 : t84 = t81*t83
2238 7448304 : t85 = F2*t34
2239 7448304 : t87 = F1 + t85*t68
2240 7448304 : t88 = t15*t87
2241 7448304 : t90 = t11*t88 + r1
2242 7448304 : t91 = f12*t90
2243 7448304 : t92 = t82**2
2244 7448304 : t93 = 0.1e1_dp/t92
2245 7448304 : t94 = C*t93
2246 7448304 : t95 = t91*t94
2247 7448304 : t96 = f34*pi
2248 7448304 : t97 = rootpi
2249 7448304 : t99 = r6*C
2250 7448304 : t102 = r4*B
2251 7448304 : t104 = r8*A
2252 7448304 : t105 = t92*t82
2253 7448304 : t108 = t97*(r15*E + t99*t90*t82 + t102*t92 + t104*t105)
2254 7448304 : t109 = 0.1e1_dp/r16
2255 7448304 : t110 = SQRT(t82)
2256 7448304 : t111 = t110*t105
2257 7448304 : t113 = t109/t111
2258 7448304 : t115 = SQRT(A)
2259 7448304 : t116 = f94*t34
2260 7448304 : t117 = t68*t1
2261 7448304 : t118 = t116*t117
2262 7448304 : t119 = t3*t10
2263 7448304 : t120 = t15*t77
2264 7448304 : t121 = t119*t120
2265 7448304 : t123 = EXP(t118*t121)
2266 7448304 : t124 = t115*t123
2267 7448304 : t125 = f32*ndrho
2268 7448304 : t126 = 0.1e1_dp/r2
2269 7448304 : t127 = t125*t126
2270 7448304 : t128 = 0.1e1_dp/t8
2271 7448304 : t129 = 0.1e1_dp/rho
2272 7448304 : t131 = t69*t77
2273 7448304 : t132 = SQRT(t131)
2274 7448304 : t133 = sscale*t132
2275 7448304 : t136 = erfc(t127*t128*t129*t133)
2276 7448304 : t140 = 0.1e1_dp/f1516
2277 7448304 : t141 = (t96 + t108*t113 - t96*t124*t136)*t140
2278 7448304 : t142 = 0.1e1_dp/t97
2279 7448304 : t144 = 0.1e1_dp/E
2280 7448304 : t145 = t142*t111*t144
2281 7448304 : t147 = -t141*t145 + r1
2282 7448304 : t148 = t147*E
2283 7448304 : t149 = 0.1e1_dp/t105
2284 7448304 : t150 = t148*t149
2285 7448304 : t151 = f158*E
2286 7448304 : t152 = t147*t83
2287 7448304 : t154 = t71 + DD + t72*t10
2288 7448304 : t155 = t154**2
2289 7448304 : t156 = t155**2
2290 7448304 : t157 = t156*t154
2291 7448304 : t158 = SQRT(t157)
2292 7448304 : t159 = 0.1e1_dp/t158
2293 7448304 : t162 = SQRT(t154)
2294 7448304 : t163 = 0.1e1_dp/t162
2295 7448304 : t166 = f68*C
2296 7448304 : t167 = t90*t83
2297 7448304 : t168 = t155*t154
2298 7448304 : t169 = SQRT(t168)
2299 7448304 : t170 = 0.1e1_dp/t169
2300 : t174 = (-t151*t152*t159 - t81*t83*t163 - t166*t167*t170) &
2301 7448304 : *omega
2302 7448304 : t176 = f52*E
2303 7448304 : t177 = t147*t93
2304 7448304 : t180 = f12*C
2305 7448304 : t181 = t90*t93
2306 7448304 : t185 = t72*omega
2307 7448304 : t186 = (-t176*t177*t159 - t180*t181*t170)*t185
2308 7448304 : t189 = 0.1e1_dp/r3/t5
2309 7448304 : t190 = t189*t129
2310 7448304 : t192 = t72**2
2311 7448304 : t193 = t192*omega
2312 7448304 : t194 = t159*t193
2313 7448304 : t195 = t194*t44
2314 7448304 : t197 = f12*A
2315 7448304 : t198 = exei(Q)
2316 7448304 : t199 = t71 + DD + t74
2317 7448304 : t200 = 0.1e1_dp/t199
2318 7448304 : t202 = LOG(t75*t200)
2319 : t206 = (t84 + t95 + t150 + t174*t128 + t186*t190 - t150*t195 &
2320 7448304 : + t197*(t198 + t202))*Clda
2321 7448304 : e_0 = e_0 + (-t80*t206)*sx
2322 7448304 : IF (order >= 1) THEN
2323 4713253 : t208 = t44*t13
2324 4713253 : t209 = t4*t208
2325 4713253 : t210 = t14*t34
2326 4713253 : t212 = t68*r3*t5
2327 4713253 : t213 = t210*t212
2328 4713253 : t215 = 0.2e1_dp/0.3e1_dp*t209*t213
2329 4713253 : t216 = t12*rho
2330 4713253 : t217 = 0.1e1_dp/t216
2331 4713253 : t218 = t217*t14
2332 4713253 : t219 = t218*t69
2333 4713253 : t221 = 2._dp*t11*t219
2334 4713253 : t222 = t3*t44
2335 4713253 : t223 = t16*t222
2336 4713253 : t224 = t15*t6
2337 4713253 : t227 = t10*t217
2338 4713253 : t228 = t227*t14
2339 4713253 : t231 = t56*t58
2340 4713253 : t232 = t231*t12
2341 4713253 : t234 = 0.1e1_dp/t8/t232
2342 4713253 : t235 = t24*t234
2343 4713253 : t236 = t22*t235
2344 4713253 : t237 = t29*t31
2345 4713253 : t238 = t237*t6
2346 4713253 : t242 = t27*t46*t31
2347 : t245 = -0.2e1_dp/0.3e1_dp*t223*t224 - (2._dp*t17*t228) - 0.4e1_dp &
2348 4713253 : /0.3e1_dp*t236*t238 - (4._dp*t25*t242)
2349 4713253 : t246 = t245*t68
2350 4713253 : t247 = t15*t246
2351 4713253 : t248 = t11*t247
2352 4713253 : t249 = t4*t18
2353 4713253 : t250 = t67**2
2354 4713253 : t251 = 0.1e1_dp/t250
2355 4713253 : t252 = t35*t235
2356 4713253 : t258 = 0.1e1_dp/t9/t232
2357 4713253 : t259 = t41*t258
2358 4713253 : t260 = t39*t259
2359 4713253 : t261 = t46*t48
2360 4713253 : t262 = t261*t6
2361 4713253 : t266 = 0.1e1_dp/t28/t12
2362 4713253 : t268 = t44*t266*t48
2363 4713253 : t274 = t60/t61/rho*t63
2364 : t277 = -0.4e1_dp/0.3e1_dp*t252*t238 - (4._dp*t36*t242) - 0.5e1_dp &
2365 4713253 : /0.3e1_dp*t260*t262 - (5._dp*t42*t268) - (8._dp*t55*t274)
2366 4713253 : t278 = t251*t277
2367 4713253 : t279 = t210*t278
2368 4713253 : t280 = t249*t279
2369 : t285 = -t215 - t221 + t248 - t280 - 0.2e1_dp/0.3e1_dp*t73*t44*r3 &
2370 4713253 : *t5
2371 4713253 : dQrho = f94*t285*t77
2372 4713253 : t287 = ndrho*t3
2373 4713253 : t288 = t287*t10
2374 4713253 : t291 = a1*ndrho
2375 4713253 : t292 = t291*t3
2376 4713253 : t295 = t1*ndrho
2377 4713253 : t296 = a2*t295
2378 4713253 : t297 = t296*t24
2379 4713253 : t300 = 2._dp*t292*t19 + 4._dp*t297*t32
2380 4713253 : t301 = t300*t68
2381 4713253 : t302 = t15*t301
2382 4713253 : t304 = a3*t295
2383 4713253 : t305 = t304*t24
2384 4713253 : t308 = a4*t21
2385 4713253 : t309 = t308*t41
2386 4713253 : t313 = a5*t38*t54
2387 4713253 : t316 = 4._dp*t305*t32 + 5._dp*t309*t49 + 6._dp*t313*t65
2388 4713253 : t317 = t251*t316
2389 4713253 : t318 = t210*t317
2390 4713253 : t320 = 2._dp*t288*t70 + t11*t302 - t249*t318
2391 4713253 : dQndrho = f94*t320*t77
2392 4713253 : t322 = t78*f89
2393 4713253 : t325 = -t215 - t221 + t248 - t280
2394 4713253 : t328 = t14*t87
2395 4713253 : t329 = t328*t6
2396 4713253 : t332 = t218*t87
2397 4713253 : t335 = F2*t245
2398 4713253 : t338 = t335*t68 - t85*t278
2399 4713253 : t339 = t15*t338
2400 : t341 = -0.2e1_dp/0.3e1_dp*t209*t329 - (2._dp*t11*t332) + (t11 &
2401 4713253 : *t339)
2402 4713253 : t342 = f12*t341
2403 4713253 : t344 = C*t149
2404 4713253 : t345 = t344*t325
2405 4713253 : t352 = t82*t325
2406 : t359 = t97*(t99*t341*t82 + t99*t90*t325 + 2._dp*t102*t352 &
2407 4713253 : + 3._dp*t104*t92*t325)
2408 4713253 : t361 = t92**2
2409 4713253 : t364 = t109/t110/t361
2410 4713253 : t365 = t364*t325
2411 4713253 : t368 = t96*t115
2412 4713253 : t369 = f94*t245
2413 4713253 : t370 = t369*t117
2414 4713253 : t373 = t251*t1*t3
2415 4713253 : t374 = t116*t373
2416 4713253 : t375 = t14*t77
2417 4713253 : t376 = t375*t277
2418 4713253 : t377 = t18*t376
2419 4713253 : t379 = t117*t3
2420 4713253 : t380 = t116*t379
2421 4713253 : t381 = t208*t14
2422 4713253 : t383 = t77*r3*t5
2423 4713253 : t384 = t381*t383
2424 4713253 : t388 = t119*t218*t77
2425 : t391 = t370*t121 - t374*t377 - 0.2e1_dp/0.3e1_dp*t380*t384 - (2._dp &
2426 4713253 : *t118*t388)
2427 4713253 : t395 = rootpi
2428 4713253 : t396 = 0.1e1_dp/t395
2429 4713253 : t397 = t123*t396
2430 4713253 : t398 = f32**2
2431 4713253 : t399 = t398*t1
2432 4713253 : t400 = t399*t119
2433 4713253 : t401 = t15*t131
2434 4713253 : t403 = EXP(-t400*t401)
2435 4713253 : t404 = t126*t27
2436 4713253 : t405 = t125*t404
2437 4713253 : t406 = t129*sscale
2438 4713253 : t408 = t132*r3*t5
2439 4713253 : t416 = t125*t126*t128
2440 4713253 : t417 = 0.1e1_dp/t132
2441 4713253 : t418 = t246*t77
2442 4713253 : t419 = t34*t251
2443 4713253 : t420 = t77*t277
2444 4713253 : t422 = t418 - t419*t420
2445 4713253 : t423 = t417*t422
2446 4713253 : t424 = t406*t423
2447 : t428 = t403*(-t405*t406*t408/0.3e1_dp - t127*t128*t13*t133 &
2448 4713253 : + t416*t424/0.2e1_dp)
2449 4713253 : t429 = t397*t428
2450 : t433 = (t359*t113 - 0.7e1_dp/0.2e1_dp*t108*t365 - (t368*t391 &
2451 4713253 : *t123*t136) + (2._dp*t368*t429))*t140
2452 4713253 : t435 = t141*t142
2453 4713253 : t437 = t110*t92*t144
2454 4713253 : t438 = t437*t325
2455 4713253 : t441 = -t433*t145 - 0.7e1_dp/0.2e1_dp*t435*t438
2456 4713253 : t442 = t441*E
2457 4713253 : t443 = t442*t149
2458 4713253 : t444 = 0.1e1_dp/t361
2459 4713253 : t445 = t444*t325
2460 4713253 : t451 = t151*t147
2461 4713253 : t452 = t93*t159
2462 4713253 : t453 = t452*t325
2463 4713253 : t456 = 0.1e1_dp/t158/t157
2464 4713253 : t457 = t83*t456
2465 4713253 : t461 = -t215 - t221 + t248 - t280 - 0.2e1_dp/0.3e1_dp*t72*t44*t6
2466 4713253 : t462 = t156*t461
2467 4713253 : t463 = t457*t462
2468 4713253 : t466 = t93*t163
2469 4713253 : t470 = 0.1e1_dp/t162/t154
2470 4713253 : t471 = t83*t470
2471 4713253 : t478 = t166*t90
2472 4713253 : t479 = t93*t170
2473 4713253 : t480 = t479*t325
2474 4713253 : t483 = 0.1e1_dp/t169/t168
2475 4713253 : t484 = t83*t483
2476 4713253 : t485 = t155*t461
2477 4713253 : t486 = t484*t485
2478 : t490 = (-t151*t441*t83*t159 + t451*t453 + 0.5e1_dp/0.2e1_dp*t451 &
2479 : *t463 + t81*t466*t325 + t81*t471*t461/0.2e1_dp - t166 &
2480 : *t341*t83*t170 + t478*t480 + 0.3e1_dp/0.2e1_dp*t478*t486)* &
2481 4713253 : omega
2482 4713253 : t493 = t27*r3*t5
2483 4713253 : t499 = t176*t147
2484 4713253 : t500 = t149*t159
2485 4713253 : t501 = t500*t325
2486 4713253 : t504 = t93*t456
2487 4713253 : t505 = t504*t462
2488 4713253 : t511 = t180*t90
2489 4713253 : t512 = t149*t170
2490 4713253 : t513 = t512*t325
2491 4713253 : t516 = t93*t483
2492 4713253 : t517 = t516*t485
2493 : t521 = (-t176*t441*t93*t159 + (2._dp*t499*t501) + 0.5e1_dp/ &
2494 : 0.2e1_dp*(t499)*(t505) - t180*t341*t93*t170 + (2._dp &
2495 4713253 : *t511*t513) + 0.3e1_dp/0.2e1_dp*(t511)*(t517))*t185
2496 4713253 : t523 = t189*t13
2497 4713253 : t526 = t148*t444
2498 4713253 : t528 = t194*t44*t325
2499 4713253 : t531 = t149*t456
2500 4713253 : t532 = t148*t531
2501 4713253 : t533 = t193*t44
2502 4713253 : t534 = t533*t462
2503 4713253 : t537 = t148*t500
2504 4713253 : t538 = t193*t258
2505 4713253 : t539 = t538*t6
2506 4713253 : t542 = dexeirho(Q, dQrho)
2507 4713253 : t544 = t199**2
2508 4713253 : t545 = 0.1e1_dp/t544
2509 4713253 : t546 = t75*t545
2510 4713253 : t548 = t285*t200 - t546*t285
2511 4713253 : t549 = 0.1e1_dp/t75
2512 4713253 : t550 = t548*t549
2513 : t554 = -t81*t93*t325 + t342*t94 - (2._dp*t91*t345) + t443 &
2514 : - (3._dp*t148*t445) + t490*t128 - t174*t493/0.3e1_dp + t521 &
2515 : *t190 - t186*t523 - t443*t195 + (3._dp*t526*t528) + 0.5e1_dp &
2516 : /0.2e1_dp*t532*t534 + 0.5e1_dp/0.3e1_dp*t537*t539 + t197*(t542 &
2517 4713253 : + t550*t199)
2518 4713253 : t555 = t554*Clda
2519 4713253 : e_rho = e_rho + (-0.4e1_dp/0.3e1_dp*t322*t206 - t80*t555)*sx
2520 4713253 : t561 = F2*t300
2521 4713253 : t564 = t561*t68 - t85*t317
2522 4713253 : t565 = t15*t564
2523 4713253 : t567 = 2._dp*t288*t88 + t11*t565
2524 4713253 : t568 = f12*t567
2525 4713253 : t570 = t344*t320
2526 : t584 = t97*(t99*t567*t82 + t99*t90*t320 + 2._dp*t102*t82 &
2527 4713253 : *t320 + 3._dp*t104*t92*t320)
2528 4713253 : t586 = t364*t320
2529 4713253 : t589 = f94*t300
2530 4713253 : t590 = t589*t117
2531 4713253 : t592 = t375*t316
2532 4713253 : t593 = t18*t592
2533 4713253 : t595 = t68*ndrho
2534 4713253 : t596 = t116*t595
2535 4713253 : t599 = t590*t121 - t374*t593 + 2._dp*t596*t121
2536 4713253 : t603 = f32*t126
2537 4713253 : t604 = t603*t128
2538 4713253 : t607 = t301*t77
2539 4713253 : t608 = t77*t316
2540 4713253 : t610 = t607 - t419*t608
2541 4713253 : t611 = t417*t610
2542 4713253 : t612 = t406*t611
2543 4713253 : t616 = t403*(t604*t406*t132 + t416*t612/0.2e1_dp)
2544 4713253 : t617 = t397*t616
2545 : t621 = (t584*t113 - 0.7e1_dp/0.2e1_dp*t108*t586 - (t368*t599 &
2546 4713253 : *t123*t136) + (2._dp*t368*t617))*t140
2547 4713253 : t623 = t437*t320
2548 4713253 : t626 = -t621*t145 - 0.7e1_dp/0.2e1_dp*t435*t623
2549 4713253 : t627 = t626*E
2550 4713253 : t628 = t627*t149
2551 4713253 : t629 = t444*t320
2552 4713253 : t635 = t452*t320
2553 4713253 : t637 = t156*t320
2554 4713253 : t638 = t457*t637
2555 4713253 : t649 = t479*t320
2556 4713253 : t651 = t155*t320
2557 4713253 : t652 = t484*t651
2558 : t656 = (-t151*t626*t83*t159 + t451*t635 + 0.5e1_dp/0.2e1_dp*t451 &
2559 : *t638 + t81*t466*t320 + t81*t471*t320/0.2e1_dp - t166 &
2560 : *t567*t83*t170 + t478*t649 + 0.3e1_dp/0.2e1_dp*t478*t652)* &
2561 4713253 : omega
2562 4713253 : t661 = t500*t320
2563 4713253 : t664 = t504*t637
2564 4713253 : t670 = t512*t320
2565 4713253 : t673 = t516*t651
2566 : t677 = (-t176*t626*t93*t159 + (2._dp*t499*t661) + 0.5e1_dp/ &
2567 : 0.2e1_dp*(t499)*(t664) - t180*t567*t93*t170 + (2._dp &
2568 4713253 : *t511*t670) + 0.3e1_dp/0.2e1_dp*(t511)*(t673))*t185
2569 4713253 : t681 = t194*t44*t320
2570 4713253 : t684 = t533*t637
2571 4713253 : t687 = dexeindrho(Q, dQndrho)
2572 4713253 : t690 = t320*t200 - t546*t320
2573 4713253 : t691 = t690*t549
2574 : t695 = -t81*t93*t320 + t568*t94 - (2._dp*t91*t570) + t628 &
2575 : - (3._dp*t148*t629) + t656*t128 + t677*t190 - t628*t195 &
2576 : + (3._dp*t526*t681) + 0.5e1_dp/0.2e1_dp*t532*t684 + t197*(t687 &
2577 4713253 : + t691*t199)
2578 4713253 : t696 = t695*Clda
2579 4713253 : e_ndrho = e_ndrho + (-t80*t696)*sx
2580 : END IF
2581 :
2582 7448304 : END SUBROUTINE xwpbe_lda_calc_3
2583 :
2584 : ! **************************************************************************************************
2585 : !> \brief Evaluates the screened hole averaged PBE exchange functional for lda.
2586 : !> \param e_0 ...
2587 : !> \param e_rho ...
2588 : !> \param e_ndrho ...
2589 : !> \param rho , ndrho: density and norm of the density gradient
2590 : !> \param ndrho ...
2591 : !> \param omega screening parameter
2592 : !> \param sscale scaling factor to enforce Lieb-Oxford bound
2593 : !> \param sx scaling factor
2594 : !> \param order maximum degree of the derivative that should be evaluated
2595 : !> \par History
2596 : !> 05.2007 created [Manuel Guidon]
2597 : !> \author Manuel Guidon
2598 : !> \note
2599 : !> This routine evaluates the functional for omega!=0 using a simple
2600 : !> gaussian expansion for large ww and a taylor expansion for the
2601 : !> parameter G.
2602 : ! **************************************************************************************************
2603 122 : SUBROUTINE xwpbe_lda_calc_4(e_0, e_rho, e_ndrho, rho, ndrho, omega, sscale, sx, order)
2604 : REAL(KIND=dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho
2605 : REAL(KIND=dp), INTENT(IN) :: rho, ndrho, omega, sscale, sx
2606 : INTEGER, INTENT(IN) :: order
2607 :
2608 : REAL(KIND=dp) :: dQndrho, dQrho, Q, t1, t10, t100, t102, t103, t105, t106, t108, t109, t11, &
2609 : t110, t111, t113, t114, t115, t116, t117, t118, t12, t121, t122, t125, t126, t127, t128, &
2610 : t129, t13, t133, t134, t136, t137, t14, t140, t141, t145, t146, t149, t15, t151, t153, &
2611 : t154, t155, t156, t158, t159, t16, t160, t161, t163, t167, t169, t17, t170, t171, t173, &
2612 : t174, t176, t177, t178, t179, t18, t180, t182, t183, t184, t185, t188, t189, t19, t192, &
2613 : t193, t195, t196, t197, t198, t199, t2, t203, t206, t207, t208, t209, t21, t210, t211, &
2614 : t212, t213, t219, t22, t220, t221, t222, t223, t227, t229, t23
2615 : REAL(KIND=dp) :: t235, t238, t239, t24, t240, t241, t246, t248, t249, t25, t252, t253, t256, &
2616 : t257, t258, t261, t262, t263, t265, t266, t269, t27, t270, t274, t277, t278, t279, t28, &
2617 : t281, t283, t286, t289, t29, t290, t293, t296, t299, t3, t300, t302, t303, t305, t306, &
2618 : t309, t31, t310, t313, t316, t32, t321, t326, t327, t329, t330, t331, t332, t333, t334, &
2619 : t34, t340, t341, t342, t345, t346, t35, t350, t351, t352, t355, t359, t36, t360, t367, &
2620 : t368, t369, t372, t373, t374, t375, t379, t38, t382, t388, t389, t39, t390, t393, t394, &
2621 : t4, t400, t401, t402, t405, t406, t41, t410, t412, t415, t417
2622 : REAL(KIND=dp) :: t42, t420, t421, t422, t423, t426, t427, t428, t431, t433, t434, t435, &
2623 : t437, t438, t439, t44, t443, t444, t450, t453, t454, t456, t457, t459, t46, t464, t465, &
2624 : t468, t469, t472, t473, t475, t476, t477, t478, t48, t484, t486, t487, t49, t498, t5, &
2625 : t500, t501, t505, t510, t513, t519, t522, t526, t530, t533, t536, t539, t54, t540, t544, &
2626 : t545, t55, t56, t58, t6, t60, t61, t63, t65, t67, t68, t69, t7, t70, t71, t72, t73, t74, &
2627 : t75, t77, t78, t8, t80, t81, t82, t83, t84, t85, t87, t88, t9, t90, t91, t92, t93, t94, &
2628 : t95, t96, t97, t99
2629 :
2630 122 : t1 = ndrho**2
2631 122 : t2 = r2**2
2632 122 : t3 = 0.1e1_dp/t2
2633 122 : t4 = t1*t3
2634 122 : t5 = pi**2
2635 122 : t6 = r3*t5
2636 122 : t7 = t6*rho
2637 122 : t8 = t7**(0.1e1_dp/0.3e1_dp)
2638 122 : t9 = t8**2
2639 122 : t10 = 0.1e1_dp/t9
2640 122 : t11 = t4*t10
2641 122 : t12 = rho**2
2642 122 : t13 = 0.1e1_dp/t12
2643 122 : t14 = sscale**2
2644 122 : t15 = t13*t14
2645 122 : t16 = a1*t1
2646 122 : t17 = t16*t3
2647 122 : t18 = t10*t13
2648 122 : t19 = t18*t14
2649 122 : t21 = t1**2
2650 122 : t22 = a2*t21
2651 122 : t23 = t2**2
2652 122 : t24 = 0.1e1_dp/t23
2653 122 : t25 = t22*t24
2654 122 : t27 = 0.1e1_dp/t8/t7
2655 122 : t28 = t12**2
2656 122 : t29 = 0.1e1_dp/t28
2657 122 : t31 = t14**2
2658 122 : t32 = t27*t29*t31
2659 122 : t34 = t17*t19 + t25*t32
2660 122 : t35 = a3*t21
2661 122 : t36 = t35*t24
2662 122 : t38 = t21*ndrho
2663 122 : t39 = a4*t38
2664 122 : t41 = 0.1e1_dp/t23/r2
2665 122 : t42 = t39*t41
2666 122 : t44 = 0.1e1_dp/t9/t7
2667 122 : t46 = 0.1e1_dp/t28/rho
2668 122 : t48 = t31*sscale
2669 122 : t49 = t44*t46*t48
2670 122 : t54 = 0.1e1_dp/t23/t2
2671 122 : t55 = a5*t21*t1*t54
2672 122 : t56 = r3**2
2673 122 : t58 = t5**2
2674 122 : t60 = 0.1e1_dp/t56/t58
2675 122 : t61 = t28**2
2676 122 : t63 = t31*t14
2677 122 : t65 = t60/t61*t63
2678 122 : t67 = r1 + t36*t32 + t42*t49 + t55*t65
2679 122 : t68 = 0.1e1_dp/t67
2680 122 : t69 = t34*t68
2681 122 : t70 = t15*t69
2682 122 : t71 = t11*t70
2683 122 : t72 = omega**2
2684 122 : t73 = beta2*t72
2685 122 : t74 = t73*t10
2686 122 : t75 = t71 + t74
2687 122 : t77 = 0.1e1_dp/A
2688 122 : Q = f94*t75*t77
2689 122 : t78 = rho**(0.1e1_dp/0.3e1_dp)
2690 122 : t80 = t78*rho*f89
2691 122 : t81 = B*f12
2692 122 : t82 = t71 + DD
2693 122 : t83 = 0.1e1_dp/t82
2694 122 : t84 = t81*t83
2695 122 : t85 = F2*t34
2696 122 : t87 = F1 + t85*t68
2697 122 : t88 = t15*t87
2698 122 : t90 = t11*t88 + r1
2699 122 : t91 = f12*t90
2700 122 : t92 = t82**2
2701 122 : t93 = 0.1e1_dp/t92
2702 122 : t94 = C*t93
2703 122 : t95 = t91*t94
2704 122 : t96 = g2*t1
2705 122 : t97 = t96*t3
2706 122 : t99 = g3*t21
2707 122 : t100 = t99*t24
2708 122 : t102 = g1 + t97*t19 + t100*t32
2709 122 : t103 = t15*t102
2710 122 : t105 = t11*t103 + r1
2711 122 : t106 = t105*E
2712 122 : t108 = 0.1e1_dp/t92/t82
2713 122 : t109 = t106*t108
2714 122 : t110 = f158*E
2715 122 : t111 = t105*t83
2716 122 : t113 = t71 + DD + t72*t10
2717 122 : t114 = t113**2
2718 122 : t115 = t114**2
2719 122 : t116 = t115*t113
2720 122 : t117 = SQRT(t116)
2721 122 : t118 = 0.1e1_dp/t117
2722 122 : t121 = SQRT(t113)
2723 122 : t122 = 0.1e1_dp/t121
2724 122 : t125 = f68*C
2725 122 : t126 = t90*t83
2726 122 : t127 = t114*t113
2727 122 : t128 = SQRT(t127)
2728 122 : t129 = 0.1e1_dp/t128
2729 : t133 = (-t110*t111*t118 - t81*t83*t122 - t125*t126*t129) &
2730 122 : *omega
2731 122 : t134 = 0.1e1_dp/t8
2732 122 : t136 = f52*E
2733 122 : t137 = t105*t93
2734 122 : t140 = f12*C
2735 122 : t141 = t90*t93
2736 122 : t145 = t72*omega
2737 122 : t146 = (-t136*t137*t118 - t140*t141*t129)*t145
2738 122 : t149 = 0.1e1_dp/r3/t5
2739 122 : t151 = t149/rho
2740 122 : t153 = t72**2
2741 122 : t154 = t153*omega
2742 122 : t155 = t118*t154
2743 122 : t156 = t155*t44
2744 122 : t158 = f12*A
2745 122 : t159 = exei(Q)
2746 122 : t160 = t71 + DD + t74
2747 122 : t161 = 0.1e1_dp/t160
2748 122 : t163 = LOG(t75*t161)
2749 : t167 = (t84 + t95 + t109 + t133*t134 + t146*t151 - t109*t156 &
2750 122 : + t158*(t159 + t163))*Clda
2751 122 : e_0 = e_0 + (-t80*t167)*sx
2752 122 : IF (order >= 1) THEN
2753 71 : t169 = t44*t13
2754 71 : t170 = t4*t169
2755 71 : t171 = t14*t34
2756 71 : t173 = t68*r3*t5
2757 71 : t174 = t171*t173
2758 71 : t176 = 0.2e1_dp/0.3e1_dp*t170*t174
2759 71 : t177 = t12*rho
2760 71 : t178 = 0.1e1_dp/t177
2761 71 : t179 = t178*t14
2762 71 : t180 = t179*t69
2763 71 : t182 = 2._dp*t11*t180
2764 71 : t183 = t3*t44
2765 71 : t184 = t16*t183
2766 71 : t185 = t15*t6
2767 71 : t188 = t10*t178
2768 71 : t189 = t188*t14
2769 71 : t192 = t56*t58
2770 71 : t193 = t192*t12
2771 71 : t195 = 0.1e1_dp/t8/t193
2772 71 : t196 = t24*t195
2773 71 : t197 = t22*t196
2774 71 : t198 = t29*t31
2775 71 : t199 = t198*t6
2776 71 : t203 = t27*t46*t31
2777 : t206 = -0.2e1_dp/0.3e1_dp*t184*t185 - (2._dp*t17*t189) - 0.4e1_dp &
2778 71 : /0.3e1_dp*t197*t199 - (4._dp*t25*t203)
2779 71 : t207 = t206*t68
2780 71 : t208 = t15*t207
2781 71 : t209 = t11*t208
2782 71 : t210 = t4*t18
2783 71 : t211 = t67**2
2784 71 : t212 = 0.1e1_dp/t211
2785 71 : t213 = t35*t196
2786 71 : t219 = 0.1e1_dp/t9/t193
2787 71 : t220 = t41*t219
2788 71 : t221 = t39*t220
2789 71 : t222 = t46*t48
2790 71 : t223 = t222*t6
2791 71 : t227 = 0.1e1_dp/t28/t12
2792 71 : t229 = t44*t227*t48
2793 71 : t235 = t60/t61/rho*t63
2794 : t238 = -0.4e1_dp/0.3e1_dp*t213*t199 - (4._dp*t36*t203) - 0.5e1_dp &
2795 71 : /0.3e1_dp*t221*t223 - (5._dp*t42*t229) - (8._dp*t55*t235)
2796 71 : t239 = t212*t238
2797 71 : t240 = t171*t239
2798 71 : t241 = t210*t240
2799 : t246 = -t176 - t182 + t209 - t241 - 0.2e1_dp/0.3e1_dp*t73*t44*r3 &
2800 71 : *t5
2801 71 : dQrho = f94*t246*t77
2802 71 : t248 = ndrho*t3
2803 71 : t249 = t248*t10
2804 71 : t252 = a1*ndrho
2805 71 : t253 = t252*t3
2806 71 : t256 = t1*ndrho
2807 71 : t257 = a2*t256
2808 71 : t258 = t257*t24
2809 71 : t261 = 2._dp*t253*t19 + 4._dp*t258*t32
2810 71 : t262 = t261*t68
2811 71 : t263 = t15*t262
2812 71 : t265 = a3*t256
2813 71 : t266 = t265*t24
2814 71 : t269 = a4*t21
2815 71 : t270 = t269*t41
2816 71 : t274 = a5*t38*t54
2817 71 : t277 = 4._dp*t266*t32 + 5._dp*t270*t49 + 6._dp*t274*t65
2818 71 : t278 = t212*t277
2819 71 : t279 = t171*t278
2820 71 : t281 = 2._dp*t249*t70 + t11*t263 - t210*t279
2821 71 : dQndrho = f94*t281*t77
2822 71 : t283 = t78*f89
2823 71 : t286 = -t176 - t182 + t209 - t241
2824 71 : t289 = t14*t87
2825 71 : t290 = t289*t6
2826 71 : t293 = t179*t87
2827 71 : t296 = F2*t206
2828 71 : t299 = t296*t68 - t85*t239
2829 71 : t300 = t15*t299
2830 : t302 = -0.2e1_dp/0.3e1_dp*t170*t290 - (2._dp*t11*t293) + (t11 &
2831 71 : *t300)
2832 71 : t303 = f12*t302
2833 71 : t305 = C*t108
2834 71 : t306 = t305*t286
2835 71 : t309 = t14*t102
2836 71 : t310 = t309*t6
2837 71 : t313 = t179*t102
2838 71 : t316 = t96*t183
2839 71 : t321 = t99*t196
2840 : t326 = -0.2e1_dp/0.3e1_dp*t316*t185 - (2._dp*t97*t189) - 0.4e1_dp &
2841 71 : /0.3e1_dp*t321*t199 - (4._dp*t100*t203)
2842 71 : t327 = t15*t326
2843 : t329 = -0.2e1_dp/0.3e1_dp*t170*t310 - (2._dp*t11*t313) + (t11 &
2844 71 : *t327)
2845 71 : t330 = t329*E
2846 71 : t331 = t330*t108
2847 71 : t332 = t92**2
2848 71 : t333 = 0.1e1_dp/t332
2849 71 : t334 = t333*t286
2850 71 : t340 = t110*t105
2851 71 : t341 = t93*t118
2852 71 : t342 = t341*t286
2853 71 : t345 = 0.1e1_dp/t117/t116
2854 71 : t346 = t83*t345
2855 71 : t350 = -t176 - t182 + t209 - t241 - 0.2e1_dp/0.3e1_dp*t72*t44*t6
2856 71 : t351 = t115*t350
2857 71 : t352 = t346*t351
2858 71 : t355 = t93*t122
2859 71 : t359 = 0.1e1_dp/t121/t113
2860 71 : t360 = t83*t359
2861 71 : t367 = t125*t90
2862 71 : t368 = t93*t129
2863 71 : t369 = t368*t286
2864 71 : t372 = 0.1e1_dp/t128/t127
2865 71 : t373 = t83*t372
2866 71 : t374 = t114*t350
2867 71 : t375 = t373*t374
2868 : t379 = (-t110*t329*t83*t118 + t340*t342 + 0.5e1_dp/0.2e1_dp*t340 &
2869 : *t352 + t81*t355*t286 + t81*t360*t350/0.2e1_dp - t125 &
2870 : *t302*t83*t129 + t367*t369 + 0.3e1_dp/0.2e1_dp*t367*t375)* &
2871 71 : omega
2872 71 : t382 = t27*r3*t5
2873 71 : t388 = t136*t105
2874 71 : t389 = t108*t118
2875 71 : t390 = t389*t286
2876 71 : t393 = t93*t345
2877 71 : t394 = t393*t351
2878 71 : t400 = t140*t90
2879 71 : t401 = t108*t129
2880 71 : t402 = t401*t286
2881 71 : t405 = t93*t372
2882 71 : t406 = t405*t374
2883 : t410 = (-t136*t329*t93*t118 + (2._dp*t388*t390) + 0.5e1_dp/ &
2884 : 0.2e1_dp*(t388)*(t394) - t140*t302*t93*t129 + (2._dp &
2885 71 : *t400*t402) + 0.3e1_dp/0.2e1_dp*(t400)*(t406))*t145
2886 71 : t412 = t149*t13
2887 71 : t415 = t106*t333
2888 71 : t417 = t155*t44*t286
2889 71 : t420 = t108*t345
2890 71 : t421 = t106*t420
2891 71 : t422 = t154*t44
2892 71 : t423 = t422*t351
2893 71 : t426 = t106*t389
2894 71 : t427 = t154*t219
2895 71 : t428 = t427*t6
2896 71 : t431 = dexeirho(Q, dQrho)
2897 71 : t433 = t160**2
2898 71 : t434 = 0.1e1_dp/t433
2899 71 : t435 = t75*t434
2900 71 : t437 = t246*t161 - t435*t246
2901 71 : t438 = 0.1e1_dp/t75
2902 71 : t439 = t437*t438
2903 : t443 = -t81*t93*t286 + t303*t94 - (2._dp*t91*t306) + t331 &
2904 : - (3._dp*t106*t334) + t379*t134 - t133*t382/0.3e1_dp + t410 &
2905 : *t151 - t146*t412 - t331*t156 + (3._dp*t415*t417) + 0.5e1_dp &
2906 : /0.2e1_dp*t421*t423 + 0.5e1_dp/0.3e1_dp*t426*t428 + t158*(t431 &
2907 71 : + t439*t160)
2908 71 : t444 = t443*Clda
2909 71 : e_rho = e_rho + (-0.4e1_dp/0.3e1_dp*t283*t167 - t80*t444)*sx
2910 71 : t450 = F2*t261
2911 71 : t453 = t450*t68 - t85*t278
2912 71 : t454 = t15*t453
2913 71 : t456 = 2._dp*t249*t88 + t11*t454
2914 71 : t457 = f12*t456
2915 71 : t459 = t305*t281
2916 71 : t464 = g2*ndrho
2917 71 : t465 = t464*t3
2918 71 : t468 = g3*t256
2919 71 : t469 = t468*t24
2920 71 : t472 = 2._dp*t465*t19 + 4._dp*t469*t32
2921 71 : t473 = t15*t472
2922 71 : t475 = 2._dp*t249*t103 + t11*t473
2923 71 : t476 = t475*E
2924 71 : t477 = t476*t108
2925 71 : t478 = t333*t281
2926 71 : t484 = t341*t281
2927 71 : t486 = t115*t281
2928 71 : t487 = t346*t486
2929 71 : t498 = t368*t281
2930 71 : t500 = t114*t281
2931 71 : t501 = t373*t500
2932 : t505 = (-t110*t475*t83*t118 + t340*t484 + 0.5e1_dp/0.2e1_dp*t340 &
2933 : *t487 + t81*t355*t281 + t81*t360*t281/0.2e1_dp - t125 &
2934 : *t456*t83*t129 + t367*t498 + 0.3e1_dp/0.2e1_dp*t367*t501)* &
2935 71 : omega
2936 71 : t510 = t389*t281
2937 71 : t513 = t393*t486
2938 71 : t519 = t401*t281
2939 71 : t522 = t405*t500
2940 : t526 = (-t136*t475*t93*t118 + (2._dp*t388*t510) + 0.5e1_dp/ &
2941 : 0.2e1_dp*(t388)*(t513) - t140*t456*t93*t129 + (2._dp &
2942 71 : *t400*t519) + 0.3e1_dp/0.2e1_dp*(t400)*(t522))*t145
2943 71 : t530 = t155*t44*t281
2944 71 : t533 = t422*t486
2945 71 : t536 = dexeindrho(Q, dQndrho)
2946 71 : t539 = t281*t161 - t435*t281
2947 71 : t540 = t539*t438
2948 : t544 = -t81*t93*t281 + t457*t94 - (2._dp*t91*t459) + t477 &
2949 : - (3._dp*t106*t478) + t505*t134 + t526*t151 - t477*t156 &
2950 : + (3._dp*t415*t530) + 0.5e1_dp/0.2e1_dp*t421*t533 + t158*(t536 &
2951 71 : + t540*t160)
2952 71 : t545 = t544*Clda
2953 71 : e_ndrho = e_ndrho + (-t80*t545)*sx
2954 : END IF
2955 :
2956 122 : END SUBROUTINE xwpbe_lda_calc_4
2957 :
2958 : ! **************************************************************************************************
2959 : !> \brief return various information on the functional
2960 : !> \param reference string with the reference of the actual functional
2961 : !> \param shortform string with the shortform of the functional name
2962 : !> \param needs the components needed by this functional are set to
2963 : !> true (does not set the unneeded components to false)
2964 : !> \param max_deriv ...
2965 : !> \par History
2966 : !> 05.2007 created [Manuel Guidon]
2967 : !> \author Manuel Guidon
2968 : ! **************************************************************************************************
2969 1164 : SUBROUTINE xwpbe_lsd_info(reference, shortform, needs, max_deriv)
2970 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
2971 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
2972 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
2973 :
2974 1164 : IF (PRESENT(reference)) THEN
2975 2 : reference = "Jochen Heyd and Gustavo E. Scuseria, J. Chem. Phys., 120, 7274 {LSD version}"
2976 : END IF
2977 1164 : IF (PRESENT(shortform)) THEN
2978 2 : shortform = "shortrange part of PBE exchange {LSD}"
2979 : END IF
2980 1164 : IF (PRESENT(needs)) THEN
2981 1162 : needs%rho_spin = .TRUE.
2982 1162 : needs%norm_drho_spin = .TRUE.
2983 : END IF
2984 1164 : IF (PRESENT(max_deriv)) max_deriv = 1
2985 1164 : END SUBROUTINE xwpbe_lsd_info
2986 :
2987 : ! **************************************************************************************************
2988 : !> \brief evaluates the screened hole averaged PBE exchange functional for lsd
2989 : !> \param rho_set the density where you want to evaluate the functional
2990 : !> \param deriv_set place where to store the functional derivatives (they are
2991 : !> added to the derivatives)
2992 : !> \param order maximum degree of the derivative that should be evaluated
2993 : !> \param xwpbe_params input parameters (scaling,omega)
2994 : !> \par History
2995 : !> 05.2007 created [Manuel Guidon]
2996 : !> \author Manuel Guidon
2997 : !> \note
2998 : !> The current version provides code for derivatives up to second order.
2999 : !> Using the maple sheet in cp2k/doc it is straightforward to produce routines
3000 : !> for higher derivatives.
3001 : ! **************************************************************************************************
3002 6620 : SUBROUTINE xwpbe_lsd_eval(rho_set, deriv_set, order, xwpbe_params)
3003 :
3004 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
3005 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
3006 : INTEGER, INTENT(IN) :: order
3007 : TYPE(section_vals_type), POINTER :: xwpbe_params
3008 :
3009 : CHARACTER(len=*), PARAMETER :: routineN = 'xwpbe_lsd_eval'
3010 :
3011 : INTEGER :: handle, npoints
3012 : INTEGER, DIMENSION(2, 3) :: bo
3013 : REAL(kind=dp) :: epsilon_norm_drho, epsilon_rho, omega, &
3014 : sx, sx0
3015 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
3016 1324 : POINTER :: dummy, e_0, e_ndrhoa, e_ndrhob, e_rhoa, &
3017 1324 : e_rhob, norm_drhoa, norm_drhob, rhoa, &
3018 1324 : rhob
3019 : TYPE(xc_derivative_type), POINTER :: deriv
3020 :
3021 1324 : CALL timeset(routineN, handle)
3022 :
3023 1324 : CALL cite_reference(Heyd2004)
3024 :
3025 1324 : CALL section_vals_val_get(xwpbe_params, "SCALE_X", r_val=sx)
3026 1324 : CALL section_vals_val_get(xwpbe_params, "SCALE_X0", r_val=sx0)
3027 1324 : CALL section_vals_val_get(xwpbe_params, "OMEGA", r_val=omega)
3028 :
3029 : CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
3030 : norm_drhob=norm_drhob, local_bounds=bo, rho_cutoff=epsilon_rho, &
3031 1324 : drho_cutoff=epsilon_norm_drho)
3032 1324 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
3033 :
3034 1324 : dummy => rhoa
3035 :
3036 1324 : e_0 => dummy
3037 1324 : e_rhoa => dummy
3038 1324 : e_rhob => dummy
3039 1324 : e_ndrhoa => dummy
3040 1324 : e_ndrhob => dummy
3041 :
3042 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
3043 1324 : allocate_deriv=.TRUE.)
3044 1324 : CALL xc_derivative_get(deriv, deriv_data=e_0)
3045 1324 : IF (order >= 1 .OR. order == -1) THEN
3046 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
3047 1144 : allocate_deriv=.TRUE.)
3048 1144 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
3049 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
3050 1144 : allocate_deriv=.TRUE.)
3051 1144 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
3052 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
3053 1144 : allocate_deriv=.TRUE.)
3054 1144 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
3055 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
3056 1144 : allocate_deriv=.TRUE.)
3057 1144 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
3058 : END IF
3059 1324 : IF (order > 1) THEN
3060 0 : CPABORT("derivatives bigger than 1 not implemented")
3061 : END IF
3062 :
3063 : !$OMP PARALLEL DEFAULT(NONE) &
3064 : !$OMP SHARED(npoints, order, rhoa, norm_drhoa, e_0, e_rhoa, e_ndrhoa) &
3065 : !$OMP SHARED(epsilon_rho, sx, sx0, omega) &
3066 1324 : !$OMP SHARED(rhob, norm_drhob, e_rhob, e_ndrhob)
3067 :
3068 : !Call lsd_calc for alpha - and beta - spins
3069 :
3070 : CALL xwpbe_lsd_calc(npoints, order, rho=rhoa, norm_drho=norm_drhoa, &
3071 : e_0=e_0, e_rho=e_rhoa, e_ndrho=e_ndrhoa, &
3072 : epsilon_rho=epsilon_rho, &
3073 : sx=sx, sx0=sx0, omega=omega)
3074 : CALL xwpbe_lsd_calc(npoints, order, rho=rhob, norm_drho=norm_drhob, &
3075 : e_0=e_0, e_rho=e_rhob, e_ndrho=e_ndrhob, &
3076 : epsilon_rho=epsilon_rho, &
3077 : sx=sx, sx0=sx0, omega=omega)
3078 :
3079 : !$OMP END PARALLEL
3080 :
3081 1324 : CALL timestop(handle)
3082 :
3083 1324 : END SUBROUTINE xwpbe_lsd_eval
3084 :
3085 : ! **************************************************************************************************
3086 : !> \brief evaluates the screened hole averaged PBE exchange functional for lsd.
3087 : !> \param npoints ...
3088 : !> \param order degree of the derivative that should be evaluated,
3089 : !> if positive all the derivatives up to the given degree are evaluated,
3090 : !> if negative only the given degree is calculated
3091 : !> \param rho , ndrho: density and norm of the density gradient
3092 : !> \param norm_drho ...
3093 : !> \param e_0 ...
3094 : !> \param e_rho ...
3095 : !> \param e_ndrho ...
3096 : !> \param epsilon_rho ...
3097 : !> \param sx , sx0: scaling factor for omega!=0 and omega=0
3098 : !> \param sx0 ...
3099 : !> \param omega screening parameter
3100 : !> \par History
3101 : !> 05.2007 created [Manuel Guidon]
3102 : !> \author Manuel Guidon
3103 : !> \note
3104 : !> - The lsd part is calculated using the spin-scaling relations for the
3105 : !> exchange energy:
3106 : !>
3107 : !> Ex[na,nb] = 0.5 * Ex[2*na] + 0.5 * Ex[2*nb].
3108 : !>
3109 : !> - In order to avoid numerical instabilities, this routine calls different
3110 : !> subroutines. There are 4 routines for the case omega!=0 and 2 routines
3111 : !> for omega=0.
3112 : ! **************************************************************************************************
3113 2648 : SUBROUTINE xwpbe_lsd_calc(npoints, order, rho, norm_drho, e_0, e_rho, e_ndrho, &
3114 : epsilon_rho, sx, sx0, omega)
3115 :
3116 : INTEGER, INTENT(in) :: npoints, order
3117 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, norm_drho, e_0, e_rho, e_ndrho
3118 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, sx0, omega
3119 :
3120 : INTEGER :: ip
3121 : REAL(dp) :: e_0_temp, gradient_before, &
3122 : gradient_increment, my_ndrho, my_rho, &
3123 : pullback
3124 : REAL(KIND=dp) :: ss, ss2, sscale, t1, t2, t3, t4, t5, t6, &
3125 : t7, t8, ww
3126 :
3127 : !$OMP DO
3128 :
3129 : DO ip = 1, npoints
3130 : !According to spin-scaling relation, we need twice the density and its gradient
3131 34895932 : my_rho = 2.0_dp*MAX(rho(ip), 0.0_dp)
3132 34895932 : IF (my_rho > epsilon_rho) THEN
3133 34522233 : gradient_before = e_ndrho(ip)
3134 34522233 : my_ndrho = 2.0_dp*MAX(norm_drho(ip), 0.0_dp)
3135 :
3136 : !Do some precalculation in order to catch the correct branch afterwards
3137 34522233 : sscale = 1.0_dp
3138 34522233 : t1 = pi**2
3139 34522233 : t2 = t1*my_rho
3140 34522233 : t3 = t2**(0.1e1_dp/0.3e1_dp)
3141 34522233 : t4 = 0.1e1_dp/t3
3142 34522233 : t5 = omega*t4
3143 34522233 : ww = 0.6933612743506347048433524e0_dp*t5
3144 34522233 : t6 = my_ndrho*t4
3145 34522233 : t7 = 0.1e1_dp/my_rho
3146 34522233 : t8 = t7*sscale
3147 34522233 : ss = 0.3466806371753173524216762e0_dp*t6*t8
3148 34522233 : IF (ss > scutoff) THEN
3149 17547329 : ss2 = ss*ss
3150 17547329 : sscale = ((smax)*ss2 - (sconst))/(ss2*ss)
3151 : END IF
3152 34522233 : e_0_temp = 0.0_dp
3153 34522233 : IF (sx0 /= 0.0_dp) THEN
3154 : !original PBE hole
3155 28282855 : IF (ss*sscale > gcutoff) THEN
3156 : CALL xwpbe_lda_calc_0(e_0_temp, e_rho(ip), e_ndrho(ip), my_rho, &
3157 28255119 : my_ndrho, sscale, sx0, order)
3158 : ELSE
3159 : CALL xwpbe_lda_calc_01(e_0_temp, e_rho(ip), e_ndrho(ip), my_rho, &
3160 27736 : my_ndrho, sscale, sx0, order)
3161 : END IF
3162 : !According to spin-scaling relation, we need only half of the energy
3163 28282855 : e_0(ip) = e_0(ip) + 0.5_dp*e_0_temp
3164 : END IF
3165 34522233 : e_0_temp = 0.0_dp
3166 34522233 : IF (sx /= 0.0_dp) THEN
3167 7628414 : IF (omega == 0.0_dp) THEN
3168 116 : IF (ss*sscale > gcutoff) THEN
3169 : CALL xwpbe_lda_calc_0(e_0_temp, e_rho(ip), e_ndrho(ip), my_rho, &
3170 88 : my_ndrho, sscale, sx, order)
3171 : ELSE
3172 : CALL xwpbe_lda_calc_01(e_0_temp, e_rho(ip), e_ndrho(ip), my_rho, &
3173 28 : my_ndrho, sscale, sx, order)
3174 : END IF
3175 7628298 : ELSE IF (ww < wcutoff .AND. ss*sscale > gcutoff) THEN
3176 : CALL xwpbe_lda_calc_1(e_0_temp, e_rho(ip), e_ndrho(ip), my_rho, &
3177 6607485 : my_ndrho, omega, sscale, sx, order)
3178 6357 : ELSE IF (ww < wcutoff .AND. ss*sscale <= gcutoff) THEN
3179 : CALL xwpbe_lda_calc_2(e_0_temp, e_rho(ip), e_ndrho(ip), my_rho, &
3180 6357 : my_ndrho, omega, sscale, sx, order)
3181 1014456 : ELSE IF (ww >= wcutoff .AND. ss*sscale > gcutoff) THEN
3182 : CALL xwpbe_lda_calc_3(e_0_temp, e_rho(ip), e_ndrho(ip), my_rho, &
3183 1014456 : my_ndrho, omega, sscale, sx, order)
3184 : ELSE
3185 : CALL xwpbe_lda_calc_4(e_0_temp, e_rho(ip), e_ndrho(ip), my_rho, &
3186 0 : my_ndrho, omega, sscale, sx, order)
3187 : END IF
3188 : !According to spin-scaling relation, we need only half of the energy
3189 : END IF
3190 34522233 : e_0(ip) = e_0(ip) + 0.5_dp*e_0_temp
3191 34522233 : IF (order >= 1 .AND. ss > scutoff) THEN
3192 : ! Pull back g_eff = g*f(s), s = g/(2*k_F*rho), through the large-s scaling.
3193 : ! The generated expressions differentiate at fixed f(s).
3194 11865498 : pullback = (-smax*ss2 + 3.0_dp*sconst)/(smax*ss2 - sconst)
3195 11865498 : gradient_increment = e_ndrho(ip) - gradient_before
3196 11865498 : e_rho(ip) = e_rho(ip) - gradient_increment*(4.0_dp*my_ndrho/(3.0_dp*my_rho))*pullback
3197 11865498 : e_ndrho(ip) = gradient_before + gradient_increment*(1.0_dp + pullback)
3198 : END IF
3199 : END IF
3200 : END DO
3201 :
3202 : !$OMP END DO
3203 :
3204 2648 : END SUBROUTINE xwpbe_lsd_calc
3205 :
3206 : ! **************************************************************************************************
3207 : !> \brief These functions evaluate products exp(x)*Ei(x) and pi*exp(x)*erfc(sqrt(x)),
3208 : !> as well as their derivatives with respect to various combinations of
3209 : !> rho and norm_drho.
3210 : !> \param Q , dQrho, dQndrho, d2Qrhondrho :
3211 : !> Argument Q and derivatives with respect to various combinations of
3212 : !> rho and norm_drho
3213 : !> \param exp_q optional precomputed EXP(Q)
3214 : !> \return ...
3215 : !> \par History
3216 : !> 05.2007 created [Manuel Guidon]
3217 : !> \author Manuel Guidon
3218 : !> \note
3219 : !> - In order to avoid numerical instabilities, these routines use Taylor-
3220 : !> expansions for the above core-products for large arguments.
3221 : !> - When adapting this module for higher order derivatives, appropriate
3222 : !> functions have to be provided!
3223 : ! **************************************************************************************************
3224 160920083 : FUNCTION exei(Q, exp_q)
3225 : REAL(dp), INTENT(IN) :: Q
3226 : REAL(dp), INTENT(IN), OPTIONAL :: exp_q
3227 : REAL(dp) :: exei
3228 :
3229 160920083 : exei = 0.0_dp
3230 160920083 : IF (Q < expcutoff) THEN
3231 : !Use exact product
3232 144045009 : IF (PRESENT(exp_q)) THEN
3233 87518835 : exei = exp_q*expint(1, Q, 1.0_dp/exp_q)
3234 : ELSE
3235 56526174 : exei = EXP(Q)*expint(1, Q)
3236 : END IF
3237 : ELSE
3238 : !Use approximation
3239 16875074 : exei = (1._dp/Q)*(Q*Q + exei1*Q + exei2)/(Q*Q + exei3*Q + exei4)
3240 : END IF
3241 160920083 : END FUNCTION exei
3242 :
3243 : ! **************************************************************************************************
3244 : !> \brief ...
3245 : !> \param Q ...
3246 : !> \return ...
3247 : ! **************************************************************************************************
3248 56219081 : FUNCTION exer(Q)
3249 : REAL(dp), INTENT(IN) :: Q
3250 : REAL(dp) :: exer
3251 :
3252 : REAL(dp) :: Q3, Q5
3253 :
3254 56219081 : exer = 0.0_dp
3255 56219081 : IF (Q < expcutoff) THEN
3256 : !Use exact expression
3257 56219081 : exer = pi*EXP(Q)*erfc(SQRT(Q))
3258 : ELSE
3259 : !Use approximation
3260 0 : Q3 = Q*Q*Q
3261 0 : Q5 = Q3*Q*Q
3262 0 : exer = pi*(1.0_dp/SQRT(Q*pi) - 1.0_dp/(2.0_dp*SQRT(pi*Q3)) + 3.0_dp/(4.0_dp*(SQRT(pi*Q5))))
3263 : END IF
3264 56219081 : END FUNCTION exer
3265 :
3266 : ! **************************************************************************************************
3267 : !> \brief ...
3268 : !> \param Q ...
3269 : !> \param dQrho ...
3270 : !> \param exei_q optional precomputed exei(Q)
3271 : !> \return ...
3272 : ! **************************************************************************************************
3273 85670823 : FUNCTION dexeirho(Q, dQrho, exei_q)
3274 : REAL(dp), INTENT(IN) :: Q, dQrho
3275 : REAL(dp), INTENT(IN), OPTIONAL :: exei_q
3276 : REAL(dp) :: dexeirho
3277 :
3278 85670823 : IF (PRESENT(exei_q)) THEN
3279 65194786 : dexeirho = dQrho*(exei_q - 1.0_dp/Q)
3280 : ELSE
3281 20476037 : dexeirho = dQrho*(exei(Q) - 1.0_dp/Q)
3282 : END IF
3283 85670823 : END FUNCTION dexeirho
3284 :
3285 : ! **************************************************************************************************
3286 : !> \brief ...
3287 : !> \param Q ...
3288 : !> \param dQndrho ...
3289 : !> \param exei_q optional precomputed exei(Q)
3290 : !> \return ...
3291 : ! **************************************************************************************************
3292 85670823 : FUNCTION dexeindrho(Q, dQndrho, exei_q)
3293 : REAL(dp), INTENT(IN) :: Q, dQndrho
3294 : REAL(dp), INTENT(IN), OPTIONAL :: exei_q
3295 : REAL(dp) :: dexeindrho
3296 :
3297 85670823 : IF (PRESENT(exei_q)) THEN
3298 65194786 : dexeindrho = dQndrho*(exei_q - 1.0_dp/Q)
3299 : ELSE
3300 20476037 : dexeindrho = dQndrho*(exei(Q) - 1.0_dp/Q)
3301 : END IF
3302 85670823 : END FUNCTION dexeindrho
3303 :
3304 : ! **************************************************************************************************
3305 : !> \brief ...
3306 : !> \param Q ...
3307 : !> \param dQrho ...
3308 : !> \return ...
3309 : ! **************************************************************************************************
3310 15762713 : FUNCTION dexerrho(Q, dQrho)
3311 : REAL(dp), INTENT(IN) :: Q, dQrho
3312 : REAL(dp) :: dexerrho
3313 :
3314 15762713 : dexerrho = dQrho*exer(Q) - dQrho*rootpi/SQRT(Q)
3315 15762713 : END FUNCTION dexerrho
3316 :
3317 : ! **************************************************************************************************
3318 : !> \brief ...
3319 : !> \param Q ...
3320 : !> \param dQndrho ...
3321 : !> \return ...
3322 : ! **************************************************************************************************
3323 15762713 : FUNCTION dexerndrho(Q, dQndrho)
3324 : REAL(dp), INTENT(IN) :: Q, dQndrho
3325 : REAL(dp) :: dexerndrho
3326 :
3327 15762713 : dexerndrho = dQndrho*exer(Q) - dQndrho*rootpi/SQRT(Q)
3328 15762713 : END FUNCTION dexerndrho
3329 :
3330 : ! **************************************************************************************************
3331 : !> \brief ...
3332 : !> \param Q ...
3333 : !> \param dQrho ...
3334 : !> \param d2Qrhorho ...
3335 : !> \return ...
3336 : ! **************************************************************************************************
3337 0 : FUNCTION d2exeirhorho(Q, dQrho, d2Qrhorho)
3338 : REAL(dp), INTENT(IN) :: Q, dQrho, d2Qrhorho
3339 : REAL(dp) :: d2exeirhorho
3340 :
3341 : d2exeirhorho = exei(Q)*(d2Qrhorho + dQrho*dQrho) + &
3342 0 : 1.0_dp/(Q*Q)*(-Q*dQrho*dQrho - Q*d2Qrhorho + dQrho*dQrho)
3343 0 : END FUNCTION d2exeirhorho
3344 :
3345 : ! **************************************************************************************************
3346 : !> \brief ...
3347 : !> \param Q ...
3348 : !> \param dQrho ...
3349 : !> \param d2Qrhorho ...
3350 : !> \return ...
3351 : ! **************************************************************************************************
3352 0 : FUNCTION d2exerrhorho(Q, dQrho, d2Qrhorho)
3353 : REAL(dp), INTENT(IN) :: Q, dQrho, d2Qrhorho
3354 : REAL(dp) :: d2exerrhorho
3355 :
3356 : REAL(dp) :: pi12, Q12
3357 :
3358 0 : Q12 = SQRT(Q)
3359 0 : pi12 = rootpi
3360 :
3361 : d2exerrhorho = exer(Q)*(d2Qrhorho + dQrho*dQrho) - dQrho*dQrho/(pi12*Q12) + &
3362 0 : 0.5_dp*dQrho*dQrho/(pi12*Q*Q12) - d2Qrhorho/(pi12*Q12)
3363 0 : END FUNCTION d2exerrhorho
3364 :
3365 : ! **************************************************************************************************
3366 : !> \brief ...
3367 : !> \param Q ...
3368 : !> \param dQrho ...
3369 : !> \param dQndrho ...
3370 : !> \param d2Qrhondrho ...
3371 : !> \return ...
3372 : ! **************************************************************************************************
3373 0 : FUNCTION d2exeirhondrho(Q, dQrho, dQndrho, d2Qrhondrho)
3374 : REAL(dp), INTENT(IN) :: Q, dQrho, dQndrho, d2Qrhondrho
3375 : REAL(dp) :: d2exeirhondrho
3376 :
3377 : d2exeirhondrho = exei(Q)*(d2Qrhondrho + dQrho*dQndrho) - &
3378 0 : 1.0_dp/Q*(dQrho*dQndrho + d2Qrhondrho) + 1.0_dp/(Q*Q)*dQrho*dQndrho
3379 0 : END FUNCTION d2exeirhondrho
3380 :
3381 : ! **************************************************************************************************
3382 : !> \brief ...
3383 : !> \param Q ...
3384 : !> \param dQrho ...
3385 : !> \param dQndrho ...
3386 : !> \param d2Qrhondrho ...
3387 : !> \return ...
3388 : ! **************************************************************************************************
3389 0 : FUNCTION d2exerrhondrho(Q, dQrho, dQndrho, d2Qrhondrho)
3390 : REAL(dp), INTENT(IN) :: Q, dQrho, dQndrho, d2Qrhondrho
3391 : REAL(dp) :: d2exerrhondrho
3392 :
3393 : REAL(dp) :: pi12, Q12
3394 :
3395 0 : Q12 = SQRT(Q)
3396 0 : pi12 = rootpi
3397 :
3398 : d2exerrhondrho = exer(Q)*(d2Qrhondrho + dQrho*dQndrho) - 1.0_dp/(pi12*Q12)*dQrho*dQndrho &
3399 0 : + 0.5_dp/(pi12*Q12*Q)*dQrho*dQndrho - 1.0_dp/(pi12*Q12)*d2Qrhondrho
3400 0 : END FUNCTION d2exerrhondrho
3401 :
3402 : ! **************************************************************************************************
3403 : !> \brief ...
3404 : !> \param Q ...
3405 : !> \param dQndrho ...
3406 : !> \param d2Qndrhondrho ...
3407 : !> \return ...
3408 : ! **************************************************************************************************
3409 0 : FUNCTION d2exeindrhondrho(Q, dQndrho, d2Qndrhondrho)
3410 : REAL(dp), INTENT(IN) :: Q, dQndrho, d2Qndrhondrho
3411 : REAL(dp) :: d2exeindrhondrho
3412 :
3413 : d2exeindrhondrho = exei(Q)*(d2Qndrhondrho + dQndrho*dQndrho) + &
3414 0 : 1.0_dp/(Q*Q)*(-Q*dQndrho*dQndrho - Q*d2Qndrhondrho + dQndrho*dQndrho)
3415 0 : END FUNCTION d2exeindrhondrho
3416 :
3417 : ! **************************************************************************************************
3418 : !> \brief ...
3419 : !> \param Q ...
3420 : !> \param dQndrho ...
3421 : !> \param d2Qndrhondrho ...
3422 : !> \return ...
3423 : ! **************************************************************************************************
3424 0 : FUNCTION d2exerndrhondrho(Q, dQndrho, d2Qndrhondrho)
3425 : REAL(dp), INTENT(IN) :: Q, dQndrho, d2Qndrhondrho
3426 : REAL(dp) :: d2exerndrhondrho
3427 :
3428 : REAL(dp) :: pi12, Q12
3429 :
3430 0 : Q12 = SQRT(Q)
3431 0 : pi12 = rootpi
3432 :
3433 : d2exerndrhondrho = exer(Q)*(d2Qndrhondrho + dQndrho*dQndrho) - dQndrho*dQndrho/(pi12*Q12) &
3434 0 : + 0.5_dp*dQndrho*dQndrho/(pi12*Q*Q12) - d2Qndrhondrho/(pi12*Q12)
3435 0 : END FUNCTION d2exerndrhondrho
3436 :
3437 3297089165 : @:define_expint(impure=0)
3438 :
3439 : END MODULE xc_xwpbe
3440 :
|
