(git:e7e05ae)
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-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  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 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  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  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
134  e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
135  e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho
136  TYPE(xc_derivative_type), POINTER :: deriv
137 
138  CALL timeset(routinen, handle)
139 
140  CALL section_vals_val_get(xb88_lr_params, "scale_x", r_val=sx)
141  CALL section_vals_val_get(xb88_lr_params, "omega", r_val=omega)
142 
143  CALL cite_reference(becke1988)
144 
145  CALL xc_rho_set_get(rho_set, rho=rho, &
146  norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
147  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  dummy => rho
152 
153  e_0 => dummy
154  e_rho => dummy
155  e_ndrho => dummy
156  e_rho_rho => dummy
157  e_ndrho_rho => dummy
158  e_ndrho_ndrho => dummy
159  e_rho_rho_rho => dummy
160  e_ndrho_rho_rho => dummy
161  e_ndrho_ndrho_rho => dummy
162  e_ndrho_ndrho_ndrho => dummy
163 
164  IF (grad_deriv >= 0) THEN
165  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
166  allocate_deriv=.true.)
167  CALL xc_derivative_get(deriv, deriv_data=e_0)
168  END IF
169  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
170  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
171  allocate_deriv=.true.)
172  CALL xc_derivative_get(deriv, deriv_data=e_rho)
173  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
174  allocate_deriv=.true.)
175  CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
176  END IF
177  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
178  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
179  allocate_deriv=.true.)
180  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
181  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
182  allocate_deriv=.true.)
183  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
184  deriv => xc_dset_get_derivative(deriv_set, &
185  [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
186  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
187  END IF
188  IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
189  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
190  allocate_deriv=.true.)
191  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
192  deriv => xc_dset_get_derivative(deriv_set, &
193  [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.true.)
194  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
195  deriv => xc_dset_get_derivative(deriv_set, &
196  [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.true.)
197  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
198  deriv => xc_dset_get_derivative(deriv_set, &
199  [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
200  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
201  END IF
202  IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
203  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 !$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  CALL timestop(handle)
225  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  SUBROUTINE xb88_lr_lda_calc(rho, norm_drho, &
255  e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
256  e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
257  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  my_epsilon_rho = epsilon_rho
291  epsilon_rho43 = my_epsilon_rho**(4.0_dp/3.0_dp)
292  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  my_rho = rho(ii)*0.5_dp
298  my_ndrho = norm_drho(ii)*0.5_dp
299  IF (my_rho > my_epsilon_rho) THEN
300  IF (grad_deriv >= 0) THEN
301  t1 = my_rho**(0.1e1_dp/0.3e1_dp)
302  t2 = t1*my_rho
303  t3 = 0.1e1_dp/t2
304  xx = my_ndrho*max(t3, epsilon_rho43)
305  t5 = my_ndrho**2
306  t6 = beta*t5
307  t7 = my_rho**2
308  t8 = t1**2
309  t10 = 0.1e1_dp/t8/t7
310  t11 = beta*my_ndrho
311  t12 = log(xx + sqrt(xx**0.2e1_dp + 0.1e1_dp))
312  t16 = 0.10e1_dp + 0.60e1_dp*t11*t3*t12
313  t17 = 0.1e1_dp/t16
314  t18 = t10*t17
315  t21 = 0.20e1_dp*cx + 0.20e1_dp*t6*t18
316  t22 = sqrt(t21)
317  t23 = t22*t21
318  t24 = my_rho*t23
319  t26 = rootpi
320  t27 = 0.1e1_dp/t26
321  t28 = 0.1e1_dp/omega
322  t29 = 0.1e1_dp/t22
323  t30 = t28*t29
324  t31 = t26*t1
325  t34 = erf(0.3000000000e1_dp*t30*t31)
326  t36 = omega*t22
327  t37 = 0.1e1_dp/t1
328  t38 = t27*t37
329  t39 = omega**2
330  t40 = 0.1e1_dp/t39
331  t41 = 0.1e1_dp/t21
332  t42 = t40*t41
333  t43 = pi*t8
334  t44 = t42*t43
335  t46 = exp(-0.8999999998e1_dp*t44)
336  t47 = t39*t21
337  t48 = 0.1e1_dp/pi
338  t49 = 0.1e1_dp/t8
339  t51 = t46 - 0.10e1_dp
340  t52 = t48*t49*t51
341  t55 = t46 - 0.15e1_dp - 0.5555555558e-1_dp*t47*t52
342  t59 = t26*t34 + 0.3333333334e0_dp*t36*t38*t55
343  t60 = t27*t59
344 
345  !! Multiply with 2.0 because it code comes from LSD
346  e_0(ii) = e_0(ii) - 0.2222222224e0_dp*t24*omega*t60*sx*2.0_dp
347 
348  END IF
349  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
350  t64 = t23*omega
351  t68 = my_rho*t22*omega
352  t69 = t7*my_rho
353  t71 = 0.1e1_dp/t8/t69
354  t72 = t71*t17
355  t75 = t16**2
356  t76 = 0.1e1_dp/t75
357  t77 = t10*t76
358  t79 = 0.1e1_dp/t1/t7
359  t84 = 1 + t5*t10
360  t85 = sqrt(t84)
361  t86 = 0.1e1_dp/t85
362  t87 = t71*t86
363  t90 = -0.8000000000e1_dp*t11*t79*t12 - 0.8000000000e1_dp*t6*t87
364  t91 = t77*t90
365  t94 = -0.5333333333e1_dp*t6*t72 - 0.20e1_dp*t6*t91
366  t95 = t60*t94
367  t98 = omega*t27
368  t99 = sqrt(0.3141592654e1_dp)
369  t100 = 0.1e1_dp/t99
370  t101 = t26*t100
371  t103 = exp(-0.9000000000e1_dp*t44)
372  t104 = 0.1e1_dp/t23
373  t105 = t28*t104
374  t109 = t26*t49
375  t112 = -0.1500000000e1_dp*t105*t31*t94 + 0.1000000000e1_dp*t30* &
376  t109
377  t113 = t103*t112
378  t116 = omega*t29
379  t117 = t116*t27
380  t118 = t37*t55
381  t122 = t27*t3
382  t126 = t21**2
383  t127 = 0.1e1_dp/t126
384  t128 = t40*t127
385  t130 = t128*t43*t94
386  t132 = pi*t37
387  t133 = t42*t132
388  t135 = 0.8999999998e1_dp*t130 - 0.5999999999e1_dp*t133
389  t136 = t135*t46
390  t137 = t39*t94
391  t140 = t8*my_rho
392  t141 = 0.1e1_dp/t140
393  t143 = t48*t141*t51
394  t146 = t47*t48
395  t147 = t49*t135
396  t148 = t147*t46
397  t151 = t136 - 0.5555555558e-1_dp*t137*t52 + 0.3703703705e-1_dp*t47 &
398  *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  t151
402 
403  e_rho(ii) = e_rho(ii) + (-0.2222222224e0_dp*t64*t60 - 0.3333333336e0_dp*t68*t95 - &
404  0.2222222224e0_dp*t24*t98*t155)*sx
405 
406  t168 = 0.60e1_dp*beta*t3*t12 + 0.60e1_dp*t11*t10*t86
407  t169 = t77*t168
408  t172 = 0.40e1_dp*t11*t18 - 0.20e1_dp*t6*t169
409  t173 = t60*t172
410  t176 = pi*t100
411  t177 = t176*t103
412  t185 = t128*pi
413  t186 = t8*t172
414  t190 = t39*t172
415  t196 = 0.8999999998e1_dp*t185*t186*t46 - 0.5555555558e-1_dp*t190 &
416  *t52 - 0.5000000001e0_dp*t41*t172*t46
417  t200 = -0.3000000000e1_dp*t177*t105*t1*t172 + 0.1666666667e0_dp* &
418  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  *t200)*sx
422 
423  END IF
424  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
425  t207 = t27*t155
426  t211 = my_rho*t29*omega
427  t212 = t94**2
428  t213 = t60*t212
429  t216 = t207*t94
430  t219 = t7**2
431  t221 = 0.1e1_dp/t8/t219
432  t222 = t221*t17
433  t225 = t71*t76
434  t226 = t225*t90
435  t230 = 0.1e1_dp/t75/t16
436  t231 = t10*t230
437  t232 = t90**2
438  t233 = t231*t232
439  t237 = 0.1e1_dp/t1/t69
440  t241 = t221*t86
441  t244 = t5**2
442  t245 = beta*t244
443  t250 = 0.1e1_dp/t85/t84
444  t251 = 0.1e1_dp/t1/t219/t69*t250
445  t254 = 0.1866666667e2_dp*t11*t237*t12 + 0.4000000000e2_dp*t6*t241 &
446  - 0.1066666667e2_dp*t245*t251
447  t255 = t77*t254
448  t258 = 0.1955555555e2_dp*t6*t222 + 0.1066666667e2_dp*t6*t226 + 0.40e1_dp &
449  *t6*t233 - 0.20e1_dp*t6*t255
450  t259 = t60*t258
451  t264 = 0.9000000000e1_dp*t130 - 0.6000000000e1_dp*t133
452  t265 = t264*t103
453  t270 = 0.1e1_dp/t22/t126
454  t271 = t28*t270
455  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  *t30*t281
459  t285 = t103*t284
460  t289 = omega*t104*t27
461  t293 = t3*t55
462  t297 = t37*t151
463  t304 = t27*t79
464  t311 = t126*t21
465  t312 = 0.1e1_dp/t311
466  t313 = t40*t312
467  t316 = 0.1800000000e2_dp*t313*t43*t212
468  t319 = 0.1200000000e2_dp*t128*t132*t94
469  t321 = t128*t43*t258
470  t323 = pi*t3
471  t325 = 0.2000000000e1_dp*t42*t323
472  t326 = -t316 + t319 + 0.8999999998e1_dp*t321 + t325
473  t328 = t135**2
474  t330 = t39*t258
475  t335 = t137*t48
476  t339 = t48*t10*t51
477  t343 = t141*t135*t46
478  t346 = t49*t326
479  t347 = t346*t46
480  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  *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  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  *t68*t259 - 0.2222222224e0_dp*t24*t98*t358)*sx
494 
495  t365 = t27*t200
496  t368 = t94*t172
497  t372 = t365*t94
498  t377 = t225*t168
499  t382 = t6*t10
500  t383 = t230*t90
501  t384 = t383*t168
502  t393 = beta*t5*my_ndrho
503  t397 = 0.1e1_dp/t1/t219/t7*t250
504  t400 = -0.8000000000e1_dp*beta*t79*t12 - 0.2400000000e2_dp*t11* &
505  t87 + 0.8000000000e1_dp*t393*t397
506  t401 = t77*t400
507  t404 = -0.1066666667e2_dp*t11*t72 + 0.5333333333e1_dp*t6*t377 - 0.40e1_dp &
508  *t11*t91 + 0.40e1_dp*t382*t384 - 0.20e1_dp*t6*t401
509  t405 = t60*t404
510  t408 = t207*t172
511  t412 = t26*pi*t100
512  t413 = t412*t128
513  t417 = t271*t26
514  t418 = t1*t94
515  t428 = 0.2250000000e1_dp*t417*t418*t172 - 0.1500000000e1_dp*t105 &
516  *t31*t404 - 0.5000000000e0_dp*t105*t109*t172
517  t429 = t103*t428
518  t435 = t37*t196
519  t451 = t313*pi
520  t452 = t8*t94
521  t455 = 0.1800000000e2_dp*t451*t452*t172
522  t457 = t128*t43*t404
523  t460 = t128*t132*t172
524  t462 = -t455 + 0.8999999998e1_dp*t457 + 0.5999999999e1_dp*t460
525  t463 = t462*t46
526  t464 = t135*t40
527  t465 = t464*t127
528  t466 = t172*t46
529  t467 = t43*t466
530  t470 = t39*t404
531  t473 = t94*t127
532  t478 = 0.1e1_dp/my_rho
533  t479 = t41*t478
534  t482 = t190*t48
535  t486 = t49*t462*t46
536  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  - 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  *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  *t24*t98*t496)*sx
551 
552  t501 = t172**2
553  t502 = t60*t501
554  t505 = t365*t172
555  t508 = beta*t10
556  t513 = t168**2
557  t514 = t231*t513
558  t519 = t219*my_rho
559  t521 = 0.1e1_dp/t1/t519
560  t525 = 0.120e2_dp*t508*t86 - 0.60e1_dp*t6*t521*t250
561  t526 = t77*t525
562  t529 = 0.40e1_dp*t508*t17 - 0.80e1_dp*t11*t169 + 0.40e1_dp*t6*t514 &
563  - 0.20e1_dp*t6*t526
564  t530 = t60*t529
565  t533 = pi**2
566  t534 = t533*t100
567  t536 = 0.1e1_dp/t39/omega
568  t537 = t534*t536
569  t539 = 0.1e1_dp/t22/t311
570  t562 = t8*t501
571  t566 = t8*t529
572  t570 = t39**2
573  t571 = 0.1e1_dp/t570
574  t572 = t126**2
575  t573 = 0.1e1_dp/t572
576  t574 = t571*t573
577  t575 = t574*t533
578  t576 = t2*t501
579  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  *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  *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  *sx
592 
593  END IF
594  IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
595  t601 = t27*t358
596  t605 = my_rho*t104*omega
597  t606 = t212*t94
598  t613 = t94*t258
599  t624 = 0.1e1_dp/t8/t519
600  t628 = t221*t76
601  t632 = t71*t230
602  t639 = t75**2
603  t640 = 0.1e1_dp/t639
604  t641 = t10*t640
605  t657 = t219**2
606  t667 = t84**2
607  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  /t69*t669)
615  t687 = t28*t539
616  t716 = t264**2
617  t722 = omega*t270*t27
618  t735 = t79*t55
619  t739 = t3*t151
620  t746 = t37*t354
621  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  *t128*t43*t677 - 0.2666666667e1_dp*t42*pi*t79
626  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  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.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  t68*t60*t677 - 0.2222222224e0_dp*t24*t98*t842
660 
661  e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + t846*sx
662 
663  t857 = t27*t496
664  t860 = t212*t172
665  t867 = t94*t404
666  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  *t244*my_ndrho/t657/t7*t669)
675  t961 = t687*t26
676  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.1111111111e0_dp*t117*t293*t404
689  t1013 = t37*t492
690  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  *t128*t323*t172
696  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  *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  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  + 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  t601*t172 - 0.2222222224e0_dp*t24*t98*(t1009 + t1141)
727 
728  e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) + t1146*sx
729 
730  t1153 = t27*t590
731  t1156 = t94*t501
732  t1163 = t404*t172
733  t1167 = t94*t529
734  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  *t397 - 0.2400000000e2_dp*t245/t657/my_rho*t669)
742  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.5999999999e1_dp*t128*t132*t529
747  t1341 = t501*t46
748  t1345 = t529*t46
749  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  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  *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  t98*t1400
786 
787  e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) + t1404*sx
788 
789  t1405 = t501*t172
790  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  *t250*my_ndrho + 0.180e2_dp*t393/t657*t669)
795  t1456 = t1405*t103
796  t1467 = t533*pi
797  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  *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  *sx
819 
820  END IF
821  END IF
822  END DO
823 !$OMP END DO
824 
825  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  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  POINTER :: dummy, e_0
850  TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
851  e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, e_rho_rho_rho, &
852  norm_drho, rho
853  TYPE(xc_derivative_type), POINTER :: deriv
854 
855  CALL timeset(routinen, handle)
856 
857  CALL cite_reference(becke1988)
858 
859  NULLIFY (deriv)
860  DO i = 1, 2
861  NULLIFY (norm_drho(i)%array, rho(i)%array)
862  END DO
863 
864  CALL section_vals_val_get(xb88_lr_params, "scale_x", r_val=sx)
865  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  local_bounds=bo)
872  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  dummy => rho(1)%array
877 
878  e_0 => dummy
879  DO i = 1, 2
880  e_rho(i)%array => dummy
881  e_ndrho(i)%array => dummy
882  e_rho_rho(i)%array => dummy
883  e_ndrho_rho(i)%array => dummy
884  e_ndrho_ndrho(i)%array => dummy
885  e_rho_rho_rho(i)%array => dummy
886  e_ndrho_rho_rho(i)%array => dummy
887  e_ndrho_ndrho_rho(i)%array => dummy
888  e_ndrho_ndrho_ndrho(i)%array => dummy
889  END DO
890 
891  IF (grad_deriv >= 0) THEN
892  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
893  allocate_deriv=.true.)
894  CALL xc_derivative_get(deriv, deriv_data=e_0)
895  END IF
896  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
897  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
898  allocate_deriv=.true.)
899  CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
900  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
901  allocate_deriv=.true.)
902  CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
903  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
904  allocate_deriv=.true.)
905  CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
906  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
907  allocate_deriv=.true.)
908  CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
909  END IF
910  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
911  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
912  allocate_deriv=.true.)
913  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  allocate_deriv=.true.)
916  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  allocate_deriv=.true.)
919  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  allocate_deriv=.true.)
922  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(2)%array)
923  deriv => xc_dset_get_derivative(deriv_set, &
924  [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.true.)
925  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(1)%array)
926  deriv => xc_dset_get_derivative(deriv_set, &
927  [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.true.)
928  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(2)%array)
929  END IF
930  IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
931  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
932  allocate_deriv=.true.)
933  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  allocate_deriv=.true.)
936  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(2)%array)
937  deriv => xc_dset_get_derivative(deriv_set, &
938  [deriv_norm_drhoa, deriv_rhoa, deriv_rhoa], allocate_deriv=.true.)
939  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(1)%array)
940  deriv => xc_dset_get_derivative(deriv_set, &
941  [deriv_norm_drhob, deriv_rhob, deriv_rhob], allocate_deriv=.true.)
942  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(2)%array)
943  deriv => xc_dset_get_derivative(deriv_set, &
944  [deriv_norm_drhoa, deriv_norm_drhoa, deriv_rhoa], allocate_deriv=.true.)
945  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(1)%array)
946  deriv => xc_dset_get_derivative(deriv_set, &
947  [deriv_norm_drhob, deriv_norm_drhob, deriv_rhob], allocate_deriv=.true.)
948  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(2)%array)
949  deriv => xc_dset_get_derivative(deriv_set, &
950  [deriv_norm_drhoa, deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.true.)
951  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(1)%array)
952  deriv => xc_dset_get_derivative(deriv_set, &
953  [deriv_norm_drhob, deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.true.)
954  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(2)%array)
955  END IF
956  IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
957  cpabort("derivatives bigger than 3 not implemented")
958  END IF
959 
960  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 !$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  CALL timestop(handle)
992 
993  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  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  my_epsilon_rho = 0.5_dp*epsilon_rho
1059  epsilon_rho43 = my_epsilon_rho**(4.0_dp/3.0_dp)
1060  cx = 1.5_dp*(3.0_dp/4.0_dp/pi)**(1.0_dp/3.0_dp)
1061 
1062 !$OMP DO
1063  DO ii = 1, npoints
1064  rho = rho_spin(ii)
1065  ndrho = norm_drho_spin(ii)
1066  IF (rho > my_epsilon_rho) THEN
1067  IF (grad_deriv >= 0) THEN
1068  t1 = rho**(0.1e1_dp/0.3e1_dp)
1069  t2 = t1*rho
1070  t3 = 0.1e1_dp/t2
1071  xx = ndrho*max(t3, epsilon_rho43)
1072  t5 = ndrho**2
1073  t6 = beta*t5
1074  t7 = rho**2
1075  t8 = t1**2
1076  t10 = 0.1e1_dp/t8/t7
1077  t11 = beta*ndrho
1078  t12 = log(xx + sqrt(xx**0.2e1_dp + 0.1e1_dp))
1079  t16 = 0.10e1_dp + 0.60e1_dp*t11*t3*t12
1080  t17 = 0.1e1_dp/t16
1081  t18 = t10*t17
1082  t21 = 0.20e1_dp*cx + 0.20e1_dp*t6*t18
1083  t22 = sqrt(t21)
1084  t23 = t22*t21
1085  t24 = rho*t23
1086  t26 = rootpi
1087  t27 = 0.1e1_dp/t26
1088  t28 = 0.1e1_dp/omega
1089  t29 = 0.1e1_dp/t22
1090  t30 = t28*t29
1091  t31 = t26*t1
1092  t34 = erf(0.3000000000e1_dp*t30*t31)
1093  t36 = omega*t22
1094  t37 = 0.1e1_dp/t1
1095  t38 = t27*t37
1096  t39 = omega**2
1097  t40 = 0.1e1_dp/t39
1098  t41 = 0.1e1_dp/t21
1099  t42 = t40*t41
1100  t43 = pi*t8
1101  t44 = t42*t43
1102  t46 = exp(-0.8999999998e1_dp*t44)
1103  t47 = t39*t21
1104  t48 = 0.1e1_dp/pi
1105  t49 = 0.1e1_dp/t8
1106  t51 = t46 - 0.10e1_dp
1107  t52 = t48*t49*t51
1108  t55 = t46 - 0.15e1_dp - 0.5555555558e-1_dp*t47*t52
1109  t59 = t26*t34 + 0.3333333334e0_dp*t36*t38*t55
1110  t60 = t27*t59
1111 
1112  e_0(ii) = e_0(ii) - 0.2222222224e0_dp*t24*omega*t60*sx
1113 
1114  END IF
1115  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1116  t64 = t23*omega
1117  t68 = rho*t22*omega
1118  t69 = t7*rho
1119  t71 = 0.1e1_dp/t8/t69
1120  t72 = t71*t17
1121  t75 = t16**2
1122  t76 = 0.1e1_dp/t75
1123  t77 = t10*t76
1124  t79 = 0.1e1_dp/t1/t7
1125  t84 = 1 + t5*t10
1126  t85 = sqrt(t84)
1127  t86 = 0.1e1_dp/t85
1128  t87 = t71*t86
1129  t90 = -0.8000000000e1_dp*t11*t79*t12 - 0.8000000000e1_dp*t6*t87
1130  t91 = t77*t90
1131  t94 = -0.5333333333e1_dp*t6*t72 - 0.20e1_dp*t6*t91
1132  t95 = t60*t94
1133  t98 = omega*t27
1134  t99 = sqrt(0.3141592654e1_dp)
1135  t100 = 0.1e1_dp/t99
1136  t101 = t26*t100
1137  t103 = exp(-0.9000000000e1_dp*t44)
1138  t104 = 0.1e1_dp/t23
1139  t105 = t28*t104
1140  t109 = t26*t49
1141  t112 = -0.1500000000e1_dp*t105*t31*t94 + 0.1000000000e1_dp*t30* &
1142  t109
1143  t113 = t103*t112
1144  t116 = omega*t29
1145  t117 = t116*t27
1146  t118 = t37*t55
1147  t122 = t27*t3
1148  t126 = t21**2
1149  t127 = 0.1e1_dp/t126
1150  t128 = t40*t127
1151  t130 = t128*t43*t94
1152  t132 = pi*t37
1153  t133 = t42*t132
1154  t135 = 0.8999999998e1_dp*t130 - 0.5999999999e1_dp*t133
1155  t136 = t135*t46
1156  t137 = t39*t94
1157  t140 = t8*rho
1158  t141 = 0.1e1_dp/t140
1159  t143 = t48*t141*t51
1160  t146 = t47*t48
1161  t147 = t49*t135
1162  t148 = t147*t46
1163  t151 = t136 - 0.5555555558e-1_dp*t137*t52 + 0.3703703705e-1_dp*t47 &
1164  *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  t151
1168 
1169  e_rho_spin(ii) = e_rho_spin(ii) + (-0.2222222224e0_dp*t64*t60 - 0.3333333336e0_dp*t68*t95 - &
1170  0.2222222224e0_dp*t24*t98*t155)*sx
1171 
1172  t168 = 0.60e1_dp*beta*t3*t12 + 0.60e1_dp*t11*t10*t86
1173  t169 = t77*t168
1174  t172 = 0.40e1_dp*t11*t18 - 0.20e1_dp*t6*t169
1175  t173 = t60*t172
1176  t176 = pi*t100
1177  t177 = t176*t103
1178  t185 = t128*pi
1179  t186 = t8*t172
1180  t190 = t39*t172
1181  t196 = 0.8999999998e1_dp*t185*t186*t46 - 0.5555555558e-1_dp*t190 &
1182  *t52 - 0.5000000001e0_dp*t41*t172*t46
1183  t200 = -0.3000000000e1_dp*t177*t105*t1*t172 + 0.1666666667e0_dp* &
1184  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  *t200)*sx
1188 
1189  END IF
1190  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1191  t207 = t27*t155
1192  t211 = rho*t29*omega
1193  t212 = t94**2
1194  t213 = t60*t212
1195  t216 = t207*t94
1196  t219 = t7**2
1197  t221 = 0.1e1_dp/t8/t219
1198  t222 = t221*t17
1199  t225 = t71*t76
1200  t226 = t225*t90
1201  t230 = 0.1e1_dp/t75/t16
1202  t231 = t10*t230
1203  t232 = t90**2
1204  t233 = t231*t232
1205  t237 = 0.1e1_dp/t1/t69
1206  t241 = t221*t86
1207  t244 = t5**2
1208  t245 = beta*t244
1209  t250 = 0.1e1_dp/t85/t84
1210  t251 = 0.1e1_dp/t1/t219/t69*t250
1211  t254 = 0.1866666667e2_dp*t11*t237*t12 + 0.4000000000e2_dp*t6*t241 &
1212  - 0.1066666667e2_dp*t245*t251
1213  t255 = t77*t254
1214  t258 = 0.1955555555e2_dp*t6*t222 + 0.1066666667e2_dp*t6*t226 + 0.40e1_dp &
1215  *t6*t233 - 0.20e1_dp*t6*t255
1216  t259 = t60*t258
1217  t264 = 0.9000000000e1_dp*t130 - 0.6000000000e1_dp*t133
1218  t265 = t264*t103
1219  t270 = 0.1e1_dp/t22/t126
1220  t271 = t28*t270
1221  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  *t30*t281
1225  t285 = t103*t284
1226  t289 = omega*t104*t27
1227  t293 = t3*t55
1228  t297 = t37*t151
1229  t304 = t27*t79
1230  t311 = t126*t21
1231  t312 = 0.1e1_dp/t311
1232  t313 = t40*t312
1233  t316 = 0.1800000000e2_dp*t313*t43*t212
1234  t319 = 0.1200000000e2_dp*t128*t132*t94
1235  t321 = t128*t43*t258
1236  t323 = pi*t3
1237  t325 = 0.2000000000e1_dp*t42*t323
1238  t326 = -t316 + t319 + 0.8999999998e1_dp*t321 + t325
1239  t328 = t135**2
1240  t330 = t39*t258
1241  t335 = t137*t48
1242  t339 = t48*t10*t51
1243  t343 = t141*t135*t46
1244  t346 = t49*t326
1245  t347 = t346*t46
1246  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  *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  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  *t68*t259 - 0.2222222224e0_dp*t24*t98*t358)*sx
1260 
1261  t365 = t27*t200
1262  t368 = t94*t172
1263  t372 = t365*t94
1264  t377 = t225*t168
1265  t382 = t6*t10
1266  t383 = t230*t90
1267  t384 = t383*t168
1268  t393 = beta*t5*ndrho
1269  t397 = 0.1e1_dp/t1/t219/t7*t250
1270  t400 = -0.8000000000e1_dp*beta*t79*t12 - 0.2400000000e2_dp*t11* &
1271  t87 + 0.8000000000e1_dp*t393*t397
1272  t401 = t77*t400
1273  t404 = -0.1066666667e2_dp*t11*t72 + 0.5333333333e1_dp*t6*t377 - 0.40e1_dp &
1274  *t11*t91 + 0.40e1_dp*t382*t384 - 0.20e1_dp*t6*t401
1275  t405 = t60*t404
1276  t408 = t207*t172
1277  t412 = t26*pi*t100
1278  t413 = t412*t128
1279  t417 = t271*t26
1280  t418 = t1*t94
1281  t428 = 0.2250000000e1_dp*t417*t418*t172 - 0.1500000000e1_dp*t105 &
1282  *t31*t404 - 0.5000000000e0_dp*t105*t109*t172
1283  t429 = t103*t428
1284  t435 = t37*t196
1285  t451 = t313*pi
1286  t452 = t8*t94
1287  t455 = 0.1800000000e2_dp*t451*t452*t172
1288  t457 = t128*t43*t404
1289  t460 = t128*t132*t172
1290  t462 = -t455 + 0.8999999998e1_dp*t457 + 0.5999999999e1_dp*t460
1291  t463 = t462*t46
1292  t464 = t135*t40
1293  t465 = t464*t127
1294  t466 = t172*t46
1295  t467 = t43*t466
1296  t470 = t39*t404
1297  t473 = t94*t127
1298  t478 = 0.1e1_dp/rho
1299  t479 = t41*t478
1300  t482 = t190*t48
1301  t486 = t49*t462*t46
1302  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.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  *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  *t24*t98*t496)*sx
1317 
1318  t501 = t172**2
1319  t502 = t60*t501
1320  t505 = t365*t172
1321  t508 = beta*t10
1322  t513 = t168**2
1323  t514 = t231*t513
1324  t519 = t219*rho
1325  t521 = 0.1e1_dp/t1/t519
1326  t525 = 0.120e2_dp*t508*t86 - 0.60e1_dp*t6*t521*t250
1327  t526 = t77*t525
1328  t529 = 0.40e1_dp*t508*t17 - 0.80e1_dp*t11*t169 + 0.40e1_dp*t6*t514 &
1329  - 0.20e1_dp*t6*t526
1330  t530 = t60*t529
1331  t533 = pi**2
1332  t534 = t533*t100
1333  t536 = 0.1e1_dp/t39/omega
1334  t537 = t534*t536
1335  t539 = 0.1e1_dp/t22/t311
1336  t562 = t8*t501
1337  t566 = t8*t529
1338  t570 = t39**2
1339  t571 = 0.1e1_dp/t570
1340  t572 = t126**2
1341  t573 = 0.1e1_dp/t572
1342  t574 = t571*t573
1343  t575 = t574*t533
1344  t576 = t2*t501
1345  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  *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  *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  *sx
1358 
1359  END IF
1360  IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
1361  t601 = t27*t358
1362  t605 = rho*t104*omega
1363  t606 = t212*t94
1364  t613 = t94*t258
1365  t624 = 0.1e1_dp/t8/t519
1366  t628 = t221*t76
1367  t632 = t71*t230
1368  t639 = t75**2
1369  t640 = 0.1e1_dp/t639
1370  t641 = t10*t640
1371  t657 = t219**2
1372  t667 = t84**2
1373  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  /t69*t669)
1381  t687 = t28*t539
1382  t716 = t264**2
1383  t722 = omega*t270*t27
1384  t735 = t79*t55
1385  t739 = t3*t151
1386  t746 = t37*t354
1387  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  *t128*t43*t677 - 0.2666666667e1_dp*t42*pi*t79
1392  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  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.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  t68*t60*t677 - 0.2222222224e0_dp*t24*t98*t842
1426 
1427  e_rho_rho_rho_spin(ii) = e_rho_rho_rho_spin(ii) + t846*sx
1428 
1429  t857 = t27*t496
1430  t860 = t212*t172
1431  t867 = t94*t404
1432  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  *t244*ndrho/t657/t7*t669)
1441  t961 = t687*t26
1442  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.1111111111e0_dp*t117*t293*t404
1455  t1013 = t37*t492
1456  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  *t128*t323*t172
1462  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  *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  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  + 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  t601*t172 - 0.2222222224e0_dp*t24*t98*(t1009 + t1141)
1493 
1494  e_ndrho_rho_rho_spin(ii) = e_ndrho_rho_rho_spin(ii) + t1146*sx
1495 
1496  t1153 = t27*t590
1497  t1156 = t94*t501
1498  t1163 = t404*t172
1499  t1167 = t94*t529
1500  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  *t397 - 0.2400000000e2_dp*t245/t657/rho*t669)
1508  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.5999999999e1_dp*t128*t132*t529
1513  t1341 = t501*t46
1514  t1345 = t529*t46
1515  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  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  *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  t98*t1400
1552 
1553  e_ndrho_ndrho_rho_spin(ii) = e_ndrho_ndrho_rho_spin(ii) + t1404*sx
1554 
1555  t1405 = t501*t172
1556  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  *t250*ndrho + 0.180e2_dp*t393/t657*t669)
1561  t1456 = t1405*t103
1562  t1467 = t533*pi
1563  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  *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  *sx
1585 
1586  END IF
1587  END IF
1588  END DO
1589 !$OMP END DO
1590 
1591  END SUBROUTINE xb88_lr_lsd_calc
1592 
1593 END MODULE xc_xbecke88_long_range
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::becke1988
integer, save, public becke1988
Definition: bibliography.F:43
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:1040
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:66
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:123
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:838