(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
14MODULE 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
45CONTAINS
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
1011END 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
cp_array_utils::cp_3d_r_cp_type
represent a pointer to a contiguous 3d array
Definition cp_array_utils.F:149
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:49
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78