Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief calculates the longrange part of Becke 88 exchange functional
10 : !> \par History
11 : !> 01.2008 created [mguidon]
12 : !> \author Manuel Guidon
13 : ! **************************************************************************************************
14 : MODULE xc_xbecke88_long_range
15 : USE bibliography, ONLY: Becke1988,&
16 : cite_reference
17 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
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 :
38 : IMPLICIT NONE
39 : PRIVATE
40 :
41 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbecke88_long_range'
43 : REAL(kind=dp), PARAMETER :: beta = 0.0042_dp
44 :
45 : PUBLIC :: xb88_lr_lda_info, xb88_lr_lsd_info, xb88_lr_lda_eval, xb88_lr_lsd_eval
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief return various information on the functional
50 : !> \param reference string with the reference of the actual functional
51 : !> \param shortform string with the shortform of the functional name
52 : !> \param needs the components needed by this functional are set to
53 : !> true (does not set the unneeded components to false)
54 : !> \param max_deriv ...
55 : !> \par History
56 : !> 01.2008 created [mguidon]
57 : !> \author Manuel Guidon
58 : ! **************************************************************************************************
59 1298 : SUBROUTINE xb88_lr_lda_info(reference, shortform, needs, max_deriv)
60 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
61 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
62 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
63 :
64 1298 : IF (PRESENT(reference)) THEN
65 8 : reference = "A. Becke, Phys. Rev. A 38, 3098 (1988) {LDA version}"
66 : END IF
67 1298 : IF (PRESENT(shortform)) THEN
68 8 : shortform = "Becke 1988 Exchange Functional (LDA)"
69 : END IF
70 1298 : IF (PRESENT(needs)) THEN
71 1290 : needs%rho = .TRUE.
72 1290 : needs%norm_drho = .TRUE.
73 : END IF
74 1298 : IF (PRESENT(max_deriv)) max_deriv = 3
75 :
76 1298 : END SUBROUTINE xb88_lr_lda_info
77 :
78 : ! **************************************************************************************************
79 : !> \brief return various information on the functional
80 : !> \param reference string with the reference of the actual functional
81 : !> \param shortform string with the shortform of the functional name
82 : !> \param needs the components needed by this functional are set to
83 : !> true (does not set the unneeded components to false)
84 : !> \param max_deriv ...
85 : !> \par History
86 : !> 01.2008 created [mguidon]
87 : !> \author Manuel Guidon
88 : ! **************************************************************************************************
89 57 : SUBROUTINE xb88_lr_lsd_info(reference, shortform, needs, max_deriv)
90 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
91 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
92 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
93 :
94 57 : IF (PRESENT(reference)) THEN
95 1 : reference = "A. Becke, Phys. Rev. A 38, 3098 (1988) {LSD version}"
96 : END IF
97 57 : IF (PRESENT(shortform)) THEN
98 1 : shortform = "Becke 1988 Exchange Functional (LSD)"
99 : END IF
100 57 : IF (PRESENT(needs)) THEN
101 56 : needs%rho_spin = .TRUE.
102 56 : needs%rho_spin_1_3 = .TRUE.
103 56 : needs%norm_drho_spin = .TRUE.
104 : END IF
105 57 : IF (PRESENT(max_deriv)) max_deriv = 3
106 :
107 57 : END SUBROUTINE xb88_lr_lsd_info
108 :
109 : ! **************************************************************************************************
110 : !> \brief evaluates the becke 88 longrange exchange functional for lda
111 : !> \param rho_set the density where you want to evaluate the functional
112 : !> \param deriv_set place where to store the functional derivatives (they are
113 : !> added to the derivatives)
114 : !> \param grad_deriv degree of the derivative that should be evaluated,
115 : !> if positive all the derivatives up to the given degree are evaluated,
116 : !> if negative only the given degree is calculated
117 : !> \param xb88_lr_params input parameters (scaling)
118 : !> \par History
119 : !> 01.2008 created [mguidon]
120 : !> \author Manuel Guidon
121 : ! **************************************************************************************************
122 4896 : SUBROUTINE xb88_lr_lda_eval(rho_set, deriv_set, grad_deriv, xb88_lr_params)
123 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
124 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
125 : INTEGER, INTENT(in) :: grad_deriv
126 : TYPE(section_vals_type), POINTER :: xb88_lr_params
127 :
128 : CHARACTER(len=*), PARAMETER :: routineN = 'xb88_lr_lda_eval'
129 :
130 : INTEGER :: handle, npoints
131 : INTEGER, DIMENSION(2, 3) :: bo
132 : REAL(kind=dp) :: epsilon_rho, omega, sx
133 1224 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
134 1224 : e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
135 1224 : e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho
136 : TYPE(xc_derivative_type), POINTER :: deriv
137 :
138 1224 : CALL timeset(routineN, handle)
139 :
140 1224 : CALL section_vals_val_get(xb88_lr_params, "scale_x", r_val=sx)
141 1224 : CALL section_vals_val_get(xb88_lr_params, "omega", r_val=omega)
142 :
143 1224 : CALL cite_reference(Becke1988)
144 :
145 : CALL xc_rho_set_get(rho_set, rho=rho, &
146 1224 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
147 1224 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
148 :
149 : ! meaningful default for the arrays we don't need: let us make compiler
150 : ! and debugger happy...
151 1224 : dummy => rho
152 :
153 1224 : e_0 => dummy
154 1224 : e_rho => dummy
155 1224 : e_ndrho => dummy
156 1224 : e_rho_rho => dummy
157 1224 : e_ndrho_rho => dummy
158 1224 : e_ndrho_ndrho => dummy
159 1224 : e_rho_rho_rho => dummy
160 1224 : e_ndrho_rho_rho => dummy
161 1224 : e_ndrho_ndrho_rho => dummy
162 1224 : e_ndrho_ndrho_ndrho => dummy
163 :
164 1224 : IF (grad_deriv >= 0) THEN
165 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
166 1224 : allocate_deriv=.TRUE.)
167 1224 : CALL xc_derivative_get(deriv, deriv_data=e_0)
168 : END IF
169 1224 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
170 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
171 1216 : allocate_deriv=.TRUE.)
172 1216 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
173 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
174 1216 : allocate_deriv=.TRUE.)
175 1216 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
176 : END IF
177 1224 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
178 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
179 378 : allocate_deriv=.TRUE.)
180 378 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
181 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
182 378 : allocate_deriv=.TRUE.)
183 378 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
184 : deriv => xc_dset_get_derivative(deriv_set, &
185 378 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
186 378 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
187 : END IF
188 1224 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
189 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
190 0 : allocate_deriv=.TRUE.)
191 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
192 : deriv => xc_dset_get_derivative(deriv_set, &
193 0 : [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.TRUE.)
194 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
195 : deriv => xc_dset_get_derivative(deriv_set, &
196 0 : [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.TRUE.)
197 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
198 : deriv => xc_dset_get_derivative(deriv_set, &
199 0 : [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
200 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
201 : END IF
202 1224 : IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
203 0 : CPABORT("derivatives bigger than 3 not implemented")
204 : END IF
205 :
206 : !$OMP PARALLEL DEFAULT(NONE) &
207 : !$OMP SHARED(rho, norm_drho, e_0, e_rho, e_ndrho, e_rho_rho) &
208 : !$OMP SHARED(e_ndrho_rho, e_ndrho_ndrho, e_rho_rho_rho) &
209 : !$OMP SHARED(e_ndrho_rho_rho, e_ndrho_ndrho_rho) &
210 : !$OMP SHARED(e_ndrho_ndrho_ndrho, grad_deriv, npoints) &
211 1224 : !$OMP SHARED(epsilon_rho, sx, omega)
212 :
213 : CALL xb88_lr_lda_calc(rho=rho, norm_drho=norm_drho, &
214 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
215 : e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
216 : e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
217 : e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, &
218 : e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, grad_deriv=grad_deriv, &
219 : npoints=npoints, epsilon_rho=epsilon_rho, &
220 : sx=sx, omega=omega)
221 :
222 : !$OMP END PARALLEL
223 :
224 1224 : CALL timestop(handle)
225 1224 : END SUBROUTINE xb88_lr_lda_eval
226 :
227 : ! **************************************************************************************************
228 : !> \brief low level calculation of the becke 88 exchange functional for lda
229 : !> \param rho alpha or beta spin density
230 : !> \param norm_drho || grad rho ||
231 : !> \param e_0 adds to it the local value of the functional
232 : !> \param e_rho derivative of the functional wrt. to the variables
233 : !> named where the * is.
234 : !> \param e_ndrho ...
235 : !> \param e_rho_rho ...
236 : !> \param e_ndrho_rho ...
237 : !> \param e_ndrho_ndrho ...
238 : !> \param e_rho_rho_rho ...
239 : !> \param e_ndrho_rho_rho ...
240 : !> \param e_ndrho_ndrho_rho ...
241 : !> \param e_ndrho_ndrho_ndrho ...
242 : !> \param grad_deriv ...
243 : !> \param npoints ...
244 : !> \param epsilon_rho ...
245 : !> \param sx scaling-parameter for exchange
246 : !> \param omega parameter for erfc
247 : !> \par History
248 : !> 01.2008 created [mguidon]
249 : !> \author Manuel Guidon
250 : !> \note
251 : !> - Just took the lsd code and scaled rho and ndrho by 1/2 (e_0 with 2.0)
252 : !> - Derivatives higher than 1 not tested
253 : ! **************************************************************************************************
254 1224 : SUBROUTINE xb88_lr_lda_calc(rho, norm_drho, &
255 1224 : e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
256 1224 : e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
257 1224 : e_ndrho_ndrho_ndrho, grad_deriv, npoints, epsilon_rho, sx, &
258 : omega)
259 : INTEGER, INTENT(in) :: npoints, grad_deriv
260 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_ndrho, &
261 : e_ndrho_ndrho_rho, e_ndrho_rho_rho, e_rho_rho_rho, e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, &
262 : e_ndrho, e_rho, e_0
263 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho
264 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, omega
265 :
266 : INTEGER :: ii
267 : REAL(kind=dp) :: Cx, epsilon_rho43, my_epsilon_rho, my_ndrho, my_rho, t1, t10, t100, t1002, &
268 : t1009, t101, t1013, t103, t104, t1044, t105, t1069, t109, t1091, t11, t1102, t112, t113, &
269 : t1136, t1141, t1146, t1153, t1156, t116, t1163, t1167, t117, t1177, t118, t12, t122, &
270 : t1220, t126, t127, t128, t1284, t130, t132, t133, t1334, t1341, t1345, t135, t136, t137, &
271 : t1370, t1396, t140, t1400, t1404, t1405, t141, t1412, t143, t1449, t1456, t146, t1467, &
272 : t147, t1472, t148, t151, t155, t1553, t16, t168, t169, t17, t172, t173, t176, t177, t18, &
273 : t185, t186, t190, t196, t2, t200, t207, t21, t211, t212, t213
274 : REAL(kind=dp) :: t216, t219, t22, t221, t222, t225, t226, t23, t230, t231, t232, t233, t237, &
275 : t24, t241, t244, t245, t250, t251, t254, t255, t258, t259, t26, t264, t265, t27, t270, &
276 : t271, t28, t281, t284, t285, t289, t29, t293, t297, t3, t30, t304, t31, t311, t312, t313, &
277 : t316, t319, t321, t323, t325, t326, t328, t330, t335, t339, t34, t343, t346, t347, t351, &
278 : t354, t358, t36, t365, t368, t37, t372, t377, t38, t382, t383, t384, t39, t393, t397, &
279 : t40, t400, t401, t404, t405, t408, t41, t412, t413, t417, t418, t42, t428, t429, t43, &
280 : t435, t44, t451, t452, t455, t457, t46, t460, t462, t463, t464
281 : REAL(kind=dp) :: t465, t466, t467, t47, t470, t473, t478, t479, t48, t482, t486, t489, t49, &
282 : t492, t496, t5, t501, t502, t505, t508, t51, t513, t514, t519, t52, t521, t525, t526, &
283 : t529, t530, t533, t534, t536, t537, t539, t55, t562, t566, t570, t571, t572, t573, t574, &
284 : t575, t576, t580, t586, t59, t590, t6, t60, t601, t605, t606, t613, t624, t628, t632, &
285 : t639, t64, t640, t641, t657, t667, t669, t677, t68, t687, t69, t7, t71, t716, t72, t722, &
286 : t735, t739, t746, t75, t76, t769, t77, t79, t791, t793, t8, t838, t84, t842, t846, t85, &
287 : t857, t86, t860, t867, t87, t880, t90, t91, t933, t94, t95, t961
288 : REAL(kind=dp) :: t98, t99, xx
289 :
290 1224 : my_epsilon_rho = epsilon_rho
291 1224 : epsilon_rho43 = my_epsilon_rho**(4.0_dp/3.0_dp)
292 1224 : Cx = 1.5_dp*(3.0_dp/4.0_dp/pi)**(1.0_dp/3.0_dp)
293 :
294 : !$OMP DO
295 : DO ii = 1, npoints
296 : !! scale densities with 0.5 because code comes from LSD
297 22177632 : my_rho = rho(ii)*0.5_dp
298 22177632 : my_ndrho = norm_drho(ii)*0.5_dp
299 22177632 : IF (my_rho > my_epsilon_rho) THEN
300 21158517 : IF (grad_deriv >= 0) THEN
301 21158517 : t1 = my_rho**(0.1e1_dp/0.3e1_dp)
302 21158517 : t2 = t1*my_rho
303 21158517 : t3 = 0.1e1_dp/t2
304 21158517 : xx = my_ndrho*MAX(t3, epsilon_rho43)
305 21158517 : t5 = my_ndrho**2
306 21158517 : t6 = beta*t5
307 21158517 : t7 = my_rho**2
308 21158517 : t8 = t1**2
309 21158517 : t10 = 0.1e1_dp/t8/t7
310 21158517 : t11 = beta*my_ndrho
311 21158517 : t12 = LOG(xx + SQRT(xx**0.2e1_dp + 0.1e1_dp))
312 21158517 : t16 = 0.10e1_dp + 0.60e1_dp*t11*t3*t12
313 21158517 : t17 = 0.1e1_dp/t16
314 21158517 : t18 = t10*t17
315 21158517 : t21 = 0.20e1_dp*Cx + 0.20e1_dp*t6*t18
316 21158517 : t22 = SQRT(t21)
317 21158517 : t23 = t22*t21
318 21158517 : t24 = my_rho*t23
319 21158517 : t26 = rootpi
320 21158517 : t27 = 0.1e1_dp/t26
321 21158517 : t28 = 0.1e1_dp/omega
322 21158517 : t29 = 0.1e1_dp/t22
323 21158517 : t30 = t28*t29
324 21158517 : t31 = t26*t1
325 21158517 : t34 = erf(0.3000000000e1_dp*t30*t31)
326 21158517 : t36 = omega*t22
327 21158517 : t37 = 0.1e1_dp/t1
328 21158517 : t38 = t27*t37
329 21158517 : t39 = omega**2
330 21158517 : t40 = 0.1e1_dp/t39
331 21158517 : t41 = 0.1e1_dp/t21
332 21158517 : t42 = t40*t41
333 21158517 : t43 = pi*t8
334 21158517 : t44 = t42*t43
335 21158517 : t46 = EXP(-0.8999999998e1_dp*t44)
336 21158517 : t47 = t39*t21
337 21158517 : t48 = 0.1e1_dp/pi
338 21158517 : t49 = 0.1e1_dp/t8
339 21158517 : t51 = t46 - 0.10e1_dp
340 21158517 : t52 = t48*t49*t51
341 21158517 : t55 = t46 - 0.15e1_dp - 0.5555555558e-1_dp*t47*t52
342 21158517 : t59 = t26*t34 + 0.3333333334e0_dp*t36*t38*t55
343 21158517 : t60 = t27*t59
344 :
345 : !! Multiply with 2.0 because it code comes from LSD
346 21158517 : e_0(ii) = e_0(ii) - 0.2222222224e0_dp*t24*omega*t60*sx*2.0_dp
347 :
348 : END IF
349 21158517 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
350 20913217 : t64 = t23*omega
351 20913217 : t68 = my_rho*t22*omega
352 20913217 : t69 = t7*my_rho
353 20913217 : t71 = 0.1e1_dp/t8/t69
354 20913217 : t72 = t71*t17
355 20913217 : t75 = t16**2
356 20913217 : t76 = 0.1e1_dp/t75
357 20913217 : t77 = t10*t76
358 20913217 : t79 = 0.1e1_dp/t1/t7
359 20913217 : t84 = 1 + t5*t10
360 20913217 : t85 = SQRT(t84)
361 20913217 : t86 = 0.1e1_dp/t85
362 20913217 : t87 = t71*t86
363 20913217 : t90 = -0.8000000000e1_dp*t11*t79*t12 - 0.8000000000e1_dp*t6*t87
364 20913217 : t91 = t77*t90
365 20913217 : t94 = -0.5333333333e1_dp*t6*t72 - 0.20e1_dp*t6*t91
366 20913217 : t95 = t60*t94
367 20913217 : t98 = omega*t27
368 20913217 : t99 = SQRT(0.3141592654e1_dp)
369 20913217 : t100 = 0.1e1_dp/t99
370 20913217 : t101 = t26*t100
371 20913217 : t103 = EXP(-0.9000000000e1_dp*t44)
372 20913217 : t104 = 0.1e1_dp/t23
373 20913217 : t105 = t28*t104
374 20913217 : t109 = t26*t49
375 : t112 = -0.1500000000e1_dp*t105*t31*t94 + 0.1000000000e1_dp*t30* &
376 20913217 : t109
377 20913217 : t113 = t103*t112
378 20913217 : t116 = omega*t29
379 20913217 : t117 = t116*t27
380 20913217 : t118 = t37*t55
381 20913217 : t122 = t27*t3
382 20913217 : t126 = t21**2
383 20913217 : t127 = 0.1e1_dp/t126
384 20913217 : t128 = t40*t127
385 20913217 : t130 = t128*t43*t94
386 20913217 : t132 = pi*t37
387 20913217 : t133 = t42*t132
388 20913217 : t135 = 0.8999999998e1_dp*t130 - 0.5999999999e1_dp*t133
389 20913217 : t136 = t135*t46
390 20913217 : t137 = t39*t94
391 20913217 : t140 = t8*my_rho
392 20913217 : t141 = 0.1e1_dp/t140
393 20913217 : t143 = t48*t141*t51
394 20913217 : t146 = t47*t48
395 20913217 : t147 = t49*t135
396 20913217 : t148 = t147*t46
397 : t151 = t136 - 0.5555555558e-1_dp*t137*t52 + 0.3703703705e-1_dp*t47 &
398 20913217 : *t143 - 0.5555555558e-1_dp*t146*t148
399 : t155 = REAL(2*t101*t113, KIND=dp) + 0.1666666667e0_dp*t117*t118*t94 - &
400 : 0.1111111111e0_dp*t36*t122*t55 + 0.3333333334e0_dp*t36*t38* &
401 20913217 : t151
402 :
403 : e_rho(ii) = e_rho(ii) + (-0.2222222224e0_dp*t64*t60 - 0.3333333336e0_dp*t68*t95 - &
404 20913217 : 0.2222222224e0_dp*t24*t98*t155)*sx
405 :
406 20913217 : t168 = 0.60e1_dp*beta*t3*t12 + 0.60e1_dp*t11*t10*t86
407 20913217 : t169 = t77*t168
408 20913217 : t172 = 0.40e1_dp*t11*t18 - 0.20e1_dp*t6*t169
409 20913217 : t173 = t60*t172
410 20913217 : t176 = pi*t100
411 20913217 : t177 = t176*t103
412 20913217 : t185 = t128*pi
413 20913217 : t186 = t8*t172
414 20913217 : t190 = t39*t172
415 : t196 = 0.8999999998e1_dp*t185*t186*t46 - 0.5555555558e-1_dp*t190 &
416 20913217 : *t52 - 0.5000000001e0_dp*t41*t172*t46
417 : t200 = -0.3000000000e1_dp*t177*t105*t1*t172 + 0.1666666667e0_dp* &
418 20913217 : t117*t118*t172 + 0.3333333334e0_dp*t36*t38*t196
419 :
420 : e_ndrho(ii) = e_ndrho(ii) + (-0.3333333336e0_dp*t68*t173 - 0.2222222224e0_dp*t24*t98 &
421 20913217 : *t200)*sx
422 :
423 : END IF
424 21158517 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
425 3723264 : t207 = t27*t155
426 3723264 : t211 = my_rho*t29*omega
427 3723264 : t212 = t94**2
428 3723264 : t213 = t60*t212
429 3723264 : t216 = t207*t94
430 3723264 : t219 = t7**2
431 3723264 : t221 = 0.1e1_dp/t8/t219
432 3723264 : t222 = t221*t17
433 3723264 : t225 = t71*t76
434 3723264 : t226 = t225*t90
435 3723264 : t230 = 0.1e1_dp/t75/t16
436 3723264 : t231 = t10*t230
437 3723264 : t232 = t90**2
438 3723264 : t233 = t231*t232
439 3723264 : t237 = 0.1e1_dp/t1/t69
440 3723264 : t241 = t221*t86
441 3723264 : t244 = t5**2
442 3723264 : t245 = beta*t244
443 3723264 : t250 = 0.1e1_dp/t85/t84
444 3723264 : t251 = 0.1e1_dp/t1/t219/t69*t250
445 : t254 = 0.1866666667e2_dp*t11*t237*t12 + 0.4000000000e2_dp*t6*t241 &
446 3723264 : - 0.1066666667e2_dp*t245*t251
447 3723264 : t255 = t77*t254
448 : t258 = 0.1955555555e2_dp*t6*t222 + 0.1066666667e2_dp*t6*t226 + 0.40e1_dp &
449 3723264 : *t6*t233 - 0.20e1_dp*t6*t255
450 3723264 : t259 = t60*t258
451 3723264 : t264 = 0.9000000000e1_dp*t130 - 0.6000000000e1_dp*t133
452 3723264 : t265 = t264*t103
453 3723264 : t270 = 0.1e1_dp/t22/t126
454 3723264 : t271 = t28*t270
455 3723264 : t281 = t26*t141
456 : t284 = 0.2250000000e1_dp*t271*t31*t212 - 0.1000000000e1_dp*t105* &
457 : t109*t94 - 0.1500000000e1_dp*t105*t31*t258 - 0.6666666667e0_dp &
458 3723264 : *t30*t281
459 3723264 : t285 = t103*t284
460 3723264 : t289 = omega*t104*t27
461 3723264 : t293 = t3*t55
462 3723264 : t297 = t37*t151
463 3723264 : t304 = t27*t79
464 3723264 : t311 = t126*t21
465 3723264 : t312 = 0.1e1_dp/t311
466 3723264 : t313 = t40*t312
467 3723264 : t316 = 0.1800000000e2_dp*t313*t43*t212
468 3723264 : t319 = 0.1200000000e2_dp*t128*t132*t94
469 3723264 : t321 = t128*t43*t258
470 3723264 : t323 = pi*t3
471 3723264 : t325 = 0.2000000000e1_dp*t42*t323
472 3723264 : t326 = -t316 + t319 + 0.8999999998e1_dp*t321 + t325
473 3723264 : t328 = t135**2
474 3723264 : t330 = t39*t258
475 3723264 : t335 = t137*t48
476 3723264 : t339 = t48*t10*t51
477 3723264 : t343 = t141*t135*t46
478 3723264 : t346 = t49*t326
479 3723264 : t347 = t346*t46
480 3723264 : t351 = t49*t328*t46
481 : t354 = t326*t46 + t328*t46 - 0.5555555558e-1_dp*t330*t52 + 0.7407407410e-1_dp &
482 : *t137*t143 - 0.1111111112e0_dp*t335*t148 - 0.6172839508e-1_dp &
483 : *t47*t339 + 0.7407407410e-1_dp*t146*t343 - 0.5555555558e-1_dp &
484 3723264 : *t146*t347 - 0.5555555558e-1_dp*t146*t351
485 : t358 = REAL(2*t101*t265*t112, KIND=dp) + REAL(2*t101*t285, KIND=dp) - 0.8333333335e-1_dp &
486 : *t289*t118*t212 - 0.1111111111e0_dp*t117*t293*t94 &
487 : + 0.3333333334e0_dp*t117*t297*t94 + 0.1666666667e0_dp*t117* &
488 : t118*t258 + 0.1481481481e0_dp*t36*t304*t55 - 0.2222222222e0_dp* &
489 3723264 : t36*t122*t151 + 0.3333333334e0_dp*t36*t38*t354
490 :
491 : e_rho_rho(ii) = e_rho_rho(ii) + (-0.6666666672e0_dp*t36*t95 - 0.4444444448e0_dp*t64*t207 - &
492 : 0.1666666668e0_dp*t211*t213 - 0.6666666672e0_dp*t68*t216 - 0.3333333336e0_dp &
493 3723264 : *t68*t259 - 0.2222222224e0_dp*t24*t98*t358)*sx
494 :
495 3723264 : t365 = t27*t200
496 3723264 : t368 = t94*t172
497 3723264 : t372 = t365*t94
498 3723264 : t377 = t225*t168
499 3723264 : t382 = t6*t10
500 3723264 : t383 = t230*t90
501 3723264 : t384 = t383*t168
502 3723264 : t393 = beta*t5*my_ndrho
503 3723264 : t397 = 0.1e1_dp/t1/t219/t7*t250
504 : t400 = -0.8000000000e1_dp*beta*t79*t12 - 0.2400000000e2_dp*t11* &
505 3723264 : t87 + 0.8000000000e1_dp*t393*t397
506 3723264 : t401 = t77*t400
507 : t404 = -0.1066666667e2_dp*t11*t72 + 0.5333333333e1_dp*t6*t377 - 0.40e1_dp &
508 3723264 : *t11*t91 + 0.40e1_dp*t382*t384 - 0.20e1_dp*t6*t401
509 3723264 : t405 = t60*t404
510 3723264 : t408 = t207*t172
511 3723264 : t412 = t26*pi*t100
512 3723264 : t413 = t412*t128
513 3723264 : t417 = t271*t26
514 3723264 : t418 = t1*t94
515 : t428 = 0.2250000000e1_dp*t417*t418*t172 - 0.1500000000e1_dp*t105 &
516 3723264 : *t31*t404 - 0.5000000000e0_dp*t105*t109*t172
517 3723264 : t429 = t103*t428
518 3723264 : t435 = t37*t196
519 3723264 : t451 = t313*pi
520 3723264 : t452 = t8*t94
521 3723264 : t455 = 0.1800000000e2_dp*t451*t452*t172
522 3723264 : t457 = t128*t43*t404
523 3723264 : t460 = t128*t132*t172
524 3723264 : t462 = -t455 + 0.8999999998e1_dp*t457 + 0.5999999999e1_dp*t460
525 3723264 : t463 = t462*t46
526 3723264 : t464 = t135*t40
527 3723264 : t465 = t464*t127
528 3723264 : t466 = t172*t46
529 3723264 : t467 = t43*t466
530 3723264 : t470 = t39*t404
531 3723264 : t473 = t94*t127
532 3723264 : t478 = 0.1e1_dp/my_rho
533 3723264 : t479 = t41*t478
534 3723264 : t482 = t190*t48
535 3723264 : t486 = t49*t462*t46
536 3723264 : t489 = t41*t135
537 : t492 = t463 + 0.8999999998e1_dp*t465*t467 - 0.5555555558e-1_dp*t470 &
538 : *t52 - 0.5000000001e0_dp*t473*t466 + 0.3703703705e-1_dp*t190*t143 &
539 : + 0.3333333334e0_dp*t479*t466 - 0.5555555558e-1_dp*t482*t148 &
540 3723264 : - 0.5555555558e-1_dp*t146*t486 - 0.5000000001e0_dp*t489*t466
541 : t496 = 0.1800000000e2_dp*t413*t186*t113 + REAL(2*t101*t429, KIND=dp) &
542 : - 0.8333333335e-1_dp*t289*t118*t368 + 0.1666666667e0_dp*t117*t435 &
543 : *t94 + 0.1666666667e0_dp*t117*t118*t404 - 0.5555555555e-1_dp &
544 : *t117*t293*t172 - 0.1111111111e0_dp*t36*t122*t196 + 0.1666666667e0_dp &
545 3723264 : *t117*t297*t172 + 0.3333333334e0_dp*t36*t38*t492
546 :
547 : e_ndrho_rho(ii) = e_ndrho_rho(ii) + (-0.3333333336e0_dp*t36*t173 - 0.2222222224e0_dp*t64*t365 - &
548 : 0.1666666668e0_dp*t211*t60*t368 - 0.3333333336e0_dp*t68*t372 &
549 : - 0.3333333336e0_dp*t68*t405 - 0.3333333336e0_dp*t68*t408 - 0.2222222224e0_dp &
550 3723264 : *t24*t98*t496)*sx
551 :
552 3723264 : t501 = t172**2
553 3723264 : t502 = t60*t501
554 3723264 : t505 = t365*t172
555 3723264 : t508 = beta*t10
556 3723264 : t513 = t168**2
557 3723264 : t514 = t231*t513
558 3723264 : t519 = t219*my_rho
559 3723264 : t521 = 0.1e1_dp/t1/t519
560 3723264 : t525 = 0.120e2_dp*t508*t86 - 0.60e1_dp*t6*t521*t250
561 3723264 : t526 = t77*t525
562 : t529 = 0.40e1_dp*t508*t17 - 0.80e1_dp*t11*t169 + 0.40e1_dp*t6*t514 &
563 3723264 : - 0.20e1_dp*t6*t526
564 3723264 : t530 = t60*t529
565 3723264 : t533 = pi**2
566 3723264 : t534 = t533*t100
567 3723264 : t536 = 0.1e1_dp/t39/omega
568 3723264 : t537 = t534*t536
569 3723264 : t539 = 0.1e1_dp/t22/t311
570 3723264 : t562 = t8*t501
571 3723264 : t566 = t8*t529
572 3723264 : t570 = t39**2
573 3723264 : t571 = 0.1e1_dp/t570
574 3723264 : t572 = t126**2
575 3723264 : t573 = 0.1e1_dp/t572
576 3723264 : t574 = t571*t573
577 3723264 : t575 = t574*t533
578 3723264 : t576 = t2*t501
579 3723264 : t580 = t39*t529
580 : t586 = -0.2250000000e2_dp*t451*t562*t46 + 0.8999999998e1_dp*t185 &
581 : *t566*t46 + 0.8099999996e2_dp*t575*t576*t46 - 0.5555555558e-1_dp &
582 3723264 : *t580*t52 - 0.5000000001e0_dp*t41*t529*t46
583 : t590 = -0.2700000000e2_dp*t537*t539*my_rho*t501*t103 + 0.4500000000e1_dp &
584 : *t177*t271*t1*t501 - 0.3000000000e1_dp*t177*t105* &
585 : t1*t529 - 0.8333333335e-1_dp*t289*t118*t501 + 0.3333333334e0_dp &
586 : *t117*t435*t172 + 0.1666666667e0_dp*t117*t118*t529 + 0.3333333334e0_dp &
587 3723264 : *t36*t38*t586
588 :
589 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + (-0.1666666668e0_dp*t211*t502 - 0.6666666672e0_dp*t68*t505 &
590 : - 0.3333333336e0_dp*t68*t530 - 0.2222222224e0_dp*t24*t98*t590) &
591 3723264 : *sx
592 :
593 : END IF
594 21158517 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
595 0 : t601 = t27*t358
596 0 : t605 = my_rho*t104*omega
597 0 : t606 = t212*t94
598 0 : t613 = t94*t258
599 0 : t624 = 0.1e1_dp/t8/t519
600 0 : t628 = t221*t76
601 0 : t632 = t71*t230
602 0 : t639 = t75**2
603 0 : t640 = 0.1e1_dp/t639
604 0 : t641 = t10*t640
605 0 : t657 = t219**2
606 0 : t667 = t84**2
607 0 : t669 = 0.1e1_dp/t85/t667
608 : t677 = -0.9125925923e2_dp*t6*t624*t17 - 0.5866666667e2_dp*t6*t628 &
609 : *t90 - 0.3200000001e2_dp*t6*t632*t232 + 0.1600000000e2_dp*t6 &
610 : *t225*t254 - 0.120e2_dp*t6*t641*t232*t90 + 0.120e2_dp*t382 &
611 : *t383*t254 - 0.20e1_dp*t6*t77*(-0.6222222223e2_dp*t11/t1/ &
612 : t219*t12 - 0.2115555556e3_dp*t6*t624*t86 + 0.1315555556e3_dp* &
613 : t245/t1/t657*t250 - 0.4266666668e2_dp*beta*t244*t5/t657 &
614 0 : /t69*t669)
615 0 : t687 = t28*t539
616 0 : t716 = t264**2
617 0 : t722 = omega*t270*t27
618 0 : t735 = t79*t55
619 0 : t739 = t3*t151
620 0 : t746 = t37*t354
621 0 : t769 = t40*t573
622 : t791 = 0.5400000000e2_dp*t769*t43*t606 - 0.3600000000e2_dp*t313* &
623 : t132*t212 - 0.5400000000e2_dp*t451*t452*t258 - 0.6000000000e1_dp &
624 : *t128*t323*t94 + 0.1800000000e2_dp*t128*t132*t258 + 0.8999999998e1_dp &
625 0 : *t128*t43*t677 - 0.2666666667e1_dp*t42*pi*t79
626 0 : t793 = t328*t135
627 : t838 = REAL(3*t326*t135*t46, KIND=dp) + REAL(t791*t46, KIND=dp) + REAL(t793* &
628 : t46, KIND=dp) - 0.5555555558e-1_dp*t39*t677*t52 + 0.1111111112e0_dp*t330 &
629 : *t143 - 0.1666666668e0_dp*t330*t48*t148 - 0.1851851853e0_dp*t137 &
630 : *t339 + 0.2222222223e0_dp*t335*t343 - 0.1666666668e0_dp*t335* &
631 : t347 - 0.1666666668e0_dp*t335*t351 + 0.1646090535e0_dp*t47*t48 &
632 : *t71*t51 - 0.1851851853e0_dp*REAL(t146, KIND=dp)*REAL(t10, KIND=dp)*REAL(t135, KIND=dp) &
633 : *REAL(t46, KIND=dp) + 0.1111111112e0_dp*REAL(t146, KIND=dp)*REAL(t141, KIND=dp)*REAL(t326, KIND=dp) &
634 : *REAL(t46, KIND=dp) + 0.1111111112e0_dp*REAL(t146, KIND=dp)*REAL(t141, KIND=dp)*REAL(t328, KIND=dp) &
635 : *REAL(t46, KIND=dp) - 0.5555555558e-1_dp*REAL(t146, KIND=dp)*REAL(t49, KIND=dp)*REAL(t791, KIND=dp) &
636 : *REAL(t46, KIND=dp) - 0.1666666668e0_dp*REAL(t146, KIND=dp)*REAL(t346, KIND=dp)*REAL(t136, KIND=dp) &
637 : - 0.5555555558e-1_dp*REAL(t146, KIND=dp)*REAL(t49, KIND=dp)*REAL(t793, KIND=dp)* &
638 0 : REAL(t46, KIND=dp)
639 : t842 = 0.2e1_dp*t101*(-t316 + t319 + 0.9000000000e1_dp*t321 + t325) &
640 : *t103*t112 + 0.2e1_dp*t101*t103*(-0.5625000000e1_dp*t687*t31 &
641 : *t606 + 0.2250000000e1_dp*t271*t109*t212 + 0.6750000000e1_dp* &
642 : t417*t418*t258 + 0.1000000000e1_dp*t105*t281*t94 - 0.1500000000e1_dp &
643 : *t105*t109*t258 - 0.1500000000e1_dp*t105*t31*t677 &
644 : + 0.1111111111e1_dp*t30*t26*t10) + 0.4e1_dp*t101*t265*t284 + &
645 : 0.2e1_dp*t101*t716*t103*t112 + 0.1250000000e0_dp*t722*t118 &
646 : *t606 + 0.8333333333e-1_dp*t289*t293*t212 - 0.2500000000e0_dp*t289 &
647 : *t297*t212 - 0.2500000000e0_dp*t289*t118*t613 + 0.2222222222e0_dp &
648 : *t117*t735*t94 - 0.3333333333e0_dp*t117*t739*t94 - &
649 : 0.1666666667e0_dp*t117*t293*t258 + 0.5000000001e0_dp*t117*t746 &
650 : *t94 + 0.5000000001e0_dp*t117*t297*t258 + 0.1666666667e0_dp*t117 &
651 : *t118*t677 - 0.3456790122e0_dp*t36*t27*t237*t55 + 0.4444444444e0_dp &
652 : *t36*t304*t151 - 0.3333333333e0_dp*t36*t122*t354 &
653 0 : + 0.3333333334e0_dp*t36*t38*t838
654 : t846 = -0.5000000004e0_dp*t116*t213 - 0.2000000001e1_dp*t36*t216 &
655 : - 0.1000000001e1_dp*t36*t259 - 0.6666666672e0_dp*t64*t601 + 0.8333333340e-1_dp &
656 : *t605*t60*t606 - 0.5000000004e0_dp*t211*t207*t212 &
657 : - 0.5000000004e0_dp*t211*t60*t613 - 0.1000000001e1_dp*t68* &
658 : t601*t94 - 0.1000000001e1_dp*t68*t207*t258 - 0.3333333336e0_dp* &
659 0 : t68*t60*t677 - 0.2222222224e0_dp*t24*t98*t842
660 :
661 0 : e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + t846*sx
662 :
663 0 : t857 = t27*t496
664 0 : t860 = t212*t172
665 0 : t867 = t94*t404
666 0 : t880 = t258*t172
667 : t933 = 0.3911111110e2_dp*t11*t222 - 0.1955555555e2_dp*t6*t628*t168 &
668 : + 0.2133333334e2_dp*t11*t226 - 0.2133333334e2_dp*t6*t71*t384 &
669 : + 0.1066666667e2_dp*t6*t225*t400 + 0.80e1_dp*t11*t233 - 0.120e2_dp &
670 : *t382*t640*t232*t168 + 0.80e1_dp*t382*t383*t400 - 0.40e1_dp &
671 : *t11*t255 + 0.40e1_dp*t382*t230*t254*t168 - 0.20e1_dp* &
672 : t6*t77*(0.1866666667e2_dp*beta*t237*t12 + 0.9866666667e2_dp* &
673 : t11*t241 - 0.8266666668e2_dp*t393*t251 + 0.3200000001e2_dp*beta &
674 0 : *t244*my_ndrho/t657/t7*t669)
675 0 : t961 = t687*t26
676 0 : t1002 = t3*t196
677 : t1009 = 0.2e1_dp*t101*(-t455 + 0.9000000000e1_dp*t457 + 0.6000000000e1_dp &
678 : *t460)*t103*t112 + 0.1800000000e2_dp*t412*t264*t40*t127 &
679 : *t8*t172*t103*t112 + 0.2e1_dp*t101*t265*t428 + 0.1800000000e2_dp &
680 : *t413*t186*t285 + 0.2e1_dp*t101*t103*(-0.5625000000e1_dp &
681 : *t961*t1*t212*t172 + 0.4500000000e1_dp*t417*t418*t404 &
682 : + 0.1500000000e1_dp*t417*t49*t94*t172 - 0.1000000000e1_dp* &
683 : t105*t109*t404 + 0.2250000000e1_dp*t417*t1*t258*t172 - 0.1500000000e1_dp &
684 : *t105*t31*t933 + 0.3333333334e0_dp*t105*t281* &
685 : t172) + 0.1250000000e0_dp*t722*t118*t860 - 0.8333333335e-1_dp*t289 &
686 : *t435*t212 - 0.1666666667e0_dp*t289*t118*t867 + 0.5555555555e-1_dp &
687 : *t289*t293*t368 - 0.1111111111e0_dp*t117*t1002*t94 &
688 0 : - 0.1111111111e0_dp*t117*t293*t404
689 0 : t1013 = t37*t492
690 0 : t1044 = t769*pi
691 : t1069 = 0.5400000000e2_dp*t1044*t8*t212*t172 - 0.3600000000e2_dp &
692 : *t451*t452*t404 - 0.2400000000e2_dp*t451*t37*t94*t172 + &
693 : 0.1200000000e2_dp*t128*t132*t404 - 0.1800000000e2_dp*t451*t8* &
694 : t258*t172 + 0.8999999998e1_dp*t128*t43*t933 - 0.2000000000e1_dp &
695 0 : *t128*t323*t172
696 0 : t1091 = t127*t172*t46
697 : t1102 = t1069*t46 + 0.8999999998e1_dp*t326*t40*t127*t467 + REAL(2 &
698 : *t136*t462, KIND=dp) + 0.8999999998e1_dp*t328*t40*t127*t467 - &
699 : 0.5555555558e-1_dp*t39*t933*t52 - 0.5000000001e0_dp*t258*t127 &
700 : *t466 + 0.7407407410e-1_dp*t470*t143 + 0.6666666668e0_dp*t94*t478 &
701 : *t1091 - 0.1111111112e0_dp*t470*t48*t148 - 0.1111111112e0_dp &
702 0 : *t335*t486 - 0.1000000001e1_dp*t94*t135*t1091
703 : t1136 = -0.6172839508e-1_dp*t190*t339 - 0.5555555556e0_dp*t41/t7 &
704 : *t466 + 0.7407407410e-1_dp*t482*t343 + 0.7407407410e-1_dp*t146* &
705 : t141*t462*t46 + 0.6666666668e0_dp*t479*t135*t172*t46 - 0.5555555558e-1_dp &
706 : *t482*t347 - 0.5555555558e-1_dp*t146*t49*t1069 &
707 : *t46 - 0.5000000001e0_dp*t41*t326*t466 - 0.5555555558e-1_dp*t482 &
708 : *t351 - 0.1111111112e0_dp*t146*t147*t463 - 0.5000000001e0_dp* &
709 0 : t41*t328*t466
710 : t1141 = -0.1666666667e0_dp*t289*t297*t368 + 0.3333333334e0_dp*t117 &
711 : *t1013*t94 + 0.3333333334e0_dp*t117*t297*t404 - 0.8333333335e-1_dp &
712 : *t289*t118*t880 + 0.1666666667e0_dp*t117*t435*t258 + &
713 : 0.1666666667e0_dp*t117*t118*t933 + 0.7407407405e-1_dp*t117*t735 &
714 : *t172 + 0.1481481481e0_dp*t36*t304*t196 - 0.1111111111e0_dp* &
715 : t117*t739*t172 - 0.2222222222e0_dp*t36*t122*t492 + 0.1666666667e0_dp &
716 : *t117*t746*t172 + 0.3333333334e0_dp*t36*t38*(t1102 &
717 0 : + t1136)
718 : t1146 = -0.3333333336e0_dp*t117*t59*t94*t172 - 0.6666666672e0_dp &
719 : *t36*t372 - 0.6666666672e0_dp*t36*t405 - 0.6666666672e0_dp*t36 &
720 : *t408 - 0.4444444448e0_dp*t64*t857 + 0.8333333340e-1_dp*t605*t60 &
721 : *t860 - 0.1666666668e0_dp*t211*t365*t212 - 0.3333333336e0_dp* &
722 : t211*t60*t867 - 0.3333333336e0_dp*t211*t207*t368 - 0.6666666672e0_dp &
723 : *t68*t857*t94 - 0.6666666672e0_dp*t68*t207*t404 - 0.1666666668e0_dp &
724 : *t211*t60*t880 - 0.3333333336e0_dp*t68*t365* &
725 : t258 - 0.3333333336e0_dp*t68*t60*t933 - 0.3333333336e0_dp*t68* &
726 0 : t601*t172 - 0.2222222224e0_dp*t24*t98*(t1009 + t1141)
727 :
728 0 : e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) + t1146*sx
729 :
730 0 : t1153 = t27*t590
731 0 : t1156 = t94*t501
732 0 : t1163 = t404*t172
733 0 : t1167 = t94*t529
734 0 : t1177 = beta*t71
735 : t1220 = -0.1066666667e2_dp*t1177*t17 + 0.2133333334e2_dp*t11*t377 &
736 : - 0.1066666667e2_dp*t6*t632*t513 + 0.5333333333e1_dp*t6*t225 &
737 : *t525 - 0.40e1_dp*t508*t76*t90 + 0.160e2_dp*t11*t10*t384 - &
738 : 0.80e1_dp*t11*t401 - 0.120e2_dp*t382*t640*t90*t513 + 0.80e1_dp &
739 : *t382*t230*t400*t168 + 0.40e1_dp*t382*t383*t525 - 0.20e1_dp &
740 : *t6*t77*(-0.3200000000e2_dp*t1177*t86 + 0.4800000000e2_dp*t6 &
741 0 : *t397 - 0.2400000000e2_dp*t245/t657/my_rho*t669)
742 0 : t1284 = t37*t586
743 : t1334 = 0.5400000000e2_dp*t1044*t452*t501 - 0.3600000000e2_dp*t451 &
744 : *t8*t404*t172 - 0.1800000000e2_dp*t451*t452*t529 + 0.8999999998e1_dp &
745 : *t128*t43*t1220 - 0.1200000000e2_dp*t313*t132*t501 &
746 0 : + 0.5999999999e1_dp*t128*t132*t529
747 0 : t1341 = t501*t46
748 0 : t1345 = t529*t46
749 0 : t1370 = t40*pi
750 : t1396 = t1334*t46 + 0.1800000000e2_dp*t462*t40*t127*t467 - 0.2250000000e2_dp &
751 : *t464*t312*t43*t1341 + 0.8999999998e1_dp*t465 &
752 : *t43*t1345 - 0.1000000000e1_dp*t404*t127*t466 + 0.8099999996e2_dp &
753 : *t135*t571*t573*t533*t2*t1341 - 0.5555555558e-1_dp*t39 &
754 : *t1220*t52 - 0.5000000001e0_dp*t473*t1345 + 0.3333333334e0_dp* &
755 : t479*t1345 + 0.1000000000e1_dp*t94*t312*t1341 - 0.4500000000e1_dp &
756 : *t94*t573*t501*t1370*t8*t46 + 0.3703703705e-1_dp*t580 &
757 : *t143 + 0.3000000000e1_dp*t312*t37*t501*t1370*t46 - 0.5555555558e-1_dp &
758 : *t580*t48*t148 - 0.1111111112e0_dp*t482*t486 - 0.5555555558e-1_dp &
759 : *t146*t49*t1334*t46 - 0.1000000000e1_dp*t41* &
760 0 : t462*t466 - 0.5000000001e0_dp*t489*t1345
761 : t1400 = -0.3600000000e2_dp*t412*t313*t562*t113 + 0.1800000000e2_dp &
762 : *t413*t566*t113 + 0.3600000000e2_dp*t413*t186*t429 + 0.1620000000e3_dp &
763 : *t26*t533*t100*t574*t576*t113 + 0.2e1_dp*t101 &
764 : *t103*(-0.5625000000e1_dp*t961*t418*t501 + 0.4500000000e1_dp &
765 : *t417*t1*t404*t172 + 0.2250000000e1_dp*t417*t418*t529 - &
766 : 0.1500000000e1_dp*t105*t31*t1220 + 0.7500000000e0_dp*t271*t109 &
767 : *t501 - 0.5000000000e0_dp*t105*t109*t529) + 0.1250000000e0_dp* &
768 : t722*t118*t1156 - 0.1666666667e0_dp*t289*t435*t368 - 0.1666666667e0_dp &
769 : *t289*t118*t1163 - 0.8333333335e-1_dp*t289*t118*t1167 &
770 : + 0.1666666667e0_dp*t117*t1284*t94 + 0.3333333334e0_dp*t117 &
771 : *t435*t404 + 0.1666666667e0_dp*t117*t118*t1220 + 0.2777777778e-1_dp &
772 : *t289*t293*t501 - 0.1111111111e0_dp*t117*t1002*t172 &
773 : - 0.5555555555e-1_dp*t117*t293*t529 - 0.1111111111e0_dp*t36*t122 &
774 : *t586 - 0.8333333335e-1_dp*t289*t297*t501 + 0.3333333334e0_dp &
775 : *t117*t1013*t172 + 0.1666666667e0_dp*t117*t297*t529 + 0.3333333334e0_dp &
776 0 : *t36*t38*t1396
777 : t1404 = -0.1666666668e0_dp*t116*t502 - 0.6666666672e0_dp*t36*t505 &
778 : - 0.3333333336e0_dp*t36*t530 - 0.2222222224e0_dp*t64*t1153 + 0.8333333340e-1_dp &
779 : *t605*t60*t1156 - 0.3333333336e0_dp*t211*t365 &
780 : *t368 - 0.3333333336e0_dp*t211*t60*t1163 - 0.1666666668e0_dp*t211 &
781 : *t60*t1167 - 0.3333333336e0_dp*t68*t1153*t94 - 0.6666666672e0_dp &
782 : *t68*t365*t404 - 0.3333333336e0_dp*t68*t60*t1220 - 0.1666666668e0_dp &
783 : *t211*t207*t501 - 0.6666666672e0_dp*t68*t857* &
784 : t172 - 0.3333333336e0_dp*t68*t207*t529 - 0.2222222224e0_dp*t24* &
785 0 : t98*t1400
786 :
787 0 : e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) + t1404*sx
788 :
789 0 : t1405 = t501*t172
790 0 : t1412 = t172*t529
791 : t1449 = -0.120e2_dp*t508*t76*t168 + 0.240e2_dp*t11*t514 - 0.120e2_dp &
792 : *t11*t526 - 0.120e2_dp*t6*t641*t513*t168 + 0.120e2_dp*t382 &
793 : *t230*t168*t525 - 0.20e1_dp*t6*t77*(-0.240e2_dp*beta*t521 &
794 0 : *t250*my_ndrho + 0.180e2_dp*t393/t657*t669)
795 0 : t1456 = t1405*t103
796 0 : t1467 = t533*pi
797 0 : t1472 = t572*t21
798 : t1553 = 0.1350000000e3_dp*t537/t22/t572*my_rho*t1456 - 0.8100000000e2_dp &
799 : *t534*t536*t539*my_rho*t172*t103*t529 - 0.2430000000e3_dp &
800 : *t1467*t100/t570/omega/t22/t1472*t140*t1456 - &
801 : 0.1125000000e2_dp*t177*t687*t1*t1405 + 0.1350000000e2_dp*t176 &
802 : *t103*t28*t270*t1*t1412 - 0.3000000000e1_dp*t177*t105* &
803 : t1*t1449 + 0.1250000000e0_dp*t722*t118*t1405 - 0.2500000000e0_dp &
804 : *t289*t435*t501 - 0.2500000000e0_dp*t289*t118*t1412 + 0.5000000001e0_dp &
805 : *t117*t1284*t172 + 0.5000000001e0_dp*t117*t435 &
806 : *t529 + 0.1666666667e0_dp*t117*t118*t1449 + 0.3333333334e0_dp*t36 &
807 : *t38*(0.6750000000e2_dp*t1044*t8*t1405*t46 - 0.6750000000e2_dp &
808 : *t451*t186*t1345 - 0.5264999998e3_dp*t571/t1472*t533 &
809 : *t2*t1405*t46 + 0.8999999998e1_dp*t185*t8*t1449*t46 + 0.2429999999e3_dp &
810 : *t575*t2*t529*t466 + 0.7289999995e3_dp/t570/t39 &
811 : /t572/t126*t1467*t7*t1405*t46 - 0.5555555558e-1_dp*t39 &
812 0 : *t1449*t52 - 0.5000000001e0_dp*t41*t1449*t46)
813 :
814 : e_ndrho_ndrho_ndrho(ii) = e_ndrho_ndrho_ndrho(ii) + (0.8333333340e-1_dp*t605*t60*t1405 - &
815 : 0.5000000004e0_dp*t211*t365*t501 - 0.5000000004e0_dp*t211*t60*t1412 - 0.1000000001e1_dp &
816 : *t68*t1153*t172 - 0.1000000001e1_dp*t68*t365*t529 - 0.3333333336e0_dp &
817 : *t68*t60*t1449 - 0.2222222224e0_dp*t24*t98*t1553) &
818 0 : *sx
819 :
820 : END IF
821 : END IF
822 : END DO
823 : !$OMP END DO
824 :
825 1224 : END SUBROUTINE xb88_lr_lda_calc
826 :
827 : ! **************************************************************************************************
828 : !> \brief evaluates the becke 88 longrange exchange functional for lsd
829 : !> \param rho_set ...
830 : !> \param deriv_set ...
831 : !> \param grad_deriv ...
832 : !> \param xb88_lr_params ...
833 : !> \par History
834 : !> 01.2008 created [mguidon]
835 : !> \author Manuel Guidon
836 : ! **************************************************************************************************
837 260 : SUBROUTINE xb88_lr_lsd_eval(rho_set, deriv_set, grad_deriv, xb88_lr_params)
838 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
839 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
840 : INTEGER, INTENT(in) :: grad_deriv
841 : TYPE(section_vals_type), POINTER :: xb88_lr_params
842 :
843 : CHARACTER(len=*), PARAMETER :: routineN = 'xb88_lr_lsd_eval'
844 :
845 : INTEGER :: handle, i, ispin, npoints
846 : INTEGER, DIMENSION(2, 3) :: bo
847 : REAL(kind=dp) :: epsilon_rho, omega, sx
848 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
849 52 : POINTER :: dummy, e_0
850 468 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
851 936 : e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, e_rho_rho_rho, &
852 312 : norm_drho, rho
853 : TYPE(xc_derivative_type), POINTER :: deriv
854 :
855 52 : CALL timeset(routineN, handle)
856 :
857 52 : CALL cite_reference(Becke1988)
858 :
859 52 : NULLIFY (deriv)
860 156 : DO i = 1, 2
861 156 : NULLIFY (norm_drho(i)%array, rho(i)%array)
862 : END DO
863 :
864 52 : CALL section_vals_val_get(xb88_lr_params, "scale_x", r_val=sx)
865 52 : CALL section_vals_val_get(xb88_lr_params, "omega", r_val=omega)
866 :
867 : CALL xc_rho_set_get(rho_set, &
868 : rhoa=rho(1)%array, &
869 : rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
870 : norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
871 52 : local_bounds=bo)
872 52 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
873 :
874 : ! meaningful default for the arrays we don't need: let us make compiler
875 : ! and debugger happy...
876 52 : dummy => rho(1)%array
877 :
878 52 : e_0 => dummy
879 156 : DO i = 1, 2
880 104 : e_rho(i)%array => dummy
881 104 : e_ndrho(i)%array => dummy
882 104 : e_rho_rho(i)%array => dummy
883 104 : e_ndrho_rho(i)%array => dummy
884 104 : e_ndrho_ndrho(i)%array => dummy
885 104 : e_rho_rho_rho(i)%array => dummy
886 104 : e_ndrho_rho_rho(i)%array => dummy
887 104 : e_ndrho_ndrho_rho(i)%array => dummy
888 156 : e_ndrho_ndrho_ndrho(i)%array => dummy
889 : END DO
890 :
891 52 : IF (grad_deriv >= 0) THEN
892 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
893 52 : allocate_deriv=.TRUE.)
894 52 : CALL xc_derivative_get(deriv, deriv_data=e_0)
895 : END IF
896 52 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
897 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
898 52 : allocate_deriv=.TRUE.)
899 52 : CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
900 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
901 52 : allocate_deriv=.TRUE.)
902 52 : CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
903 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
904 52 : allocate_deriv=.TRUE.)
905 52 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
906 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
907 52 : allocate_deriv=.TRUE.)
908 52 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
909 : END IF
910 52 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
911 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
912 0 : allocate_deriv=.TRUE.)
913 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(1)%array)
914 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
915 0 : allocate_deriv=.TRUE.)
916 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(2)%array)
917 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
918 0 : allocate_deriv=.TRUE.)
919 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(1)%array)
920 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
921 0 : allocate_deriv=.TRUE.)
922 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(2)%array)
923 : deriv => xc_dset_get_derivative(deriv_set, &
924 0 : [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
925 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(1)%array)
926 : deriv => xc_dset_get_derivative(deriv_set, &
927 0 : [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
928 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(2)%array)
929 : END IF
930 52 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
931 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
932 0 : allocate_deriv=.TRUE.)
933 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(1)%array)
934 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
935 0 : allocate_deriv=.TRUE.)
936 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(2)%array)
937 : deriv => xc_dset_get_derivative(deriv_set, &
938 0 : [deriv_norm_drhoa, deriv_rhoa, deriv_rhoa], allocate_deriv=.TRUE.)
939 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(1)%array)
940 : deriv => xc_dset_get_derivative(deriv_set, &
941 0 : [deriv_norm_drhob, deriv_rhob, deriv_rhob], allocate_deriv=.TRUE.)
942 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(2)%array)
943 : deriv => xc_dset_get_derivative(deriv_set, &
944 0 : [deriv_norm_drhoa, deriv_norm_drhoa, deriv_rhoa], allocate_deriv=.TRUE.)
945 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(1)%array)
946 : deriv => xc_dset_get_derivative(deriv_set, &
947 0 : [deriv_norm_drhob, deriv_norm_drhob, deriv_rhob], allocate_deriv=.TRUE.)
948 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(2)%array)
949 : deriv => xc_dset_get_derivative(deriv_set, &
950 0 : [deriv_norm_drhoa, deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
951 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(1)%array)
952 : deriv => xc_dset_get_derivative(deriv_set, &
953 0 : [deriv_norm_drhob, deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
954 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(2)%array)
955 : END IF
956 52 : IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
957 0 : CPABORT("derivatives bigger than 3 not implemented")
958 : END IF
959 :
960 156 : DO ispin = 1, 2
961 :
962 : !$OMP PARALLEL DEFAULT(NONE) &
963 : !$OMP SHARED(rho, norm_drho, e_0, e_rho, e_ndrho) &
964 : !$OMP SHARED(e_rho_rho, e_ndrho_rho, e_ndrho_ndrho) &
965 : !$OMP SHARED(e_rho_rho_rho, e_ndrho_rho_rho) &
966 : !$OMP SHARED(e_ndrho_ndrho_rho, e_ndrho_ndrho_ndrho) &
967 : !$OMP SHARED(grad_deriv, npoints, epsilon_rho, sx, omega) &
968 156 : !$OMP SHARED(ispin)
969 :
970 : CALL xb88_lr_lsd_calc( &
971 : rho_spin=rho(ispin)%array, &
972 : norm_drho_spin=norm_drho(ispin)%array, &
973 : e_0=e_0, &
974 : e_rho_spin=e_rho(ispin)%array, &
975 : e_ndrho_spin=e_ndrho(ispin)%array, &
976 : e_rho_rho_spin=e_rho_rho(ispin)%array, &
977 : e_ndrho_rho_spin=e_ndrho_rho(ispin)%array, &
978 : e_ndrho_ndrho_spin=e_ndrho_ndrho(ispin)%array, &
979 : e_rho_rho_rho_spin=e_rho_rho_rho(ispin)%array, &
980 : e_ndrho_rho_rho_spin=e_ndrho_rho_rho(ispin)%array, &
981 : e_ndrho_ndrho_rho_spin=e_ndrho_ndrho_rho(ispin)%array, &
982 : e_ndrho_ndrho_ndrho_spin=e_ndrho_ndrho_ndrho(ispin)%array, &
983 : grad_deriv=grad_deriv, npoints=npoints, &
984 : epsilon_rho=epsilon_rho, &
985 : sx=sx, omega=omega)
986 :
987 : !$OMP END PARALLEL
988 :
989 : END DO
990 :
991 52 : CALL timestop(handle)
992 :
993 52 : END SUBROUTINE xb88_lr_lsd_eval
994 :
995 : ! **************************************************************************************************
996 : !> \brief low level calculation of the becke 88 exchange functional for lsd
997 : !> \param rho_spin alpha or beta spin density
998 : !> \param norm_drho_spin || grad rho_spin ||
999 : !> \param e_0 adds to it the local value of the functional
1000 : !> \param e_rho_spin e_*_spin derivative of the functional wrt. to the variables
1001 : !> named where the * is. Everything wrt. to the spin of the arguments.
1002 : !> \param e_ndrho_spin ...
1003 : !> \param e_rho_rho_spin ...
1004 : !> \param e_ndrho_rho_spin ...
1005 : !> \param e_ndrho_ndrho_spin ...
1006 : !> \param e_rho_rho_rho_spin ...
1007 : !> \param e_ndrho_rho_rho_spin ...
1008 : !> \param e_ndrho_ndrho_rho_spin ...
1009 : !> \param e_ndrho_ndrho_ndrho_spin ...
1010 : !> \param grad_deriv ...
1011 : !> \param npoints ...
1012 : !> \param epsilon_rho ...
1013 : !> \param sx scaling-parameter for exchange
1014 : !> \param omega ...
1015 : !> \par History
1016 : !> 01.2008 created [mguidon]
1017 : !> \author Manuel Guidon
1018 : !> \note
1019 : !> - Derivatives higher than one not tested
1020 : ! **************************************************************************************************
1021 104 : SUBROUTINE xb88_lr_lsd_calc(rho_spin, norm_drho_spin, e_0, &
1022 : e_rho_spin, e_ndrho_spin, e_rho_rho_spin, e_ndrho_rho_spin, &
1023 : e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
1024 : e_ndrho_ndrho_rho_spin, &
1025 : e_ndrho_ndrho_ndrho_spin, grad_deriv, npoints, epsilon_rho, sx, &
1026 : omega)
1027 : REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho_spin, norm_drho_spin
1028 : REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho_spin, e_ndrho_spin, e_rho_rho_spin, &
1029 : e_ndrho_rho_spin, e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
1030 : e_ndrho_ndrho_rho_spin, e_ndrho_ndrho_ndrho_spin
1031 : INTEGER, INTENT(in) :: grad_deriv, npoints
1032 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, omega
1033 :
1034 : INTEGER :: ii
1035 : REAL(kind=dp) :: Cx, epsilon_rho43, my_epsilon_rho, ndrho, rho, t1, t10, t100, t1002, t1009, &
1036 : t101, t1013, t103, t104, t1044, t105, t1069, t109, t1091, t11, t1102, t112, t113, t1136, &
1037 : t1141, t1146, t1153, t1156, t116, t1163, t1167, t117, t1177, t118, t12, t122, t1220, &
1038 : t126, t127, t128, t1284, t130, t132, t133, t1334, t1341, t1345, t135, t136, t137, t1370, &
1039 : t1396, t140, t1400, t1404, t1405, t141, t1412, t143, t1449, t1456, t146, t1467, t147, &
1040 : t1472, t148, t151, t155, t1553, t16, t168, t169, t17, t172, t173, t176, t177, t18, t185, &
1041 : t186, t190, t196, t2, t200, t207, t21, t211, t212, t213, t216
1042 : REAL(kind=dp) :: t219, t22, t221, t222, t225, t226, t23, t230, t231, t232, t233, t237, t24, &
1043 : t241, t244, t245, t250, t251, t254, t255, t258, t259, t26, t264, t265, t27, t270, t271, &
1044 : t28, t281, t284, t285, t289, t29, t293, t297, t3, t30, t304, t31, t311, t312, t313, t316, &
1045 : t319, t321, t323, t325, t326, t328, t330, t335, t339, t34, t343, t346, t347, t351, t354, &
1046 : t358, t36, t365, t368, t37, t372, t377, t38, t382, t383, t384, t39, t393, t397, t40, &
1047 : t400, t401, t404, t405, t408, t41, t412, t413, t417, t418, t42, t428, t429, t43, t435, &
1048 : t44, t451, t452, t455, t457, t46, t460, t462, t463, t464, t465
1049 : REAL(kind=dp) :: t466, t467, t47, t470, t473, t478, t479, t48, t482, t486, t489, t49, t492, &
1050 : t496, t5, t501, t502, t505, t508, t51, t513, t514, t519, t52, t521, t525, t526, t529, &
1051 : t530, t533, t534, t536, t537, t539, t55, t562, t566, t570, t571, t572, t573, t574, t575, &
1052 : t576, t580, t586, t59, t590, t6, t60, t601, t605, t606, t613, t624, t628, t632, t639, &
1053 : t64, t640, t641, t657, t667, t669, t677, t68, t687, t69, t7, t71, t716, t72, t722, t735, &
1054 : t739, t746, t75, t76, t769, t77, t79, t791, t793, t8, t838, t84, t842, t846, t85, t857, &
1055 : t86, t860, t867, t87, t880, t90, t91, t933, t94, t95, t961, t98
1056 : REAL(kind=dp) :: t99, xx
1057 :
1058 104 : my_epsilon_rho = 0.5_dp*epsilon_rho
1059 104 : epsilon_rho43 = my_epsilon_rho**(4.0_dp/3.0_dp)
1060 104 : Cx = 1.5_dp*(3.0_dp/4.0_dp/pi)**(1.0_dp/3.0_dp)
1061 :
1062 104 : !$OMP DO
1063 : DO ii = 1, npoints
1064 4224064 : rho = rho_spin(ii)
1065 4224064 : ndrho = norm_drho_spin(ii)
1066 4224064 : IF (rho > my_epsilon_rho) THEN
1067 4216664 : IF (grad_deriv >= 0) THEN
1068 4216664 : t1 = rho**(0.1e1_dp/0.3e1_dp)
1069 4216664 : t2 = t1*rho
1070 4216664 : t3 = 0.1e1_dp/t2
1071 4216664 : xx = ndrho*MAX(t3, epsilon_rho43)
1072 4216664 : t5 = ndrho**2
1073 4216664 : t6 = beta*t5
1074 4216664 : t7 = rho**2
1075 4216664 : t8 = t1**2
1076 4216664 : t10 = 0.1e1_dp/t8/t7
1077 4216664 : t11 = beta*ndrho
1078 4216664 : t12 = LOG(xx + SQRT(xx**0.2e1_dp + 0.1e1_dp))
1079 4216664 : t16 = 0.10e1_dp + 0.60e1_dp*t11*t3*t12
1080 4216664 : t17 = 0.1e1_dp/t16
1081 4216664 : t18 = t10*t17
1082 4216664 : t21 = 0.20e1_dp*Cx + 0.20e1_dp*t6*t18
1083 4216664 : t22 = SQRT(t21)
1084 4216664 : t23 = t22*t21
1085 4216664 : t24 = rho*t23
1086 4216664 : t26 = rootpi
1087 4216664 : t27 = 0.1e1_dp/t26
1088 4216664 : t28 = 0.1e1_dp/omega
1089 4216664 : t29 = 0.1e1_dp/t22
1090 4216664 : t30 = t28*t29
1091 4216664 : t31 = t26*t1
1092 4216664 : t34 = erf(0.3000000000e1_dp*t30*t31)
1093 4216664 : t36 = omega*t22
1094 4216664 : t37 = 0.1e1_dp/t1
1095 4216664 : t38 = t27*t37
1096 4216664 : t39 = omega**2
1097 4216664 : t40 = 0.1e1_dp/t39
1098 4216664 : t41 = 0.1e1_dp/t21
1099 4216664 : t42 = t40*t41
1100 4216664 : t43 = pi*t8
1101 4216664 : t44 = t42*t43
1102 4216664 : t46 = EXP(-0.8999999998e1_dp*t44)
1103 4216664 : t47 = t39*t21
1104 4216664 : t48 = 0.1e1_dp/pi
1105 4216664 : t49 = 0.1e1_dp/t8
1106 4216664 : t51 = t46 - 0.10e1_dp
1107 4216664 : t52 = t48*t49*t51
1108 4216664 : t55 = t46 - 0.15e1_dp - 0.5555555558e-1_dp*t47*t52
1109 4216664 : t59 = t26*t34 + 0.3333333334e0_dp*t36*t38*t55
1110 4216664 : t60 = t27*t59
1111 :
1112 4216664 : e_0(ii) = e_0(ii) - 0.2222222224e0_dp*t24*omega*t60*sx
1113 :
1114 : END IF
1115 4216664 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1116 4216664 : t64 = t23*omega
1117 4216664 : t68 = rho*t22*omega
1118 4216664 : t69 = t7*rho
1119 4216664 : t71 = 0.1e1_dp/t8/t69
1120 4216664 : t72 = t71*t17
1121 4216664 : t75 = t16**2
1122 4216664 : t76 = 0.1e1_dp/t75
1123 4216664 : t77 = t10*t76
1124 4216664 : t79 = 0.1e1_dp/t1/t7
1125 4216664 : t84 = 1 + t5*t10
1126 4216664 : t85 = SQRT(t84)
1127 4216664 : t86 = 0.1e1_dp/t85
1128 4216664 : t87 = t71*t86
1129 4216664 : t90 = -0.8000000000e1_dp*t11*t79*t12 - 0.8000000000e1_dp*t6*t87
1130 4216664 : t91 = t77*t90
1131 4216664 : t94 = -0.5333333333e1_dp*t6*t72 - 0.20e1_dp*t6*t91
1132 4216664 : t95 = t60*t94
1133 4216664 : t98 = omega*t27
1134 4216664 : t99 = SQRT(0.3141592654e1_dp)
1135 4216664 : t100 = 0.1e1_dp/t99
1136 4216664 : t101 = t26*t100
1137 4216664 : t103 = EXP(-0.9000000000e1_dp*t44)
1138 4216664 : t104 = 0.1e1_dp/t23
1139 4216664 : t105 = t28*t104
1140 4216664 : t109 = t26*t49
1141 : t112 = -0.1500000000e1_dp*t105*t31*t94 + 0.1000000000e1_dp*t30* &
1142 4216664 : t109
1143 4216664 : t113 = t103*t112
1144 4216664 : t116 = omega*t29
1145 4216664 : t117 = t116*t27
1146 4216664 : t118 = t37*t55
1147 4216664 : t122 = t27*t3
1148 4216664 : t126 = t21**2
1149 4216664 : t127 = 0.1e1_dp/t126
1150 4216664 : t128 = t40*t127
1151 4216664 : t130 = t128*t43*t94
1152 4216664 : t132 = pi*t37
1153 4216664 : t133 = t42*t132
1154 4216664 : t135 = 0.8999999998e1_dp*t130 - 0.5999999999e1_dp*t133
1155 4216664 : t136 = t135*t46
1156 4216664 : t137 = t39*t94
1157 4216664 : t140 = t8*rho
1158 4216664 : t141 = 0.1e1_dp/t140
1159 4216664 : t143 = t48*t141*t51
1160 4216664 : t146 = t47*t48
1161 4216664 : t147 = t49*t135
1162 4216664 : t148 = t147*t46
1163 : t151 = t136 - 0.5555555558e-1_dp*t137*t52 + 0.3703703705e-1_dp*t47 &
1164 4216664 : *t143 - 0.5555555558e-1_dp*t146*t148
1165 : t155 = REAL(2*t101*t113, KIND=dp) + 0.1666666667e0_dp*t117*t118*t94 - &
1166 : 0.1111111111e0_dp*t36*t122*t55 + 0.3333333334e0_dp*t36*t38* &
1167 4216664 : t151
1168 :
1169 : e_rho_spin(ii) = e_rho_spin(ii) + (-0.2222222224e0_dp*t64*t60 - 0.3333333336e0_dp*t68*t95 - &
1170 4216664 : 0.2222222224e0_dp*t24*t98*t155)*sx
1171 :
1172 4216664 : t168 = 0.60e1_dp*beta*t3*t12 + 0.60e1_dp*t11*t10*t86
1173 4216664 : t169 = t77*t168
1174 4216664 : t172 = 0.40e1_dp*t11*t18 - 0.20e1_dp*t6*t169
1175 4216664 : t173 = t60*t172
1176 4216664 : t176 = pi*t100
1177 4216664 : t177 = t176*t103
1178 4216664 : t185 = t128*pi
1179 4216664 : t186 = t8*t172
1180 4216664 : t190 = t39*t172
1181 : t196 = 0.8999999998e1_dp*t185*t186*t46 - 0.5555555558e-1_dp*t190 &
1182 4216664 : *t52 - 0.5000000001e0_dp*t41*t172*t46
1183 : t200 = -0.3000000000e1_dp*t177*t105*t1*t172 + 0.1666666667e0_dp* &
1184 4216664 : t117*t118*t172 + 0.3333333334e0_dp*t36*t38*t196
1185 :
1186 : e_ndrho_spin(ii) = e_ndrho_spin(ii) + (-0.3333333336e0_dp*t68*t173 - 0.2222222224e0_dp*t24*t98 &
1187 4216664 : *t200)*sx
1188 :
1189 : END IF
1190 4216664 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1191 0 : t207 = t27*t155
1192 0 : t211 = rho*t29*omega
1193 0 : t212 = t94**2
1194 0 : t213 = t60*t212
1195 0 : t216 = t207*t94
1196 0 : t219 = t7**2
1197 0 : t221 = 0.1e1_dp/t8/t219
1198 0 : t222 = t221*t17
1199 0 : t225 = t71*t76
1200 0 : t226 = t225*t90
1201 0 : t230 = 0.1e1_dp/t75/t16
1202 0 : t231 = t10*t230
1203 0 : t232 = t90**2
1204 0 : t233 = t231*t232
1205 0 : t237 = 0.1e1_dp/t1/t69
1206 0 : t241 = t221*t86
1207 0 : t244 = t5**2
1208 0 : t245 = beta*t244
1209 0 : t250 = 0.1e1_dp/t85/t84
1210 0 : t251 = 0.1e1_dp/t1/t219/t69*t250
1211 : t254 = 0.1866666667e2_dp*t11*t237*t12 + 0.4000000000e2_dp*t6*t241 &
1212 0 : - 0.1066666667e2_dp*t245*t251
1213 0 : t255 = t77*t254
1214 : t258 = 0.1955555555e2_dp*t6*t222 + 0.1066666667e2_dp*t6*t226 + 0.40e1_dp &
1215 0 : *t6*t233 - 0.20e1_dp*t6*t255
1216 0 : t259 = t60*t258
1217 0 : t264 = 0.9000000000e1_dp*t130 - 0.6000000000e1_dp*t133
1218 0 : t265 = t264*t103
1219 0 : t270 = 0.1e1_dp/t22/t126
1220 0 : t271 = t28*t270
1221 0 : t281 = t26*t141
1222 : t284 = 0.2250000000e1_dp*t271*t31*t212 - 0.1000000000e1_dp*t105* &
1223 : t109*t94 - 0.1500000000e1_dp*t105*t31*t258 - 0.6666666667e0_dp &
1224 0 : *t30*t281
1225 0 : t285 = t103*t284
1226 0 : t289 = omega*t104*t27
1227 0 : t293 = t3*t55
1228 0 : t297 = t37*t151
1229 0 : t304 = t27*t79
1230 0 : t311 = t126*t21
1231 0 : t312 = 0.1e1_dp/t311
1232 0 : t313 = t40*t312
1233 0 : t316 = 0.1800000000e2_dp*t313*t43*t212
1234 0 : t319 = 0.1200000000e2_dp*t128*t132*t94
1235 0 : t321 = t128*t43*t258
1236 0 : t323 = pi*t3
1237 0 : t325 = 0.2000000000e1_dp*t42*t323
1238 0 : t326 = -t316 + t319 + 0.8999999998e1_dp*t321 + t325
1239 0 : t328 = t135**2
1240 0 : t330 = t39*t258
1241 0 : t335 = t137*t48
1242 0 : t339 = t48*t10*t51
1243 0 : t343 = t141*t135*t46
1244 0 : t346 = t49*t326
1245 0 : t347 = t346*t46
1246 0 : t351 = t49*t328*t46
1247 : t354 = t326*t46 + t328*t46 - 0.5555555558e-1_dp*t330*t52 + 0.7407407410e-1_dp &
1248 : *t137*t143 - 0.1111111112e0_dp*t335*t148 - 0.6172839508e-1_dp &
1249 : *t47*t339 + 0.7407407410e-1_dp*t146*t343 - 0.5555555558e-1_dp &
1250 0 : *t146*t347 - 0.5555555558e-1_dp*t146*t351
1251 : t358 = REAL(2*t101*t265*t112, KIND=dp) + REAL(2*t101*t285, KIND=dp) - 0.8333333335e-1_dp &
1252 : *t289*t118*t212 - 0.1111111111e0_dp*t117*t293*t94 &
1253 : + 0.3333333334e0_dp*t117*t297*t94 + 0.1666666667e0_dp*t117* &
1254 : t118*t258 + 0.1481481481e0_dp*t36*t304*t55 - 0.2222222222e0_dp* &
1255 0 : t36*t122*t151 + 0.3333333334e0_dp*t36*t38*t354
1256 :
1257 : e_rho_rho_spin(ii) = e_rho_rho_spin(ii) + (-0.6666666672e0_dp*t36*t95 - 0.4444444448e0_dp*t64*t207 - &
1258 : 0.1666666668e0_dp*t211*t213 - 0.6666666672e0_dp*t68*t216 - 0.3333333336e0_dp &
1259 0 : *t68*t259 - 0.2222222224e0_dp*t24*t98*t358)*sx
1260 :
1261 0 : t365 = t27*t200
1262 0 : t368 = t94*t172
1263 0 : t372 = t365*t94
1264 0 : t377 = t225*t168
1265 0 : t382 = t6*t10
1266 0 : t383 = t230*t90
1267 0 : t384 = t383*t168
1268 0 : t393 = beta*t5*ndrho
1269 0 : t397 = 0.1e1_dp/t1/t219/t7*t250
1270 : t400 = -0.8000000000e1_dp*beta*t79*t12 - 0.2400000000e2_dp*t11* &
1271 0 : t87 + 0.8000000000e1_dp*t393*t397
1272 0 : t401 = t77*t400
1273 : t404 = -0.1066666667e2_dp*t11*t72 + 0.5333333333e1_dp*t6*t377 - 0.40e1_dp &
1274 0 : *t11*t91 + 0.40e1_dp*t382*t384 - 0.20e1_dp*t6*t401
1275 0 : t405 = t60*t404
1276 0 : t408 = t207*t172
1277 0 : t412 = t26*pi*t100
1278 0 : t413 = t412*t128
1279 0 : t417 = t271*t26
1280 0 : t418 = t1*t94
1281 : t428 = 0.2250000000e1_dp*t417*t418*t172 - 0.1500000000e1_dp*t105 &
1282 0 : *t31*t404 - 0.5000000000e0_dp*t105*t109*t172
1283 0 : t429 = t103*t428
1284 0 : t435 = t37*t196
1285 0 : t451 = t313*pi
1286 0 : t452 = t8*t94
1287 0 : t455 = 0.1800000000e2_dp*t451*t452*t172
1288 0 : t457 = t128*t43*t404
1289 0 : t460 = t128*t132*t172
1290 0 : t462 = -t455 + 0.8999999998e1_dp*t457 + 0.5999999999e1_dp*t460
1291 0 : t463 = t462*t46
1292 0 : t464 = t135*t40
1293 0 : t465 = t464*t127
1294 0 : t466 = t172*t46
1295 0 : t467 = t43*t466
1296 0 : t470 = t39*t404
1297 0 : t473 = t94*t127
1298 0 : t478 = 0.1e1_dp/rho
1299 0 : t479 = t41*t478
1300 0 : t482 = t190*t48
1301 0 : t486 = t49*t462*t46
1302 0 : t489 = t41*t135
1303 : t492 = t463 + 0.8999999998e1_dp*t465*t467 - 0.5555555558e-1_dp*t470 &
1304 : *t52 - 0.5000000001e0_dp*t473*t466 + 0.3703703705e-1_dp*t190*t143 &
1305 : + 0.3333333334e0_dp*t479*t466 - 0.5555555558e-1_dp*t482*t148 &
1306 0 : - 0.5555555558e-1_dp*t146*t486 - 0.5000000001e0_dp*t489*t466
1307 : t496 = 0.1800000000e2_dp*t413*t186*t113 + REAL(2*t101*t429, KIND=dp) &
1308 : - 0.8333333335e-1_dp*t289*t118*t368 + 0.1666666667e0_dp*t117*t435 &
1309 : *t94 + 0.1666666667e0_dp*t117*t118*t404 - 0.5555555555e-1_dp &
1310 : *t117*t293*t172 - 0.1111111111e0_dp*t36*t122*t196 + 0.1666666667e0_dp &
1311 0 : *t117*t297*t172 + 0.3333333334e0_dp*t36*t38*t492
1312 :
1313 : e_ndrho_rho_spin(ii) = e_ndrho_rho_spin(ii) + (-0.3333333336e0_dp*t36*t173 - 0.2222222224e0_dp*t64*t365 - &
1314 : 0.1666666668e0_dp*t211*t60*t368 - 0.3333333336e0_dp*t68*t372 &
1315 : - 0.3333333336e0_dp*t68*t405 - 0.3333333336e0_dp*t68*t408 - 0.2222222224e0_dp &
1316 0 : *t24*t98*t496)*sx
1317 :
1318 0 : t501 = t172**2
1319 0 : t502 = t60*t501
1320 0 : t505 = t365*t172
1321 0 : t508 = beta*t10
1322 0 : t513 = t168**2
1323 0 : t514 = t231*t513
1324 0 : t519 = t219*rho
1325 0 : t521 = 0.1e1_dp/t1/t519
1326 0 : t525 = 0.120e2_dp*t508*t86 - 0.60e1_dp*t6*t521*t250
1327 0 : t526 = t77*t525
1328 : t529 = 0.40e1_dp*t508*t17 - 0.80e1_dp*t11*t169 + 0.40e1_dp*t6*t514 &
1329 0 : - 0.20e1_dp*t6*t526
1330 0 : t530 = t60*t529
1331 0 : t533 = pi**2
1332 0 : t534 = t533*t100
1333 0 : t536 = 0.1e1_dp/t39/omega
1334 0 : t537 = t534*t536
1335 0 : t539 = 0.1e1_dp/t22/t311
1336 0 : t562 = t8*t501
1337 0 : t566 = t8*t529
1338 0 : t570 = t39**2
1339 0 : t571 = 0.1e1_dp/t570
1340 0 : t572 = t126**2
1341 0 : t573 = 0.1e1_dp/t572
1342 0 : t574 = t571*t573
1343 0 : t575 = t574*t533
1344 0 : t576 = t2*t501
1345 0 : t580 = t39*t529
1346 : t586 = -0.2250000000e2_dp*t451*t562*t46 + 0.8999999998e1_dp*t185 &
1347 : *t566*t46 + 0.8099999996e2_dp*t575*t576*t46 - 0.5555555558e-1_dp &
1348 0 : *t580*t52 - 0.5000000001e0_dp*t41*t529*t46
1349 : t590 = -0.2700000000e2_dp*t537*t539*rho*t501*t103 + 0.4500000000e1_dp &
1350 : *t177*t271*t1*t501 - 0.3000000000e1_dp*t177*t105* &
1351 : t1*t529 - 0.8333333335e-1_dp*t289*t118*t501 + 0.3333333334e0_dp &
1352 : *t117*t435*t172 + 0.1666666667e0_dp*t117*t118*t529 + 0.3333333334e0_dp &
1353 0 : *t36*t38*t586
1354 :
1355 : e_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_spin(ii) + (-0.1666666668e0_dp*t211*t502 - 0.6666666672e0_dp*t68*t505 &
1356 : - 0.3333333336e0_dp*t68*t530 - 0.2222222224e0_dp*t24*t98*t590) &
1357 0 : *sx
1358 :
1359 : END IF
1360 4216664 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
1361 0 : t601 = t27*t358
1362 0 : t605 = rho*t104*omega
1363 0 : t606 = t212*t94
1364 0 : t613 = t94*t258
1365 0 : t624 = 0.1e1_dp/t8/t519
1366 0 : t628 = t221*t76
1367 0 : t632 = t71*t230
1368 0 : t639 = t75**2
1369 0 : t640 = 0.1e1_dp/t639
1370 0 : t641 = t10*t640
1371 0 : t657 = t219**2
1372 0 : t667 = t84**2
1373 0 : t669 = 0.1e1_dp/t85/t667
1374 : t677 = -0.9125925923e2_dp*t6*t624*t17 - 0.5866666667e2_dp*t6*t628 &
1375 : *t90 - 0.3200000001e2_dp*t6*t632*t232 + 0.1600000000e2_dp*t6 &
1376 : *t225*t254 - 0.120e2_dp*t6*t641*t232*t90 + 0.120e2_dp*t382 &
1377 : *t383*t254 - 0.20e1_dp*t6*t77*(-0.6222222223e2_dp*t11/t1/ &
1378 : t219*t12 - 0.2115555556e3_dp*t6*t624*t86 + 0.1315555556e3_dp* &
1379 : t245/t1/t657*t250 - 0.4266666668e2_dp*beta*t244*t5/t657 &
1380 0 : /t69*t669)
1381 0 : t687 = t28*t539
1382 0 : t716 = t264**2
1383 0 : t722 = omega*t270*t27
1384 0 : t735 = t79*t55
1385 0 : t739 = t3*t151
1386 0 : t746 = t37*t354
1387 0 : t769 = t40*t573
1388 : t791 = 0.5400000000e2_dp*t769*t43*t606 - 0.3600000000e2_dp*t313* &
1389 : t132*t212 - 0.5400000000e2_dp*t451*t452*t258 - 0.6000000000e1_dp &
1390 : *t128*t323*t94 + 0.1800000000e2_dp*t128*t132*t258 + 0.8999999998e1_dp &
1391 0 : *t128*t43*t677 - 0.2666666667e1_dp*t42*pi*t79
1392 0 : t793 = t328*t135
1393 : t838 = REAL(3*t326*t135*t46, KIND=dp) + REAL(t791*t46, KIND=dp) + REAL(t793* &
1394 : t46, KIND=dp) - 0.5555555558e-1_dp*t39*t677*t52 + 0.1111111112e0_dp*t330 &
1395 : *t143 - 0.1666666668e0_dp*t330*t48*t148 - 0.1851851853e0_dp*t137 &
1396 : *t339 + 0.2222222223e0_dp*t335*t343 - 0.1666666668e0_dp*t335* &
1397 : t347 - 0.1666666668e0_dp*t335*t351 + 0.1646090535e0_dp*t47*t48 &
1398 : *t71*t51 - 0.1851851853e0_dp*REAL(t146, KIND=dp)*REAL(t10, KIND=dp)*REAL(t135, KIND=dp) &
1399 : *REAL(t46, KIND=dp) + 0.1111111112e0_dp*REAL(t146, KIND=dp)*REAL(t141, KIND=dp)*REAL(t326, KIND=dp) &
1400 : *REAL(t46, KIND=dp) + 0.1111111112e0_dp*REAL(t146, KIND=dp)*REAL(t141, KIND=dp)*REAL(t328, KIND=dp) &
1401 : *REAL(t46, KIND=dp) - 0.5555555558e-1_dp*REAL(t146, KIND=dp)*REAL(t49, KIND=dp)*REAL(t791, KIND=dp) &
1402 : *REAL(t46, KIND=dp) - 0.1666666668e0_dp*REAL(t146, KIND=dp)*REAL(t346, KIND=dp)*REAL(t136, KIND=dp) &
1403 : - 0.5555555558e-1_dp*REAL(t146, KIND=dp)*REAL(t49, KIND=dp)*REAL(t793, KIND=dp)* &
1404 0 : REAL(t46, KIND=dp)
1405 : t842 = 0.2e1_dp*t101*(-t316 + t319 + 0.9000000000e1_dp*t321 + t325) &
1406 : *t103*t112 + 0.2e1_dp*t101*t103*(-0.5625000000e1_dp*t687*t31 &
1407 : *t606 + 0.2250000000e1_dp*t271*t109*t212 + 0.6750000000e1_dp* &
1408 : t417*t418*t258 + 0.1000000000e1_dp*t105*t281*t94 - 0.1500000000e1_dp &
1409 : *t105*t109*t258 - 0.1500000000e1_dp*t105*t31*t677 &
1410 : + 0.1111111111e1_dp*t30*t26*t10) + 0.4e1_dp*t101*t265*t284 + &
1411 : 0.2e1_dp*t101*t716*t103*t112 + 0.1250000000e0_dp*t722*t118 &
1412 : *t606 + 0.8333333333e-1_dp*t289*t293*t212 - 0.2500000000e0_dp*t289 &
1413 : *t297*t212 - 0.2500000000e0_dp*t289*t118*t613 + 0.2222222222e0_dp &
1414 : *t117*t735*t94 - 0.3333333333e0_dp*t117*t739*t94 - &
1415 : 0.1666666667e0_dp*t117*t293*t258 + 0.5000000001e0_dp*t117*t746 &
1416 : *t94 + 0.5000000001e0_dp*t117*t297*t258 + 0.1666666667e0_dp*t117 &
1417 : *t118*t677 - 0.3456790122e0_dp*t36*t27*t237*t55 + 0.4444444444e0_dp &
1418 : *t36*t304*t151 - 0.3333333333e0_dp*t36*t122*t354 &
1419 0 : + 0.3333333334e0_dp*t36*t38*t838
1420 : t846 = -0.5000000004e0_dp*t116*t213 - 0.2000000001e1_dp*t36*t216 &
1421 : - 0.1000000001e1_dp*t36*t259 - 0.6666666672e0_dp*t64*t601 + 0.8333333340e-1_dp &
1422 : *t605*t60*t606 - 0.5000000004e0_dp*t211*t207*t212 &
1423 : - 0.5000000004e0_dp*t211*t60*t613 - 0.1000000001e1_dp*t68* &
1424 : t601*t94 - 0.1000000001e1_dp*t68*t207*t258 - 0.3333333336e0_dp* &
1425 0 : t68*t60*t677 - 0.2222222224e0_dp*t24*t98*t842
1426 :
1427 0 : e_rho_rho_rho_spin(ii) = e_rho_rho_rho_spin(ii) + t846*sx
1428 :
1429 0 : t857 = t27*t496
1430 0 : t860 = t212*t172
1431 0 : t867 = t94*t404
1432 0 : t880 = t258*t172
1433 : t933 = 0.3911111110e2_dp*t11*t222 - 0.1955555555e2_dp*t6*t628*t168 &
1434 : + 0.2133333334e2_dp*t11*t226 - 0.2133333334e2_dp*t6*t71*t384 &
1435 : + 0.1066666667e2_dp*t6*t225*t400 + 0.80e1_dp*t11*t233 - 0.120e2_dp &
1436 : *t382*t640*t232*t168 + 0.80e1_dp*t382*t383*t400 - 0.40e1_dp &
1437 : *t11*t255 + 0.40e1_dp*t382*t230*t254*t168 - 0.20e1_dp* &
1438 : t6*t77*(0.1866666667e2_dp*beta*t237*t12 + 0.9866666667e2_dp* &
1439 : t11*t241 - 0.8266666668e2_dp*t393*t251 + 0.3200000001e2_dp*beta &
1440 0 : *t244*ndrho/t657/t7*t669)
1441 0 : t961 = t687*t26
1442 0 : t1002 = t3*t196
1443 : t1009 = 0.2e1_dp*t101*(-t455 + 0.9000000000e1_dp*t457 + 0.6000000000e1_dp &
1444 : *t460)*t103*t112 + 0.1800000000e2_dp*t412*t264*t40*t127 &
1445 : *t8*t172*t103*t112 + 0.2e1_dp*t101*t265*t428 + 0.1800000000e2_dp &
1446 : *t413*t186*t285 + 0.2e1_dp*t101*t103*(-0.5625000000e1_dp &
1447 : *t961*t1*t212*t172 + 0.4500000000e1_dp*t417*t418*t404 &
1448 : + 0.1500000000e1_dp*t417*t49*t94*t172 - 0.1000000000e1_dp* &
1449 : t105*t109*t404 + 0.2250000000e1_dp*t417*t1*t258*t172 - 0.1500000000e1_dp &
1450 : *t105*t31*t933 + 0.3333333334e0_dp*t105*t281* &
1451 : t172) + 0.1250000000e0_dp*t722*t118*t860 - 0.8333333335e-1_dp*t289 &
1452 : *t435*t212 - 0.1666666667e0_dp*t289*t118*t867 + 0.5555555555e-1_dp &
1453 : *t289*t293*t368 - 0.1111111111e0_dp*t117*t1002*t94 &
1454 0 : - 0.1111111111e0_dp*t117*t293*t404
1455 0 : t1013 = t37*t492
1456 0 : t1044 = t769*pi
1457 : t1069 = 0.5400000000e2_dp*t1044*t8*t212*t172 - 0.3600000000e2_dp &
1458 : *t451*t452*t404 - 0.2400000000e2_dp*t451*t37*t94*t172 + &
1459 : 0.1200000000e2_dp*t128*t132*t404 - 0.1800000000e2_dp*t451*t8* &
1460 : t258*t172 + 0.8999999998e1_dp*t128*t43*t933 - 0.2000000000e1_dp &
1461 0 : *t128*t323*t172
1462 0 : t1091 = t127*t172*t46
1463 : t1102 = t1069*t46 + 0.8999999998e1_dp*t326*t40*t127*t467 + REAL(2 &
1464 : *t136*t462, KIND=dp) + 0.8999999998e1_dp*t328*t40*t127*t467 - &
1465 : 0.5555555558e-1_dp*t39*t933*t52 - 0.5000000001e0_dp*t258*t127 &
1466 : *t466 + 0.7407407410e-1_dp*t470*t143 + 0.6666666668e0_dp*t94*t478 &
1467 : *t1091 - 0.1111111112e0_dp*t470*t48*t148 - 0.1111111112e0_dp &
1468 0 : *t335*t486 - 0.1000000001e1_dp*t94*t135*t1091
1469 : t1136 = -0.6172839508e-1_dp*t190*t339 - 0.5555555556e0_dp*t41/t7 &
1470 : *t466 + 0.7407407410e-1_dp*t482*t343 + 0.7407407410e-1_dp*t146* &
1471 : t141*t462*t46 + 0.6666666668e0_dp*t479*t135*t172*t46 - 0.5555555558e-1_dp &
1472 : *t482*t347 - 0.5555555558e-1_dp*t146*t49*t1069 &
1473 : *t46 - 0.5000000001e0_dp*t41*t326*t466 - 0.5555555558e-1_dp*t482 &
1474 : *t351 - 0.1111111112e0_dp*t146*t147*t463 - 0.5000000001e0_dp* &
1475 0 : t41*t328*t466
1476 : t1141 = -0.1666666667e0_dp*t289*t297*t368 + 0.3333333334e0_dp*t117 &
1477 : *t1013*t94 + 0.3333333334e0_dp*t117*t297*t404 - 0.8333333335e-1_dp &
1478 : *t289*t118*t880 + 0.1666666667e0_dp*t117*t435*t258 + &
1479 : 0.1666666667e0_dp*t117*t118*t933 + 0.7407407405e-1_dp*t117*t735 &
1480 : *t172 + 0.1481481481e0_dp*t36*t304*t196 - 0.1111111111e0_dp* &
1481 : t117*t739*t172 - 0.2222222222e0_dp*t36*t122*t492 + 0.1666666667e0_dp &
1482 : *t117*t746*t172 + 0.3333333334e0_dp*t36*t38*(t1102 &
1483 0 : + t1136)
1484 : t1146 = -0.3333333336e0_dp*t117*t59*t94*t172 - 0.6666666672e0_dp &
1485 : *t36*t372 - 0.6666666672e0_dp*t36*t405 - 0.6666666672e0_dp*t36 &
1486 : *t408 - 0.4444444448e0_dp*t64*t857 + 0.8333333340e-1_dp*t605*t60 &
1487 : *t860 - 0.1666666668e0_dp*t211*t365*t212 - 0.3333333336e0_dp* &
1488 : t211*t60*t867 - 0.3333333336e0_dp*t211*t207*t368 - 0.6666666672e0_dp &
1489 : *t68*t857*t94 - 0.6666666672e0_dp*t68*t207*t404 - 0.1666666668e0_dp &
1490 : *t211*t60*t880 - 0.3333333336e0_dp*t68*t365* &
1491 : t258 - 0.3333333336e0_dp*t68*t60*t933 - 0.3333333336e0_dp*t68* &
1492 0 : t601*t172 - 0.2222222224e0_dp*t24*t98*(t1009 + t1141)
1493 :
1494 0 : e_ndrho_rho_rho_spin(ii) = e_ndrho_rho_rho_spin(ii) + t1146*sx
1495 :
1496 0 : t1153 = t27*t590
1497 0 : t1156 = t94*t501
1498 0 : t1163 = t404*t172
1499 0 : t1167 = t94*t529
1500 0 : t1177 = beta*t71
1501 : t1220 = -0.1066666667e2_dp*t1177*t17 + 0.2133333334e2_dp*t11*t377 &
1502 : - 0.1066666667e2_dp*t6*t632*t513 + 0.5333333333e1_dp*t6*t225 &
1503 : *t525 - 0.40e1_dp*t508*t76*t90 + 0.160e2_dp*t11*t10*t384 - &
1504 : 0.80e1_dp*t11*t401 - 0.120e2_dp*t382*t640*t90*t513 + 0.80e1_dp &
1505 : *t382*t230*t400*t168 + 0.40e1_dp*t382*t383*t525 - 0.20e1_dp &
1506 : *t6*t77*(-0.3200000000e2_dp*t1177*t86 + 0.4800000000e2_dp*t6 &
1507 0 : *t397 - 0.2400000000e2_dp*t245/t657/rho*t669)
1508 0 : t1284 = t37*t586
1509 : t1334 = 0.5400000000e2_dp*t1044*t452*t501 - 0.3600000000e2_dp*t451 &
1510 : *t8*t404*t172 - 0.1800000000e2_dp*t451*t452*t529 + 0.8999999998e1_dp &
1511 : *t128*t43*t1220 - 0.1200000000e2_dp*t313*t132*t501 &
1512 0 : + 0.5999999999e1_dp*t128*t132*t529
1513 0 : t1341 = t501*t46
1514 0 : t1345 = t529*t46
1515 0 : t1370 = t40*pi
1516 : t1396 = t1334*t46 + 0.1800000000e2_dp*t462*t40*t127*t467 - 0.2250000000e2_dp &
1517 : *t464*t312*t43*t1341 + 0.8999999998e1_dp*t465 &
1518 : *t43*t1345 - 0.1000000000e1_dp*t404*t127*t466 + 0.8099999996e2_dp &
1519 : *t135*t571*t573*t533*t2*t1341 - 0.5555555558e-1_dp*t39 &
1520 : *t1220*t52 - 0.5000000001e0_dp*t473*t1345 + 0.3333333334e0_dp* &
1521 : t479*t1345 + 0.1000000000e1_dp*t94*t312*t1341 - 0.4500000000e1_dp &
1522 : *t94*t573*t501*t1370*t8*t46 + 0.3703703705e-1_dp*t580 &
1523 : *t143 + 0.3000000000e1_dp*t312*t37*t501*t1370*t46 - 0.5555555558e-1_dp &
1524 : *t580*t48*t148 - 0.1111111112e0_dp*t482*t486 - 0.5555555558e-1_dp &
1525 : *t146*t49*t1334*t46 - 0.1000000000e1_dp*t41* &
1526 0 : t462*t466 - 0.5000000001e0_dp*t489*t1345
1527 : t1400 = -0.3600000000e2_dp*t412*t313*t562*t113 + 0.1800000000e2_dp &
1528 : *t413*t566*t113 + 0.3600000000e2_dp*t413*t186*t429 + 0.1620000000e3_dp &
1529 : *t26*t533*t100*t574*t576*t113 + 0.2e1_dp*t101 &
1530 : *t103*(-0.5625000000e1_dp*t961*t418*t501 + 0.4500000000e1_dp &
1531 : *t417*t1*t404*t172 + 0.2250000000e1_dp*t417*t418*t529 - &
1532 : 0.1500000000e1_dp*t105*t31*t1220 + 0.7500000000e0_dp*t271*t109 &
1533 : *t501 - 0.5000000000e0_dp*t105*t109*t529) + 0.1250000000e0_dp* &
1534 : t722*t118*t1156 - 0.1666666667e0_dp*t289*t435*t368 - 0.1666666667e0_dp &
1535 : *t289*t118*t1163 - 0.8333333335e-1_dp*t289*t118*t1167 &
1536 : + 0.1666666667e0_dp*t117*t1284*t94 + 0.3333333334e0_dp*t117 &
1537 : *t435*t404 + 0.1666666667e0_dp*t117*t118*t1220 + 0.2777777778e-1_dp &
1538 : *t289*t293*t501 - 0.1111111111e0_dp*t117*t1002*t172 &
1539 : - 0.5555555555e-1_dp*t117*t293*t529 - 0.1111111111e0_dp*t36*t122 &
1540 : *t586 - 0.8333333335e-1_dp*t289*t297*t501 + 0.3333333334e0_dp &
1541 : *t117*t1013*t172 + 0.1666666667e0_dp*t117*t297*t529 + 0.3333333334e0_dp &
1542 0 : *t36*t38*t1396
1543 : t1404 = -0.1666666668e0_dp*t116*t502 - 0.6666666672e0_dp*t36*t505 &
1544 : - 0.3333333336e0_dp*t36*t530 - 0.2222222224e0_dp*t64*t1153 + 0.8333333340e-1_dp &
1545 : *t605*t60*t1156 - 0.3333333336e0_dp*t211*t365 &
1546 : *t368 - 0.3333333336e0_dp*t211*t60*t1163 - 0.1666666668e0_dp*t211 &
1547 : *t60*t1167 - 0.3333333336e0_dp*t68*t1153*t94 - 0.6666666672e0_dp &
1548 : *t68*t365*t404 - 0.3333333336e0_dp*t68*t60*t1220 - 0.1666666668e0_dp &
1549 : *t211*t207*t501 - 0.6666666672e0_dp*t68*t857* &
1550 : t172 - 0.3333333336e0_dp*t68*t207*t529 - 0.2222222224e0_dp*t24* &
1551 0 : t98*t1400
1552 :
1553 0 : e_ndrho_ndrho_rho_spin(ii) = e_ndrho_ndrho_rho_spin(ii) + t1404*sx
1554 :
1555 0 : t1405 = t501*t172
1556 0 : t1412 = t172*t529
1557 : t1449 = -0.120e2_dp*t508*t76*t168 + 0.240e2_dp*t11*t514 - 0.120e2_dp &
1558 : *t11*t526 - 0.120e2_dp*t6*t641*t513*t168 + 0.120e2_dp*t382 &
1559 : *t230*t168*t525 - 0.20e1_dp*t6*t77*(-0.240e2_dp*beta*t521 &
1560 0 : *t250*ndrho + 0.180e2_dp*t393/t657*t669)
1561 0 : t1456 = t1405*t103
1562 0 : t1467 = t533*pi
1563 0 : t1472 = t572*t21
1564 : t1553 = 0.1350000000e3_dp*t537/t22/t572*rho*t1456 - 0.8100000000e2_dp &
1565 : *t534*t536*t539*rho*t172*t103*t529 - 0.2430000000e3_dp &
1566 : *t1467*t100/t570/omega/t22/t1472*t140*t1456 - &
1567 : 0.1125000000e2_dp*t177*t687*t1*t1405 + 0.1350000000e2_dp*t176 &
1568 : *t103*t28*t270*t1*t1412 - 0.3000000000e1_dp*t177*t105* &
1569 : t1*t1449 + 0.1250000000e0_dp*t722*t118*t1405 - 0.2500000000e0_dp &
1570 : *t289*t435*t501 - 0.2500000000e0_dp*t289*t118*t1412 + 0.5000000001e0_dp &
1571 : *t117*t1284*t172 + 0.5000000001e0_dp*t117*t435 &
1572 : *t529 + 0.1666666667e0_dp*t117*t118*t1449 + 0.3333333334e0_dp*t36 &
1573 : *t38*(0.6750000000e2_dp*t1044*t8*t1405*t46 - 0.6750000000e2_dp &
1574 : *t451*t186*t1345 - 0.5264999998e3_dp*t571/t1472*t533 &
1575 : *t2*t1405*t46 + 0.8999999998e1_dp*t185*t8*t1449*t46 + 0.2429999999e3_dp &
1576 : *t575*t2*t529*t466 + 0.7289999995e3_dp/t570/t39 &
1577 : /t572/t126*t1467*t7*t1405*t46 - 0.5555555558e-1_dp*t39 &
1578 0 : *t1449*t52 - 0.5000000001e0_dp*t41*t1449*t46)
1579 :
1580 : e_ndrho_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_ndrho_spin(ii) + (0.8333333340e-1_dp*t605*t60*t1405 - &
1581 : 0.5000000004e0_dp*t211*t365*t501 - 0.5000000004e0_dp*t211*t60*t1412 - 0.1000000001e1_dp &
1582 : *t68*t1153*t172 - 0.1000000001e1_dp*t68*t365*t529 - 0.3333333336e0_dp &
1583 : *t68*t60*t1449 - 0.2222222224e0_dp*t24*t98*t1553) &
1584 0 : *sx
1585 :
1586 : END IF
1587 : END IF
1588 : END DO
1589 : !$OMP END DO
1590 :
1591 104 : END SUBROUTINE xb88_lr_lsd_calc
1592 :
1593 : END MODULE xc_xbecke88_long_range
|
