(git:ccc2433)
xc_xbecke88.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 Becke 88 exchange functional
10 !> \par History
11 !> 11.2003 created [fawzi]
12 !> \author fawzi
13 ! **************************************************************************************************
14 MODULE xc_xbecke88
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  USE xc_derivative_desc, ONLY: deriv_norm_drho,&
23  deriv_norm_drhoa,&
24  deriv_norm_drhob,&
25  deriv_rho,&
26  deriv_rhoa,&
27  deriv_rhob
28  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
29  xc_dset_get_derivative
30  USE xc_derivative_types, ONLY: xc_derivative_get,&
31  xc_derivative_type
32  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
33  USE xc_rho_set_types, ONLY: xc_rho_set_get,&
34  xc_rho_set_type
35 #include "../base/base_uses.f90"
36 
37  IMPLICIT NONE
38  PRIVATE
39 
40  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
41  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbecke88'
42  REAL(kind=dp), PARAMETER :: beta = 0.0042_dp
43 
44  PUBLIC :: xb88_lda_info, xb88_lsd_info, xb88_lda_eval, xb88_lsd_eval
45 CONTAINS
46 
47 ! **************************************************************************************************
48 !> \brief return various information on the functional
49 !> \param reference string with the reference of the actual functional
50 !> \param shortform string with the shortform of the functional name
51 !> \param needs the components needed by this functional are set to
52 !> true (does not set the unneeded components to false)
53 !> \param max_deriv ...
54 !> \par History
55 !> 11.2003 created [fawzi]
56 !> \author fawzi
57 ! **************************************************************************************************
58  SUBROUTINE xb88_lda_info(reference, shortform, needs, max_deriv)
59  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
60  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
61  INTEGER, INTENT(out), OPTIONAL :: max_deriv
62 
63  IF (PRESENT(reference)) THEN
64  reference = "A. Becke, Phys. Rev. A 38, 3098 (1988) {LDA version}"
65  END IF
66  IF (PRESENT(shortform)) THEN
67  shortform = "Becke 1988 Exchange Functional (LDA)"
68  END IF
69  IF (PRESENT(needs)) THEN
70  needs%rho = .true.
71  needs%rho_1_3 = .true.
72  needs%norm_drho = .true.
73  END IF
74  IF (PRESENT(max_deriv)) max_deriv = 3
75 
76  END SUBROUTINE xb88_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 !> 11.2003 created [fawzi]
87 !> \author fawzi
88 ! **************************************************************************************************
89  SUBROUTINE xb88_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_lsd_info
108 
109 ! **************************************************************************************************
110 !> \brief evaluates the becke 88 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_params input parameters (scaling)
118 !> \par History
119 !> 11.2003 created [fawzi]
120 !> 01.2007 added scaling [Manuel Guidon]
121 !> \author fawzi
122 ! **************************************************************************************************
123  SUBROUTINE xb88_lda_eval(rho_set, deriv_set, grad_deriv, xb88_params)
124  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
125  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
126  INTEGER, INTENT(in) :: grad_deriv
127  TYPE(section_vals_type), POINTER :: xb88_params
128 
129  CHARACTER(len=*), PARAMETER :: routinen = 'xb88_lda_eval'
130 
131  INTEGER :: handle, npoints
132  INTEGER, DIMENSION(2, 3) :: bo
133  REAL(kind=dp) :: epsilon_rho, sx
134  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
135  e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
136  e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho, rho_1_3
137  TYPE(xc_derivative_type), POINTER :: deriv
138 
139  CALL timeset(routinen, handle)
140 
141  CALL section_vals_val_get(xb88_params, "scale_x", r_val=sx)
142 
143  CALL cite_reference(becke1988)
144 
145  CALL xc_rho_set_get(rho_set, rho_1_3=rho_1_3, 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  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  IF (grad_deriv >= 0) THEN
163  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
164  allocate_deriv=.true.)
165  CALL xc_derivative_get(deriv, deriv_data=e_0)
166  END IF
167  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
168  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
169  allocate_deriv=.true.)
170  CALL xc_derivative_get(deriv, deriv_data=e_rho)
171  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
172  allocate_deriv=.true.)
173  CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
174  END IF
175  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
176  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
177  allocate_deriv=.true.)
178  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
179  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
180  allocate_deriv=.true.)
181  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
182  deriv => xc_dset_get_derivative(deriv_set, &
183  [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
184  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
185  END IF
186  IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
187  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
188  allocate_deriv=.true.)
189  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
190  deriv => xc_dset_get_derivative(deriv_set, &
191  [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.true.)
192  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
193  deriv => xc_dset_get_derivative(deriv_set, &
194  [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.true.)
195  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
196  deriv => xc_dset_get_derivative(deriv_set, &
197  [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
198  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
199  END IF
200  IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
201  cpabort("derivatives bigger than 3 not implemented")
202  END IF
203 
204 !$OMP PARALLEL DEFAULT(NONE) &
205 !$OMP SHARED(rho, rho_1_3, norm_drho, e_0, e_rho) &
206 !$OMP SHARED(e_ndrho, e_rho_rho, e_ndrho_rho) &
207 !$OMP SHARED(e_ndrho_ndrho, e_rho_rho_rho) &
208 !$OMP SHARED(e_ndrho_rho_rho, e_ndrho_ndrho_rho) &
209 !$OMP SHARED(e_ndrho_ndrho_ndrho, grad_deriv, npoints) &
210 !$OMP SHARED(epsilon_rho, sx)
211 
212  CALL xb88_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
213  e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
214  e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
215  e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
216  e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, &
217  e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, grad_deriv=grad_deriv, &
218  npoints=npoints, epsilon_rho=epsilon_rho, sx=sx)
219 
220 !$OMP END PARALLEL
221 
222  CALL timestop(handle)
223  END SUBROUTINE xb88_lda_eval
224 
225 ! **************************************************************************************************
226 !> \brief evaluates the becke 88 exchange functional for lda
227 !> \param rho the density where you want to evaluate the functional
228 !> \param rho_1_3 ...
229 !> \param norm_drho ...
230 !> \param e_0 ...
231 !> \param e_rho ...
232 !> \param e_ndrho ...
233 !> \param e_rho_rho ...
234 !> \param e_ndrho_rho ...
235 !> \param e_ndrho_ndrho ...
236 !> \param e_rho_rho_rho ...
237 !> \param e_ndrho_rho_rho ...
238 !> \param e_ndrho_ndrho_rho ...
239 !> \param e_ndrho_ndrho_ndrho ...
240 !> \param grad_deriv degree of the derivative that should be evaluated,
241 !> if positive all the derivatives up to the given degree are evaluated,
242 !> if negative only the given degree is calculated
243 !> \param npoints ...
244 !> \param epsilon_rho ...
245 !> \param sx scaling-parameter for exchange
246 !> \par History
247 !> 11.2003 created [fawzi]
248 !> 01.2007 added scaling [Manuel Guidon]
249 !> \author fawzi
250 ! **************************************************************************************************
251  SUBROUTINE xb88_lda_calc(rho, rho_1_3, norm_drho, &
252  e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
253  e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
254  e_ndrho_ndrho_ndrho, grad_deriv, npoints, epsilon_rho, sx)
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_1_3, rho
260  REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
261 
262  INTEGER :: ii
263  REAL(kind=dp) :: c, epsilon_rho43, my_rho, my_rho_1_3, t1, t10, t100, t104, t11, t12, t126, &
264  t13, t159, t164, t17, t170, t176, t18, t189, t2, t21, t23, t29, t3, t31, t33, t37, t39, &
265  t40, t43, t44, t49, t5, t51, t54, t6, t64, t67, t7, t72, t74, t75, t79, t83, t86, t90, &
266  t91, t98, t99, x
267 
268  c = 1.5_dp*(0.75_dp/pi)**(1.0_dp/3.0_dp)
269  t1 = 2**(0.1e1_dp/0.3e1_dp)
270  t2 = t1**2
271  t5 = beta*t2
272  t7 = beta*t1
273  epsilon_rho43 = epsilon_rho**(4./3.)
274 
275  SELECT CASE (grad_deriv)
276  CASE (0)
277 
278 !$OMP DO
279 
280  DO ii = 1, npoints
281 
282  my_rho = rho(ii)
283  IF (my_rho > epsilon_rho) THEN
284  my_rho_1_3 = rho_1_3(ii)
285  t3 = my_rho_1_3*my_rho
286  x = norm_drho(ii)/max(t3, epsilon_rho43)
287  t6 = x**2
288  t10 = t2*t6 + 0.1e1_dp
289  t11 = sqrt(t10)
290  t12 = t1*x + t11
291  t13 = log(t12)
292  t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
293  t18 = 0.1e1_dp/t17
294  t21 = -c - t5*t6*t18
295 
296  e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
297  END IF
298  END DO
299 
300 !$OMP END DO
301 
302  CASE (1)
303 
304 !$OMP DO
305 
306  DO ii = 1, npoints
307 
308  my_rho = rho(ii)
309  IF (my_rho > epsilon_rho) THEN
310  my_rho_1_3 = rho_1_3(ii)
311  t3 = my_rho_1_3*my_rho
312  x = norm_drho(ii)/t3
313  t6 = x**2
314  t10 = t2*t6 + 0.1e1_dp
315  t11 = sqrt(t10)
316  t12 = t1*x + t11
317  t13 = log(t12)
318  t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
319  t18 = 0.1e1_dp/t17
320  t21 = -c - t5*t6*t18
321 
322  e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
323 
324  t23 = t2*my_rho_1_3
325  t29 = 0.1e1_dp/t11
326  t31 = t1 + t2*x*t29
327  t33 = 0.1e1_dp/t12
328  t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
329  t39 = t17**2
330  t40 = 0.1e1_dp/t39
331  t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
332  t44 = t43*x
333 
334  e_rho(ii) = e_rho(ii) &
335  - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
336  t23*t21)*sx
337  e_ndrho(ii) = e_ndrho(ii) &
338  + (t2*t43/0.2e1_dp)*sx
339  END IF
340 
341  END DO
342 
343 !$OMP END DO
344 
345  CASE (-1)
346 
347 !$OMP DO
348 
349  DO ii = 1, npoints
350 
351  my_rho = rho(ii)
352  IF (my_rho > epsilon_rho) THEN
353  my_rho_1_3 = rho_1_3(ii)
354  t3 = my_rho_1_3*my_rho
355  x = norm_drho(ii)/t3
356  t6 = x**2
357  t10 = t2*t6 + 0.1e1_dp
358  t11 = sqrt(t10)
359  t12 = t1*x + t11
360  t13 = log(t12)
361  t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
362  t18 = 0.1e1_dp/t17
363  t21 = -c - t5*t6*t18
364 
365  t23 = t2*my_rho_1_3
366  t29 = 0.1e1_dp/t11
367  t31 = t1 + t2*x*t29
368  t33 = 0.1e1_dp/t12
369  t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
370  t39 = t17**2
371  t40 = 0.1e1_dp/t39
372  t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
373  t44 = t43*x
374 
375  e_rho(ii) = e_rho(ii) &
376  - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
377  t23*t21)*sx
378  e_ndrho(ii) = e_ndrho(ii) &
379  + (t2*t43/0.2e1_dp)*sx
380  END IF
381  END DO
382 
383 !$OMP END DO
384 
385  CASE (2)
386 
387 !$OMP DO
388 
389  DO ii = 1, npoints
390 
391  my_rho = rho(ii)
392  IF (my_rho > epsilon_rho) THEN
393  my_rho_1_3 = rho_1_3(ii)
394  t3 = my_rho_1_3*my_rho
395  x = norm_drho(ii)/t3
396  t6 = x**2
397  t10 = t2*t6 + 0.1e1_dp
398  t11 = sqrt(t10)
399  t12 = t1*x + t11
400  t13 = log(t12)
401  t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
402  t18 = 0.1e1_dp/t17
403  t21 = -c - t5*t6*t18
404 
405  e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
406 
407  t23 = t2*my_rho_1_3
408  t29 = 0.1e1_dp/t11
409  t31 = t1 + t2*x*t29
410  t33 = 0.1e1_dp/t12
411  t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
412  t39 = t17**2
413  t40 = 0.1e1_dp/t39
414  t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
415  t44 = t43*x
416 
417  e_rho(ii) = e_rho(ii) &
418  - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
419  t23*t21)*sx
420  e_ndrho(ii) = e_ndrho(ii) &
421  + (t2*t43/0.2e1_dp)*sx
422 
423  t49 = my_rho_1_3**2
424  t51 = t2/t49
425  t54 = x*t40
426  t64 = 0.1e1_dp/t11/t10
427  t67 = t2*t29 - 0.2e1_dp*t1*t6*t64
428  t72 = t31**2
429  t74 = t12**2
430  t75 = 0.1e1_dp/t74
431  t79 = 0.12e2_dp*t7*t31*t33 + 0.6e1_dp*t7*x*t67*t33 &
432  - 0.6e1_dp*t7*x*t72*t75
433  t83 = t37**2
434  t86 = 0.1e1_dp/t39/t17
435  t90 = -0.2e1_dp*t5*t18 + 0.4e1_dp*t5*t54*t37 + t5*t6* &
436  t79*t40 - 0.2e1_dp*t5*t6*t83*t86
437  t91 = t90*t6
438  t98 = 0.1e1_dp/my_rho
439  t99 = t2*t98
440  t100 = t90*x
441  t104 = 0.1e1_dp/t3
442 
443  e_rho_rho(ii) = e_rho_rho(ii) &
444  + (0.8e1_dp/0.9e1_dp*t51*t91 - 0.2e1_dp/0.9e1_dp* &
445  t51*t44 + 0.2e1_dp/0.9e1_dp*t51*t21)*sx
446  e_ndrho_rho(ii) = e_ndrho_rho(ii) &
447  - (0.2e1_dp/0.3e1_dp*t99*t100)*sx
448  e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
449  + (t2*t90*t104/0.2e1_dp)*sx
450  END IF
451  END DO
452 
453 !$OMP END DO
454 
455  CASE (-2)
456 
457 !$OMP DO
458 
459  DO ii = 1, npoints
460 
461  my_rho = rho(ii)
462  IF (my_rho > epsilon_rho) THEN
463  my_rho_1_3 = rho_1_3(ii)
464  t3 = my_rho_1_3*my_rho
465  x = norm_drho(ii)/t3
466  t6 = x**2
467  t10 = t2*t6 + 0.1e1_dp
468  t11 = sqrt(t10)
469  t12 = t1*x + t11
470  t13 = log(t12)
471  t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
472  t18 = 0.1e1_dp/t17
473  t21 = -c - t5*t6*t18
474 
475  t23 = t2*my_rho_1_3
476  t29 = 0.1e1_dp/t11
477  t31 = t1 + t2*x*t29
478  t33 = 0.1e1_dp/t12
479  t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
480  t39 = t17**2
481  t40 = 0.1e1_dp/t39
482  t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
483  t44 = t43*x
484 
485  t49 = my_rho_1_3**2
486  t51 = t2/t49
487  t54 = x*t40
488  t64 = 0.1e1_dp/t11/t10
489  t67 = t2*t29 - 0.2e1_dp*t1*t6*t64
490  t72 = t31**2
491  t74 = t12**2
492  t75 = 0.1e1_dp/t74
493  t79 = 0.12e2_dp*t7*t31*t33 + 0.6e1_dp*t7*x*t67*t33 &
494  - 0.6e1_dp*t7*x*t72*t75
495  t83 = t37**2
496  t86 = 0.1e1_dp/t39/t17
497  t90 = -0.2e1_dp*t5*t18 + 0.4e1_dp*t5*t54*t37 + t5*t6* &
498  t79*t40 - 0.2e1_dp*t5*t6*t83*t86
499  t91 = t90*t6
500  t98 = 0.1e1_dp/my_rho
501  t99 = t2*t98
502  t100 = t90*x
503  t104 = 0.1e1_dp/t3
504 
505  e_rho_rho(ii) = e_rho_rho(ii) &
506  + (0.8e1_dp/0.9e1_dp*t51*t91 - 0.2e1_dp/0.9e1_dp* &
507  t51*t44 + 0.2e1_dp/0.9e1_dp*t51*t21)*sx
508  e_ndrho_rho(ii) = e_ndrho_rho(ii) &
509  - (0.2e1_dp/0.3e1_dp*t99*t100)*sx
510  e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
511  + (t2*t90*t104/0.2e1_dp)*sx
512  END IF
513  END DO
514 
515 !$OMP END DO
516 
517  CASE default
518 
519 !$OMP DO
520 
521  DO ii = 1, npoints
522 
523  my_rho = rho(ii)
524  IF (my_rho > epsilon_rho) THEN
525  my_rho_1_3 = rho_1_3(ii)
526  t3 = my_rho_1_3*my_rho
527  x = norm_drho(ii)/t3
528  t6 = x**2
529  t10 = t2*t6 + 0.1e1_dp
530  t11 = sqrt(t10)
531  t12 = t1*x + t11
532  t13 = log(t12)
533  t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
534  t18 = 0.1e1_dp/t17
535  t21 = -c - t5*t6*t18
536 
537  IF (grad_deriv >= 0) THEN
538  e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
539  END IF
540 
541  t23 = t2*my_rho_1_3
542  t29 = 0.1e1_dp/t11
543  t31 = t1 + t2*x*t29
544  t33 = 0.1e1_dp/t12
545  t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
546  t39 = t17**2
547  t40 = 0.1e1_dp/t39
548  t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
549  t44 = t43*x
550 
551  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
552  e_rho(ii) = e_rho(ii) &
553  - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
554  t23*t21)*sx
555  e_ndrho(ii) = e_ndrho(ii) &
556  + (t2*t43/0.2e1_dp)*sx
557  END IF
558 
559  t49 = my_rho_1_3**2
560  t51 = t2/t49
561  t54 = x*t40
562  t64 = 0.1e1_dp/t11/t10
563  t67 = t2*t29 - 0.2e1_dp*t1*t6*t64
564  t72 = t31**2
565  t74 = t12**2
566  t75 = 0.1e1_dp/t74
567  t79 = 0.12e2_dp*t7*t31*t33 + 0.6e1_dp*t7*x*t67*t33 &
568  - 0.6e1_dp*t7*x*t72*t75
569  t83 = t37**2
570  t86 = 0.1e1_dp/t39/t17
571  t90 = -0.2e1_dp*t5*t18 + 0.4e1_dp*t5*t54*t37 + t5*t6* &
572  t79*t40 - 0.2e1_dp*t5*t6*t83*t86
573  t91 = t90*t6
574  t98 = 0.1e1_dp/my_rho
575  t99 = t2*t98
576  t100 = t90*x
577  t104 = 0.1e1_dp/t3
578 
579  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
580  e_rho_rho(ii) = e_rho_rho(ii) &
581  + (0.8e1_dp/0.9e1_dp*t51*t91 - 0.2e1_dp/0.9e1_dp* &
582  t51*t44 + 0.2e1_dp/0.9e1_dp*t51*t21)*sx
583  e_ndrho_rho(ii) = e_ndrho_rho(ii) &
584  - (0.2e1_dp/0.3e1_dp*t99*t100)*sx
585  e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
586  + (t2*t90*t104/0.2e1_dp)*sx
587  END IF
588 
589  t126 = t10**2
590  t159 = t39**2
591  t164 = 0.6e1_dp*t5*t40*t37 - 0.12e2_dp*t5*x*t86*t83 &
592  + 0.6e1_dp*t5*t54*t79 + t5*t6*(0.18e2_dp*t7*t67 &
593  *t33 - 0.18e2_dp*t7*t72*t75 + 0.6e1_dp*t7*x* &
594  (-0.6e1_dp*t1*t64*x + 0.12e2_dp*t6*x/t11/t126) &
595  *t33 - 0.18e2_dp*t7*x*t67*t75*t31 + 0.12e2_dp*t7 &
596  *x*t72*t31/t74/t12)*t40 - 0.6e1_dp*t5*t6*t79 &
597  *t86*t37 + 0.6e1_dp*t5*t6*t83*t37/t159
598  t170 = 0.8e1_dp/0.9e1_dp*t51*t164*t6 + 0.14e2_dp/0.9e1_dp &
599  *t51*t100
600  t176 = t2/t49/my_rho
601  t189 = my_rho**2
602 
603  IF (grad_deriv == -3 .OR. grad_deriv >= 3) THEN
604  e_rho_rho_rho(ii) = e_rho_rho_rho(ii) &
605  - (0.4e1_dp/0.3e1_dp*t170*x*t98 + 0.16e2_dp/0.27e2_dp &
606  *t176*t91 - 0.4e1_dp/0.27e2_dp*t176*t44 + &
607  0.4e1_dp/0.27e2_dp*t176*t21)*sx
608  e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) &
609  + (t170*t104)*sx
610  e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) &
611  + ((-0.2e1_dp/0.3e1_dp*t99*t164*x - &
612  0.2e1_dp/0.3e1_dp*t99*t90)*t104)*sx
613  e_ndrho_ndrho_ndrho(ii) = e_ndrho_ndrho_ndrho(ii) &
614  + (t2*t164/t49/t189/0.2e1_dp)*sx
615  END IF
616  END IF
617  END DO
618 
619 !$OMP END DO
620 
621  END SELECT
622 
623  END SUBROUTINE xb88_lda_calc
624 
625 ! **************************************************************************************************
626 !> \brief evaluates the becke 88 exchange functional for lsd
627 !> \param rho_set ...
628 !> \param deriv_set ...
629 !> \param grad_deriv ...
630 !> \param xb88_params ...
631 !> \par History
632 !> 11.2003 created [fawzi]
633 !> 01.2007 added scaling [Manuel Guidon]
634 !> \author fawzi
635 ! **************************************************************************************************
636  SUBROUTINE xb88_lsd_eval(rho_set, deriv_set, grad_deriv, xb88_params)
637  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
638  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
639  INTEGER, INTENT(in) :: grad_deriv
640  TYPE(section_vals_type), POINTER :: xb88_params
641 
642  CHARACTER(len=*), PARAMETER :: routinen = 'xb88_lsd_eval'
643 
644  INTEGER :: handle, i, ispin, npoints
645  INTEGER, DIMENSION(2, 3) :: bo
646  REAL(kind=dp) :: epsilon_rho, sx
647  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
648  POINTER :: dummy, e_0
649  TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
650  e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, e_rho_rho_rho, &
651  norm_drho, rho, rho_1_3
652  TYPE(xc_derivative_type), POINTER :: deriv
653 
654  CALL timeset(routinen, handle)
655 
656  CALL cite_reference(becke1988)
657 
658  NULLIFY (deriv)
659  DO i = 1, 2
660  NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
661  END DO
662 
663  CALL section_vals_val_get(xb88_params, "scale_x", r_val=sx)
664  CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
665  rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
666  rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
667  norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
668  local_bounds=bo)
669  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
670 
671  dummy => rho(1)%array
672 
673  e_0 => dummy
674  DO i = 1, 2
675  e_rho(i)%array => dummy
676  e_ndrho(i)%array => dummy
677  e_rho_rho(i)%array => dummy
678  e_ndrho_rho(i)%array => dummy
679  e_ndrho_ndrho(i)%array => dummy
680  e_rho_rho_rho(i)%array => dummy
681  e_ndrho_rho_rho(i)%array => dummy
682  e_ndrho_ndrho_rho(i)%array => dummy
683  e_ndrho_ndrho_ndrho(i)%array => dummy
684  END DO
685 
686  IF (grad_deriv >= 0) THEN
687  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
688  allocate_deriv=.true.)
689  CALL xc_derivative_get(deriv, deriv_data=e_0)
690  END IF
691  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
692  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
693  allocate_deriv=.true.)
694  CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
695  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
696  allocate_deriv=.true.)
697  CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
698  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
699  allocate_deriv=.true.)
700  CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
701  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
702  allocate_deriv=.true.)
703  CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
704  END IF
705  IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
706  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
707  allocate_deriv=.true.)
708  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(1)%array)
709  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
710  allocate_deriv=.true.)
711  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(2)%array)
712  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
713  allocate_deriv=.true.)
714  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(1)%array)
715  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
716  allocate_deriv=.true.)
717  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(2)%array)
718  deriv => xc_dset_get_derivative(deriv_set, &
719  [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.true.)
720  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(1)%array)
721  deriv => xc_dset_get_derivative(deriv_set, &
722  [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.true.)
723  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(2)%array)
724  END IF
725  IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
726  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
727  allocate_deriv=.true.)
728  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(1)%array)
729  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
730  allocate_deriv=.true.)
731  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(2)%array)
732  deriv => xc_dset_get_derivative(deriv_set, &
733  [deriv_norm_drhoa, deriv_rhoa, deriv_rhoa], allocate_deriv=.true.)
734  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(1)%array)
735  deriv => xc_dset_get_derivative(deriv_set, &
736  [deriv_norm_drhob, deriv_rhob, deriv_rhob], allocate_deriv=.true.)
737  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(2)%array)
738  deriv => xc_dset_get_derivative(deriv_set, &
739  [deriv_norm_drhoa, deriv_norm_drhoa, deriv_rhoa], allocate_deriv=.true.)
740  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(1)%array)
741  deriv => xc_dset_get_derivative(deriv_set, &
742  [deriv_norm_drhob, deriv_norm_drhob, deriv_rhob], allocate_deriv=.true.)
743  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(2)%array)
744  deriv => xc_dset_get_derivative(deriv_set, &
745  [deriv_norm_drhoa, deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.true.)
746  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(1)%array)
747  deriv => xc_dset_get_derivative(deriv_set, &
748  [deriv_norm_drhob, deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.true.)
749  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(2)%array)
750  END IF
751  IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
752  cpabort("derivatives bigger than 3 not implemented")
753  END IF
754 
755  DO ispin = 1, 2
756 
757 !$OMP PARALLEL DEFAULT(NONE) &
758 !$OMP SHARED(rho, ispin, rho_1_3, norm_drho, e_0) &
759 !$OMP SHARED(e_rho, e_ndrho, e_rho_rho, e_ndrho_rho) &
760 !$OMP SHARED(e_ndrho_ndrho, e_rho_rho_rho) &
761 !$OMP SHARED(e_ndrho_rho_rho, e_ndrho_ndrho_rho) &
762 !$OMP SHARED(e_ndrho_ndrho_ndrho, grad_deriv, npoints) &
763 !$OMP SHARED(epsilon_rho, sx)
764 
765  CALL xb88_lsd_calc( &
766  rho_spin=rho(ispin)%array, &
767  rho_1_3_spin=rho_1_3(ispin)%array, &
768  norm_drho_spin=norm_drho(ispin)%array, &
769  e_0=e_0, &
770  e_rho_spin=e_rho(ispin)%array, &
771  e_ndrho_spin=e_ndrho(ispin)%array, &
772  e_rho_rho_spin=e_rho_rho(ispin)%array, &
773  e_ndrho_rho_spin=e_ndrho_rho(ispin)%array, &
774  e_ndrho_ndrho_spin=e_ndrho_ndrho(ispin)%array, &
775  e_rho_rho_rho_spin=e_rho_rho_rho(ispin)%array, &
776  e_ndrho_rho_rho_spin=e_ndrho_rho_rho(ispin)%array, &
777  e_ndrho_ndrho_rho_spin=e_ndrho_ndrho_rho(ispin)%array, &
778  e_ndrho_ndrho_ndrho_spin=e_ndrho_ndrho_ndrho(ispin)%array, &
779  grad_deriv=grad_deriv, npoints=npoints, &
780  epsilon_rho=epsilon_rho, sx=sx)
781 
782 !$OMP END PARALLEL
783 
784  END DO
785 
786  CALL timestop(handle)
787 
788  END SUBROUTINE xb88_lsd_eval
789 
790 ! **************************************************************************************************
791 !> \brief low level calculation of the becke 88 exchange functional for lsd
792 !> \param rho_spin alpha or beta spin density
793 !> \param rho_1_3_spin rho_spin**(1./3.)
794 !> \param norm_drho_spin || grad rho_spin ||
795 !> \param e_0 adds to it the local value of the functional
796 !> \param e_rho_spin e_*_spin: derivative of the functional wrt. to the variables
797 !> named where the * is. Everything wrt. to the spin of the arguments.
798 !> \param e_ndrho_spin ...
799 !> \param e_rho_rho_spin ...
800 !> \param e_ndrho_rho_spin ...
801 !> \param e_ndrho_ndrho_spin ...
802 !> \param e_rho_rho_rho_spin ...
803 !> \param e_ndrho_rho_rho_spin ...
804 !> \param e_ndrho_ndrho_rho_spin ...
805 !> \param e_ndrho_ndrho_ndrho_spin ...
806 !> \param grad_deriv ...
807 !> \param npoints ...
808 !> \param epsilon_rho ...
809 !> \param sx scaling-parameter for exchange
810 !> \par History
811 !> 01.2004 created [fawzi]
812 !> 01.2007 added scaling [Manuel Guidon]
813 !> \author fawzi
814 ! **************************************************************************************************
815  SUBROUTINE xb88_lsd_calc(rho_spin, rho_1_3_spin, norm_drho_spin, e_0, &
816  e_rho_spin, e_ndrho_spin, e_rho_rho_spin, e_ndrho_rho_spin, &
817  e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
818  e_ndrho_ndrho_rho_spin, &
819  e_ndrho_ndrho_ndrho_spin, grad_deriv, npoints, epsilon_rho, sx)
820  REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho_spin, rho_1_3_spin, norm_drho_spin
821  REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho_spin, e_ndrho_spin, e_rho_rho_spin, &
822  e_ndrho_rho_spin, e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
823  e_ndrho_ndrho_rho_spin, e_ndrho_ndrho_ndrho_spin
824  INTEGER, INTENT(in) :: grad_deriv, npoints
825  REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
826 
827  INTEGER :: ii
828  REAL(kind=dp) :: c, epsilon_rho43, my_epsilon_rho, my_rho, t1, t103, t11, t12, t127, t133, &
829  t138, t14, t151, t17, t18, t2, t20, t22, t23, t27, t28, t3, t30, t35, t36, t4, t42, t43, &
830  t44, t5, t51, t53, t57, t58, t59, t6, t63, t64, t66, t67, t7, t75, t76, t79, t8, t87, x
831 
832  c = 1.5_dp*(0.75_dp/pi)**(1.0_dp/3.0_dp)
833  my_epsilon_rho = 0.5_dp*epsilon_rho
834  epsilon_rho43 = my_epsilon_rho**(4._dp/3._dp)
835 
836  SELECT CASE (grad_deriv)
837  CASE (0)
838 
839 !$OMP DO
840 
841  DO ii = 1, npoints
842  my_rho = rho_spin(ii)
843  IF (my_rho > my_epsilon_rho) THEN
844  t1 = rho_1_3_spin(ii)*my_rho
845  x = norm_drho_spin(ii)/t1
846  t2 = x**2
847  t3 = beta*t2
848  t4 = beta*x
849  t5 = t2 + 0.1e1_dp
850  t6 = sqrt(t5)
851  t7 = x + t6
852  t8 = log(t7)
853  t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
854  t12 = 0.1e1_dp/t11
855  t14 = -c - t3*t12
856 
857  e_0(ii) = e_0(ii) &
858  + (t1*t14)*sx
859  END IF
860  END DO
861 
862 !$OMP END DO
863 
864  CASE (1)
865 
866 !$OMP DO
867 
868  DO ii = 1, npoints
869  my_rho = rho_spin(ii)
870  IF (my_rho > my_epsilon_rho) THEN
871  t1 = rho_1_3_spin(ii)*my_rho
872  x = norm_drho_spin(ii)/t1
873  t2 = x**2
874  t3 = beta*t2
875  t4 = beta*x
876  t5 = t2 + 0.1e1_dp
877  t6 = sqrt(t5)
878  t7 = x + t6
879  t8 = log(t7)
880  t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
881  t12 = 0.1e1_dp/t11
882  t14 = -c - t3*t12
883 
884  e_0(ii) = e_0(ii) &
885  + (t1*t14)*sx
886 
887  t17 = t11**2
888  t18 = 0.1e1_dp/t17
889  t20 = 0.1e1_dp/t6
890  t22 = 0.1e1_dp + t20*x
891  t23 = 0.1e1_dp/t7
892  t27 = 0.6e1_dp*beta*t8 + 0.6e1_dp*t4*t22*t23
893  t28 = t18*t27
894  t30 = -0.2e1_dp*t4*t12 + t3*t28
895 
896  e_rho_spin(ii) = e_rho_spin(ii) &
897  - (0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t30*x - &
898  0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t14)*sx
899  e_ndrho_spin(ii) = e_ndrho_spin(ii) &
900  + (t30)*sx
901  END IF
902  END DO
903 
904 !$OMP END DO
905 
906  CASE default
907 
908 !$OMP DO
909 
910  DO ii = 1, npoints
911  my_rho = rho_spin(ii)
912  IF (my_rho > my_epsilon_rho) THEN
913  t1 = rho_1_3_spin(ii)*my_rho
914  x = norm_drho_spin(ii)/t1
915  t2 = x**2
916  t3 = beta*t2
917  t4 = beta*x
918  t5 = t2 + 0.1e1_dp
919  t6 = sqrt(t5)
920  t7 = x + t6
921  t8 = log(t7)
922  t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
923  t12 = 0.1e1_dp/t11
924  t14 = -c - t3*t12
925 
926  IF (grad_deriv >= 0) THEN
927  e_0(ii) = e_0(ii) &
928  + (t1*t14)*sx
929  END IF
930 
931  t17 = t11**2
932  t18 = 0.1e1_dp/t17
933  t20 = 0.1e1_dp/t6
934  t22 = 0.1e1_dp + t20*x
935  t23 = 0.1e1_dp/t7
936  t27 = 0.6e1_dp*beta*t8 + 0.6e1_dp*t4*t22*t23
937  t28 = t18*t27
938  t30 = -0.2e1_dp*t4*t12 + t3*t28
939 
940  IF (grad_deriv == -1 .OR. grad_deriv >= 1) THEN
941  e_rho_spin(ii) = e_rho_spin(ii) &
942  - (0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t30*x - &
943  0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t14)*sx
944  e_ndrho_spin(ii) = e_ndrho_spin(ii) &
945  + (t30)*sx
946  END IF
947 
948  t35 = rho_1_3_spin(ii)**2
949  t36 = 0.1e1_dp/t35
950  t42 = 0.1e1_dp/t17/t11
951  t43 = t27**2
952  t44 = t42*t43
953  t51 = 0.1e1_dp/t6/t5
954  t53 = -t51*t2 + t20
955  t57 = t22**2
956  t58 = t7**2
957  t59 = 0.1e1_dp/t58
958  t63 = 0.12e2_dp*beta*t22*t23 + 0.6e1_dp*t4*t53*t23 - 0.6e1_dp*t4*t57*t59
959  t64 = t18*t63
960  t66 = -0.2e1_dp*beta*t12 + 0.4e1_dp*t4*t28 - 0.2e1_dp*t3*t44 + t3*t64
961  t67 = t36*t66
962  t75 = 0.1e1_dp/my_rho
963  t76 = t75*t66
964  t79 = 0.1e1_dp/t1
965 
966  IF (grad_deriv == -2 .OR. grad_deriv >= 2) THEN
967  e_rho_rho_spin(ii) = e_rho_rho_spin(ii) &
968  + (0.16e2_dp/0.9e1_dp*t67*t2 - 0.4e1_dp/0.9e1_dp &
969  *t36*t30*x + 0.4e1_dp/0.9e1_dp*t36*t14)*sx
970  e_ndrho_rho_spin(ii) = e_ndrho_rho_spin(ii) &
971  - (0.4e1_dp/0.3e1_dp*t76*x)*sx
972  e_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_spin(ii) + &
973  (t66*t79)*sx
974  END IF
975 
976  t87 = t17**2
977  t103 = t5**2
978  t127 = 0.6e1_dp*beta*t18*t27 - 0.12e2_dp*t4*t44 + 0.6e1_dp &
979  *t4*t64 + 0.6e1_dp*t3/t87*t43*t27 - 0.6e1_dp*t3 &
980  *t42*t27*t63 + t3*t18*(0.18e2_dp*beta*t53*t23 &
981  - 0.18e2_dp*beta*t57*t59 + 0.6e1_dp*t4*(0.3e1_dp/ &
982  t6/t103*t2*x - 0.3e1_dp*t51*x)*t23 - 0.18e2_dp* &
983  t4*t53*t59*t22 + 0.12e2_dp*t4*t57*t22/t58/t7)
984  t133 = 0.16e2_dp/0.9e1_dp*t36*t127*t2 + 0.28e2_dp/0.9e1_dp &
985  *t67*x
986  t138 = 0.1e1_dp/t35/my_rho
987  t151 = my_rho**2
988 
989  IF (grad_deriv == -3 .OR. grad_deriv >= 3) THEN
990  e_rho_rho_rho_spin(ii) = e_rho_rho_rho_spin(ii) &
991  - (0.4e1_dp/0.3e1_dp*t133*x*t75 + 0.32e2_dp/0.27e2_dp &
992  *t138*t66*t2 - 0.8e1_dp/0.27e2_dp*t138*t30* &
993  x + 0.8e1_dp/0.27e2_dp*t138*t14)*sx
994  e_ndrho_rho_rho_spin(ii) = e_ndrho_rho_rho_spin(ii) &
995  + (t133*t79)*sx
996  e_ndrho_ndrho_rho_spin(ii) = e_ndrho_ndrho_rho_spin(ii) &
997  + ((-0.4e1_dp/0.3e1_dp*t75*t127*x - &
998  0.4e1_dp/0.3e1_dp*t76)*t79)*sx
999  e_ndrho_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_ndrho_spin(ii) &
1000  + (t127/t35/t151)*sx
1001  END IF
1002  END IF
1003  END DO
1004 
1005 !$OMP END DO
1006 
1007  END SELECT
1008 
1009  END SUBROUTINE xb88_lsd_calc
1010 
1011 END MODULE xc_xbecke88
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
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
calculates the Becke 88 exchange functional
Definition: xc_xbecke88.F:14
xc_xbecke88::xb88_lsd_info
subroutine, public xb88_lsd_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition: xc_xbecke88.F:90
xc_xbecke88::xb88_lda_eval
subroutine, public xb88_lda_eval(rho_set, deriv_set, grad_deriv, xb88_params)
evaluates the becke 88 exchange functional for lda
Definition: xc_xbecke88.F:124
xc_xbecke88::xb88_lsd_eval
subroutine, public xb88_lsd_eval(rho_set, deriv_set, grad_deriv, xb88_params)
evaluates the becke 88 exchange functional for lsd
Definition: xc_xbecke88.F:637
xc_xbecke88::xb88_lda_info
subroutine, public xb88_lda_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition: xc_xbecke88.F:59