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