(git:c641773)
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-2026 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! **************************************************************************************************
15 USE bibliography, ONLY: becke1988,&
16 cite_reference
20 USE kinds, ONLY: dp
21 USE mathconstants, ONLY: pi
25 deriv_rho,&
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
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 maximum degree of the derivative that should be evaluated
115!> \param xb88_params input parameters (scaling)
116!> \par History
117!> 11.2003 created [fawzi]
118!> 01.2007 added scaling [Manuel Guidon]
119!> \author fawzi
120! **************************************************************************************************
121 SUBROUTINE xb88_lda_eval(rho_set, deriv_set, grad_deriv, xb88_params)
122 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
123 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
124 INTEGER, INTENT(in) :: grad_deriv
125 TYPE(section_vals_type), POINTER :: xb88_params
126
127 CHARACTER(len=*), PARAMETER :: routinen = 'xb88_lda_eval'
128
129 INTEGER :: handle, npoints
130 INTEGER, DIMENSION(2, 3) :: bo
131 REAL(kind=dp) :: epsilon_rho, sx
132 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
133 e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
134 e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho, rho_1_3
135 TYPE(xc_derivative_type), POINTER :: deriv
136
137 CALL timeset(routinen, handle)
138
139 CALL section_vals_val_get(xb88_params, "scale_x", r_val=sx)
140
141 CALL cite_reference(becke1988)
142
143 CALL xc_rho_set_get(rho_set, rho_1_3=rho_1_3, rho=rho, &
144 norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
145 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
146
147 dummy => rho
148
149 e_0 => dummy
150 e_rho => dummy
151 e_ndrho => dummy
152 e_rho_rho => dummy
153 e_ndrho_rho => dummy
154 e_ndrho_ndrho => dummy
155 e_rho_rho_rho => dummy
156 e_ndrho_rho_rho => dummy
157 e_ndrho_ndrho_rho => dummy
158 e_ndrho_ndrho_ndrho => dummy
159
160 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
161 allocate_deriv=.true.)
162 CALL xc_derivative_get(deriv, deriv_data=e_0)
163 IF (grad_deriv >= 1) THEN
164 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
165 allocate_deriv=.true.)
166 CALL xc_derivative_get(deriv, deriv_data=e_rho)
167 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
168 allocate_deriv=.true.)
169 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
170 END IF
171 IF (grad_deriv >= 2) THEN
172 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
173 allocate_deriv=.true.)
174 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
175 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
176 allocate_deriv=.true.)
177 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
178 deriv => xc_dset_get_derivative(deriv_set, &
179 [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
180 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
181 END IF
182 IF (grad_deriv >= 3) THEN
183 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
184 allocate_deriv=.true.)
185 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
186 deriv => xc_dset_get_derivative(deriv_set, &
187 [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.true.)
188 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
189 deriv => xc_dset_get_derivative(deriv_set, &
190 [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.true.)
191 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
192 deriv => xc_dset_get_derivative(deriv_set, &
193 [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
194 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
195 END IF
196 IF (grad_deriv > 3) THEN
197 cpabort("derivatives bigger than 3 not implemented")
198 END IF
199
200!$OMP PARALLEL DEFAULT(NONE) &
201!$OMP SHARED(rho, rho_1_3, norm_drho, e_0, e_rho) &
202!$OMP SHARED(e_ndrho, e_rho_rho, e_ndrho_rho) &
203!$OMP SHARED(e_ndrho_ndrho, e_rho_rho_rho) &
204!$OMP SHARED(e_ndrho_rho_rho, e_ndrho_ndrho_rho) &
205!$OMP SHARED(e_ndrho_ndrho_ndrho, grad_deriv, npoints) &
206!$OMP SHARED(epsilon_rho, sx)
207
208 CALL xb88_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
209 e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
210 e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
211 e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
212 e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, &
213 e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, grad_deriv=grad_deriv, &
214 npoints=npoints, epsilon_rho=epsilon_rho, sx=sx)
215
216!$OMP END PARALLEL
217
218 CALL timestop(handle)
219 END SUBROUTINE xb88_lda_eval
220
221! **************************************************************************************************
222!> \brief evaluates the becke 88 exchange functional for lda
223!> \param rho the density where you want to evaluate the functional
224!> \param rho_1_3 ...
225!> \param norm_drho ...
226!> \param e_0 ...
227!> \param e_rho ...
228!> \param e_ndrho ...
229!> \param e_rho_rho ...
230!> \param e_ndrho_rho ...
231!> \param e_ndrho_ndrho ...
232!> \param e_rho_rho_rho ...
233!> \param e_ndrho_rho_rho ...
234!> \param e_ndrho_ndrho_rho ...
235!> \param e_ndrho_ndrho_ndrho ...
236!> \param grad_deriv maximum degree of the derivative that should be evaluated
237!> \param npoints ...
238!> \param epsilon_rho ...
239!> \param sx scaling-parameter for exchange
240!> \par History
241!> 11.2003 created [fawzi]
242!> 01.2007 added scaling [Manuel Guidon]
243!> \author fawzi
244! **************************************************************************************************
245 SUBROUTINE xb88_lda_calc(rho, rho_1_3, norm_drho, &
246 e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
247 e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
248 e_ndrho_ndrho_ndrho, grad_deriv, npoints, epsilon_rho, sx)
249 INTEGER, INTENT(in) :: npoints, grad_deriv
250 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_ndrho, &
251 e_ndrho_ndrho_rho, e_ndrho_rho_rho, e_rho_rho_rho, e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, &
252 e_ndrho, e_rho, e_0
253 REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho_1_3, rho
254 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
255
256 INTEGER :: ii
257 REAL(kind=dp) :: c, epsilon_rho43, my_rho, my_rho_1_3, t1, t10, t100, t104, t11, t12, t126, &
258 t13, t159, t164, t17, t170, t176, t18, t189, t2, t21, t23, t29, t3, t31, t33, t37, t39, &
259 t40, t43, t44, t49, t5, t51, t54, t6, t64, t67, t7, t72, t74, t75, t79, t83, t86, t90, &
260 t91, t98, t99, x
261
262 c = 1.5_dp*(0.75_dp/pi)**(1.0_dp/3.0_dp)
263 t1 = 2**(0.1e1_dp/0.3e1_dp)
264 t2 = t1**2
265 t5 = beta*t2
266 t7 = beta*t1
267 epsilon_rho43 = epsilon_rho**(4._dp/3._dp)
268
269 SELECT CASE (grad_deriv)
270 CASE (0)
271
272!$OMP DO
273
274 DO ii = 1, npoints
275
276 my_rho = rho(ii)
277 IF (my_rho > epsilon_rho) THEN
278 my_rho_1_3 = rho_1_3(ii)
279 t3 = my_rho_1_3*my_rho
280 x = norm_drho(ii)/max(t3, epsilon_rho43)
281 t6 = x**2
282 t10 = t2*t6 + 0.1e1_dp
283 t11 = sqrt(t10)
284 t12 = t1*x + t11
285 t13 = log(t12)
286 t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
287 t18 = 0.1e1_dp/t17
288 t21 = -c - t5*t6*t18
289
290 e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
291 END IF
292 END DO
293
294!$OMP END DO
295
296 CASE (1)
297
298!$OMP DO
299
300 DO ii = 1, npoints
301
302 my_rho = rho(ii)
303 IF (my_rho > epsilon_rho) THEN
304 my_rho_1_3 = rho_1_3(ii)
305 t3 = my_rho_1_3*my_rho
306 x = norm_drho(ii)/t3
307 t6 = x**2
308 t10 = t2*t6 + 0.1e1_dp
309 t11 = sqrt(t10)
310 t12 = t1*x + t11
311 t13 = log(t12)
312 t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
313 t18 = 0.1e1_dp/t17
314 t21 = -c - t5*t6*t18
315
316 e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
317
318 t23 = t2*my_rho_1_3
319 t29 = 0.1e1_dp/t11
320 t31 = t1 + t2*x*t29
321 t33 = 0.1e1_dp/t12
322 t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
323 t39 = t17**2
324 t40 = 0.1e1_dp/t39
325 t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
326 t44 = t43*x
327
328 e_rho(ii) = e_rho(ii) &
329 - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
330 t23*t21)*sx
331 e_ndrho(ii) = e_ndrho(ii) &
332 + (t2*t43/0.2e1_dp)*sx
333 END IF
334
335 END DO
336
337!$OMP END DO
338
339 CASE (2)
340
341!$OMP DO
342
343 DO ii = 1, npoints
344
345 my_rho = rho(ii)
346 IF (my_rho > epsilon_rho) THEN
347 my_rho_1_3 = rho_1_3(ii)
348 t3 = my_rho_1_3*my_rho
349 x = norm_drho(ii)/t3
350 t6 = x**2
351 t10 = t2*t6 + 0.1e1_dp
352 t11 = sqrt(t10)
353 t12 = t1*x + t11
354 t13 = log(t12)
355 t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
356 t18 = 0.1e1_dp/t17
357 t21 = -c - t5*t6*t18
358
359 e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
360
361 t23 = t2*my_rho_1_3
362 t29 = 0.1e1_dp/t11
363 t31 = t1 + t2*x*t29
364 t33 = 0.1e1_dp/t12
365 t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
366 t39 = t17**2
367 t40 = 0.1e1_dp/t39
368 t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
369 t44 = t43*x
370
371 e_rho(ii) = e_rho(ii) &
372 - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
373 t23*t21)*sx
374 e_ndrho(ii) = e_ndrho(ii) &
375 + (t2*t43/0.2e1_dp)*sx
376
377 t49 = my_rho_1_3**2
378 t51 = t2/t49
379 t54 = x*t40
380 t64 = 0.1e1_dp/t11/t10
381 t67 = t2*t29 - 0.2e1_dp*t1*t6*t64
382 t72 = t31**2
383 t74 = t12**2
384 t75 = 0.1e1_dp/t74
385 t79 = 0.12e2_dp*t7*t31*t33 + 0.6e1_dp*t7*x*t67*t33 &
386 - 0.6e1_dp*t7*x*t72*t75
387 t83 = t37**2
388 t86 = 0.1e1_dp/t39/t17
389 t90 = -0.2e1_dp*t5*t18 + 0.4e1_dp*t5*t54*t37 + t5*t6* &
390 t79*t40 - 0.2e1_dp*t5*t6*t83*t86
391 t91 = t90*t6
392 t98 = 0.1e1_dp/my_rho
393 t99 = t2*t98
394 t100 = t90*x
395 t104 = 0.1e1_dp/t3
396
397 e_rho_rho(ii) = e_rho_rho(ii) &
398 + (0.8e1_dp/0.9e1_dp*t51*t91 - 0.2e1_dp/0.9e1_dp* &
399 t51*t44 + 0.2e1_dp/0.9e1_dp*t51*t21)*sx
400 e_ndrho_rho(ii) = e_ndrho_rho(ii) &
401 - (0.2e1_dp/0.3e1_dp*t99*t100)*sx
402 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
403 + (t2*t90*t104/0.2e1_dp)*sx
404 END IF
405 END DO
406
407!$OMP END DO
408
409 CASE default
410
411!$OMP DO
412
413 DO ii = 1, npoints
414
415 my_rho = rho(ii)
416 IF (my_rho > epsilon_rho) THEN
417 my_rho_1_3 = rho_1_3(ii)
418 t3 = my_rho_1_3*my_rho
419 x = norm_drho(ii)/t3
420 t6 = x**2
421 t10 = t2*t6 + 0.1e1_dp
422 t11 = sqrt(t10)
423 t12 = t1*x + t11
424 t13 = log(t12)
425 t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
426 t18 = 0.1e1_dp/t17
427 t21 = -c - t5*t6*t18
428
429 e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
430
431 IF (grad_deriv >= 1) THEN
432 t23 = t2*my_rho_1_3
433 t29 = 0.1e1_dp/t11
434 t31 = t1 + t2*x*t29
435 t33 = 0.1e1_dp/t12
436 t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
437 t39 = t17**2
438 t40 = 0.1e1_dp/t39
439 t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
440 t44 = t43*x
441
442 e_rho(ii) = e_rho(ii) &
443 - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
444 t23*t21)*sx
445 e_ndrho(ii) = e_ndrho(ii) &
446 + (t2*t43/0.2e1_dp)*sx
447 END IF
448
449 IF (grad_deriv >= 2) THEN
450 t49 = my_rho_1_3**2
451 t51 = t2/t49
452 t54 = x*t40
453 t64 = 0.1e1_dp/t11/t10
454 t67 = t2*t29 - 0.2e1_dp*t1*t6*t64
455 t72 = t31**2
456 t74 = t12**2
457 t75 = 0.1e1_dp/t74
458 t79 = 0.12e2_dp*t7*t31*t33 + 0.6e1_dp*t7*x*t67*t33 &
459 - 0.6e1_dp*t7*x*t72*t75
460 t83 = t37**2
461 t86 = 0.1e1_dp/t39/t17
462 t90 = -0.2e1_dp*t5*t18 + 0.4e1_dp*t5*t54*t37 + t5*t6* &
463 t79*t40 - 0.2e1_dp*t5*t6*t83*t86
464 t91 = t90*t6
465 t98 = 0.1e1_dp/my_rho
466 t99 = t2*t98
467 t100 = t90*x
468 t104 = 0.1e1_dp/t3
469
470 e_rho_rho(ii) = e_rho_rho(ii) &
471 + (0.8e1_dp/0.9e1_dp*t51*t91 - 0.2e1_dp/0.9e1_dp* &
472 t51*t44 + 0.2e1_dp/0.9e1_dp*t51*t21)*sx
473 e_ndrho_rho(ii) = e_ndrho_rho(ii) &
474 - (0.2e1_dp/0.3e1_dp*t99*t100)*sx
475 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
476 + (t2*t90*t104/0.2e1_dp)*sx
477 END IF
478
479 IF (grad_deriv >= 3) THEN
480 t126 = t10**2
481 t159 = t39**2
482 t164 = 0.6e1_dp*t5*t40*t37 - 0.12e2_dp*t5*x*t86*t83 &
483 + 0.6e1_dp*t5*t54*t79 + t5*t6*(0.18e2_dp*t7*t67 &
484 *t33 - 0.18e2_dp*t7*t72*t75 + 0.6e1_dp*t7*x* &
485 (-0.6e1_dp*t1*t64*x + 0.12e2_dp*t6*x/t11/t126) &
486 *t33 - 0.18e2_dp*t7*x*t67*t75*t31 + 0.12e2_dp*t7 &
487 *x*t72*t31/t74/t12)*t40 - 0.6e1_dp*t5*t6*t79 &
488 *t86*t37 + 0.6e1_dp*t5*t6*t83*t37/t159
489 t170 = 0.8e1_dp/0.9e1_dp*t51*t164*t6 + 0.14e2_dp/0.9e1_dp &
490 *t51*t100
491 t176 = t2/t49/my_rho
492 t189 = my_rho**2
493
494 e_rho_rho_rho(ii) = e_rho_rho_rho(ii) &
495 - (0.4e1_dp/0.3e1_dp*t170*x*t98 + 0.16e2_dp/0.27e2_dp &
496 *t176*t91 - 0.4e1_dp/0.27e2_dp*t176*t44 + &
497 0.4e1_dp/0.27e2_dp*t176*t21)*sx
498 e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) &
499 + (t170*t104)*sx
500 e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) &
501 + ((-0.2e1_dp/0.3e1_dp*t99*t164*x - &
502 0.2e1_dp/0.3e1_dp*t99*t90)*t104)*sx
503 e_ndrho_ndrho_ndrho(ii) = e_ndrho_ndrho_ndrho(ii) &
504 + (t2*t164/t49/t189/0.2e1_dp)*sx
505 END IF
506 END IF
507 END DO
508
509!$OMP END DO
510
511 END SELECT
512
513 END SUBROUTINE xb88_lda_calc
514
515! **************************************************************************************************
516!> \brief evaluates the becke 88 exchange functional for lsd
517!> \param rho_set ...
518!> \param deriv_set ...
519!> \param grad_deriv ...
520!> \param xb88_params ...
521!> \par History
522!> 11.2003 created [fawzi]
523!> 01.2007 added scaling [Manuel Guidon]
524!> \author fawzi
525! **************************************************************************************************
526 SUBROUTINE xb88_lsd_eval(rho_set, deriv_set, grad_deriv, xb88_params)
527 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
528 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
529 INTEGER, INTENT(in) :: grad_deriv
530 TYPE(section_vals_type), POINTER :: xb88_params
531
532 CHARACTER(len=*), PARAMETER :: routinen = 'xb88_lsd_eval'
533
534 INTEGER :: handle, i, ispin, npoints
535 INTEGER, DIMENSION(2, 3) :: bo
536 REAL(kind=dp) :: epsilon_rho, sx
537 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
538 POINTER :: dummy, e_0
539 TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
540 e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, e_rho_rho_rho, &
541 norm_drho, rho, rho_1_3
542 TYPE(xc_derivative_type), POINTER :: deriv
543
544 CALL timeset(routinen, handle)
545
546 CALL cite_reference(becke1988)
547
548 NULLIFY (deriv)
549 DO i = 1, 2
550 NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
551 END DO
552
553 CALL section_vals_val_get(xb88_params, "scale_x", r_val=sx)
554 CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
555 rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
556 rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
557 norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
558 local_bounds=bo)
559 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
560
561 dummy => rho(1)%array
562
563 e_0 => dummy
564 DO i = 1, 2
565 e_rho(i)%array => dummy
566 e_ndrho(i)%array => dummy
567 e_rho_rho(i)%array => dummy
568 e_ndrho_rho(i)%array => dummy
569 e_ndrho_ndrho(i)%array => dummy
570 e_rho_rho_rho(i)%array => dummy
571 e_ndrho_rho_rho(i)%array => dummy
572 e_ndrho_ndrho_rho(i)%array => dummy
573 e_ndrho_ndrho_ndrho(i)%array => dummy
574 END DO
575
576 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
577 allocate_deriv=.true.)
578 CALL xc_derivative_get(deriv, deriv_data=e_0)
579 IF (grad_deriv >= 1) THEN
580 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
581 allocate_deriv=.true.)
582 CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
583 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
584 allocate_deriv=.true.)
585 CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
586 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
587 allocate_deriv=.true.)
588 CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
589 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
590 allocate_deriv=.true.)
591 CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
592 END IF
593 IF (grad_deriv >= 2) THEN
594 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
595 allocate_deriv=.true.)
596 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(1)%array)
597 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
598 allocate_deriv=.true.)
599 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(2)%array)
600 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
601 allocate_deriv=.true.)
602 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(1)%array)
603 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
604 allocate_deriv=.true.)
605 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(2)%array)
606 deriv => xc_dset_get_derivative(deriv_set, &
607 [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.true.)
608 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(1)%array)
609 deriv => xc_dset_get_derivative(deriv_set, &
610 [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.true.)
611 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(2)%array)
612 END IF
613 IF (grad_deriv >= 3) THEN
614 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
615 allocate_deriv=.true.)
616 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(1)%array)
617 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
618 allocate_deriv=.true.)
619 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(2)%array)
620 deriv => xc_dset_get_derivative(deriv_set, &
621 [deriv_norm_drhoa, deriv_rhoa, deriv_rhoa], allocate_deriv=.true.)
622 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(1)%array)
623 deriv => xc_dset_get_derivative(deriv_set, &
624 [deriv_norm_drhob, deriv_rhob, deriv_rhob], allocate_deriv=.true.)
625 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(2)%array)
626 deriv => xc_dset_get_derivative(deriv_set, &
627 [deriv_norm_drhoa, deriv_norm_drhoa, deriv_rhoa], allocate_deriv=.true.)
628 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(1)%array)
629 deriv => xc_dset_get_derivative(deriv_set, &
630 [deriv_norm_drhob, deriv_norm_drhob, deriv_rhob], allocate_deriv=.true.)
631 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(2)%array)
632 deriv => xc_dset_get_derivative(deriv_set, &
633 [deriv_norm_drhoa, deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.true.)
634 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(1)%array)
635 deriv => xc_dset_get_derivative(deriv_set, &
636 [deriv_norm_drhob, deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.true.)
637 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(2)%array)
638 END IF
639 IF (grad_deriv > 3) THEN
640 cpabort("derivatives bigger than 3 not implemented")
641 END IF
642
643 DO ispin = 1, 2
644
645!$OMP PARALLEL DEFAULT(NONE) &
646!$OMP SHARED(rho, ispin, rho_1_3, norm_drho, e_0) &
647!$OMP SHARED(e_rho, e_ndrho, e_rho_rho, e_ndrho_rho) &
648!$OMP SHARED(e_ndrho_ndrho, e_rho_rho_rho) &
649!$OMP SHARED(e_ndrho_rho_rho, e_ndrho_ndrho_rho) &
650!$OMP SHARED(e_ndrho_ndrho_ndrho, grad_deriv, npoints) &
651!$OMP SHARED(epsilon_rho, sx)
652
653 CALL xb88_lsd_calc( &
654 rho_spin=rho(ispin)%array, &
655 rho_1_3_spin=rho_1_3(ispin)%array, &
656 norm_drho_spin=norm_drho(ispin)%array, &
657 e_0=e_0, &
658 e_rho_spin=e_rho(ispin)%array, &
659 e_ndrho_spin=e_ndrho(ispin)%array, &
660 e_rho_rho_spin=e_rho_rho(ispin)%array, &
661 e_ndrho_rho_spin=e_ndrho_rho(ispin)%array, &
662 e_ndrho_ndrho_spin=e_ndrho_ndrho(ispin)%array, &
663 e_rho_rho_rho_spin=e_rho_rho_rho(ispin)%array, &
664 e_ndrho_rho_rho_spin=e_ndrho_rho_rho(ispin)%array, &
665 e_ndrho_ndrho_rho_spin=e_ndrho_ndrho_rho(ispin)%array, &
666 e_ndrho_ndrho_ndrho_spin=e_ndrho_ndrho_ndrho(ispin)%array, &
667 grad_deriv=grad_deriv, npoints=npoints, &
668 epsilon_rho=epsilon_rho, sx=sx)
669
670!$OMP END PARALLEL
671
672 END DO
673
674 CALL timestop(handle)
675
676 END SUBROUTINE xb88_lsd_eval
677
678! **************************************************************************************************
679!> \brief low level calculation of the becke 88 exchange functional for lsd
680!> \param rho_spin alpha or beta spin density
681!> \param rho_1_3_spin rho_spin**(1./3.)
682!> \param norm_drho_spin || grad rho_spin ||
683!> \param e_0 adds to it the local value of the functional
684!> \param e_rho_spin e_*_spin: derivative of the functional wrt. to the variables
685!> named where the * is. Everything wrt. to the spin of the arguments.
686!> \param e_ndrho_spin ...
687!> \param e_rho_rho_spin ...
688!> \param e_ndrho_rho_spin ...
689!> \param e_ndrho_ndrho_spin ...
690!> \param e_rho_rho_rho_spin ...
691!> \param e_ndrho_rho_rho_spin ...
692!> \param e_ndrho_ndrho_rho_spin ...
693!> \param e_ndrho_ndrho_ndrho_spin ...
694!> \param grad_deriv ...
695!> \param npoints ...
696!> \param epsilon_rho ...
697!> \param sx scaling-parameter for exchange
698!> \par History
699!> 01.2004 created [fawzi]
700!> 01.2007 added scaling [Manuel Guidon]
701!> \author fawzi
702! **************************************************************************************************
703 SUBROUTINE xb88_lsd_calc(rho_spin, rho_1_3_spin, norm_drho_spin, e_0, &
704 e_rho_spin, e_ndrho_spin, e_rho_rho_spin, e_ndrho_rho_spin, &
705 e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
706 e_ndrho_ndrho_rho_spin, &
707 e_ndrho_ndrho_ndrho_spin, grad_deriv, npoints, epsilon_rho, sx)
708 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho_spin, rho_1_3_spin, norm_drho_spin
709 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho_spin, e_ndrho_spin, e_rho_rho_spin, &
710 e_ndrho_rho_spin, e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
711 e_ndrho_ndrho_rho_spin, e_ndrho_ndrho_ndrho_spin
712 INTEGER, INTENT(in) :: grad_deriv, npoints
713 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
714
715 INTEGER :: ii
716 REAL(kind=dp) :: c, epsilon_rho43, my_epsilon_rho, my_rho, t1, t103, t11, t12, t127, t133, &
717 t138, t14, t151, t17, t18, t2, t20, t22, t23, t27, t28, t3, t30, t35, t36, t4, t42, t43, &
718 t44, t5, t51, t53, t57, t58, t59, t6, t63, t64, t66, t67, t7, t75, t76, t79, t8, t87, x
719
720 c = 1.5_dp*(0.75_dp/pi)**(1.0_dp/3.0_dp)
721 my_epsilon_rho = 0.5_dp*epsilon_rho
722 epsilon_rho43 = my_epsilon_rho**(4._dp/3._dp)
723
724 SELECT CASE (grad_deriv)
725 CASE (0)
726
727!$OMP DO
728
729 DO ii = 1, npoints
730 my_rho = rho_spin(ii)
731 IF (my_rho > my_epsilon_rho) THEN
732 t1 = rho_1_3_spin(ii)*my_rho
733 x = norm_drho_spin(ii)/t1
734 t2 = x**2
735 t3 = beta*t2
736 t4 = beta*x
737 t5 = t2 + 0.1e1_dp
738 t6 = sqrt(t5)
739 t7 = x + t6
740 t8 = log(t7)
741 t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
742 t12 = 0.1e1_dp/t11
743 t14 = -c - t3*t12
744
745 e_0(ii) = e_0(ii) &
746 + (t1*t14)*sx
747 END IF
748 END DO
749
750!$OMP END DO
751
752 CASE (1)
753
754!$OMP DO
755
756 DO ii = 1, npoints
757 my_rho = rho_spin(ii)
758 IF (my_rho > my_epsilon_rho) THEN
759 t1 = rho_1_3_spin(ii)*my_rho
760 x = norm_drho_spin(ii)/t1
761 t2 = x**2
762 t3 = beta*t2
763 t4 = beta*x
764 t5 = t2 + 0.1e1_dp
765 t6 = sqrt(t5)
766 t7 = x + t6
767 t8 = log(t7)
768 t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
769 t12 = 0.1e1_dp/t11
770 t14 = -c - t3*t12
771
772 e_0(ii) = e_0(ii) &
773 + (t1*t14)*sx
774
775 t17 = t11**2
776 t18 = 0.1e1_dp/t17
777 t20 = 0.1e1_dp/t6
778 t22 = 0.1e1_dp + t20*x
779 t23 = 0.1e1_dp/t7
780 t27 = 0.6e1_dp*beta*t8 + 0.6e1_dp*t4*t22*t23
781 t28 = t18*t27
782 t30 = -0.2e1_dp*t4*t12 + t3*t28
783
784 e_rho_spin(ii) = e_rho_spin(ii) &
785 - (0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t30*x - &
786 0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t14)*sx
787 e_ndrho_spin(ii) = e_ndrho_spin(ii) &
788 + (t30)*sx
789 END IF
790 END DO
791
792!$OMP END DO
793
794 CASE default
795
796!$OMP DO
797
798 DO ii = 1, npoints
799 my_rho = rho_spin(ii)
800 IF (my_rho > my_epsilon_rho) THEN
801 t1 = rho_1_3_spin(ii)*my_rho
802 x = norm_drho_spin(ii)/t1
803 t2 = x**2
804 t3 = beta*t2
805 t4 = beta*x
806 t5 = t2 + 0.1e1_dp
807 t6 = sqrt(t5)
808 t7 = x + t6
809 t8 = log(t7)
810 t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
811 t12 = 0.1e1_dp/t11
812 t14 = -c - t3*t12
813
814 e_0(ii) = e_0(ii) &
815 + (t1*t14)*sx
816
817 IF (grad_deriv >= 1) THEN
818 t17 = t11**2
819 t18 = 0.1e1_dp/t17
820 t20 = 0.1e1_dp/t6
821 t22 = 0.1e1_dp + t20*x
822 t23 = 0.1e1_dp/t7
823 t27 = 0.6e1_dp*beta*t8 + 0.6e1_dp*t4*t22*t23
824 t28 = t18*t27
825 t30 = -0.2e1_dp*t4*t12 + t3*t28
826
827 e_rho_spin(ii) = e_rho_spin(ii) &
828 - (0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t30*x - &
829 0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t14)*sx
830 e_ndrho_spin(ii) = e_ndrho_spin(ii) &
831 + (t30)*sx
832 END IF
833
834 IF (grad_deriv >= 2) THEN
835 t35 = rho_1_3_spin(ii)**2
836 t36 = 0.1e1_dp/t35
837 t42 = 0.1e1_dp/t17/t11
838 t43 = t27**2
839 t44 = t42*t43
840 t51 = 0.1e1_dp/t6/t5
841 t53 = -t51*t2 + t20
842 t57 = t22**2
843 t58 = t7**2
844 t59 = 0.1e1_dp/t58
845 t63 = 0.12e2_dp*beta*t22*t23 + 0.6e1_dp*t4*t53*t23 - 0.6e1_dp*t4*t57*t59
846 t64 = t18*t63
847 t66 = -0.2e1_dp*beta*t12 + 0.4e1_dp*t4*t28 - 0.2e1_dp*t3*t44 + t3*t64
848 t67 = t36*t66
849 t75 = 0.1e1_dp/my_rho
850 t76 = t75*t66
851 t79 = 0.1e1_dp/t1
852
853 e_rho_rho_spin(ii) = e_rho_rho_spin(ii) &
854 + (0.16e2_dp/0.9e1_dp*t67*t2 - 0.4e1_dp/0.9e1_dp &
855 *t36*t30*x + 0.4e1_dp/0.9e1_dp*t36*t14)*sx
856 e_ndrho_rho_spin(ii) = e_ndrho_rho_spin(ii) &
857 - (0.4e1_dp/0.3e1_dp*t76*x)*sx
858 e_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_spin(ii) + &
859 (t66*t79)*sx
860 END IF
861
862 IF (grad_deriv >= 3) THEN
863 t87 = t17**2
864 t103 = t5**2
865 t127 = 0.6e1_dp*beta*t18*t27 - 0.12e2_dp*t4*t44 + 0.6e1_dp &
866 *t4*t64 + 0.6e1_dp*t3/t87*t43*t27 - 0.6e1_dp*t3 &
867 *t42*t27*t63 + t3*t18*(0.18e2_dp*beta*t53*t23 &
868 - 0.18e2_dp*beta*t57*t59 + 0.6e1_dp*t4*(0.3e1_dp/ &
869 t6/t103*t2*x - 0.3e1_dp*t51*x)*t23 - 0.18e2_dp* &
870 t4*t53*t59*t22 + 0.12e2_dp*t4*t57*t22/t58/t7)
871 t133 = 0.16e2_dp/0.9e1_dp*t36*t127*t2 + 0.28e2_dp/0.9e1_dp &
872 *t67*x
873 t138 = 0.1e1_dp/t35/my_rho
874 t151 = my_rho**2
875
876 e_rho_rho_rho_spin(ii) = e_rho_rho_rho_spin(ii) &
877 - (0.4e1_dp/0.3e1_dp*t133*x*t75 + 0.32e2_dp/0.27e2_dp &
878 *t138*t66*t2 - 0.8e1_dp/0.27e2_dp*t138*t30* &
879 x + 0.8e1_dp/0.27e2_dp*t138*t14)*sx
880 e_ndrho_rho_rho_spin(ii) = e_ndrho_rho_rho_spin(ii) &
881 + (t133*t79)*sx
882 e_ndrho_ndrho_rho_spin(ii) = e_ndrho_ndrho_rho_spin(ii) &
883 + ((-0.4e1_dp/0.3e1_dp*t75*t127*x - &
884 0.4e1_dp/0.3e1_dp*t76)*t79)*sx
885 e_ndrho_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_ndrho_spin(ii) &
886 + (t127/t35/t151)*sx
887 END IF
888 END IF
889 END DO
890
891!$OMP END DO
892
893 END SELECT
894
895 END SUBROUTINE xb88_lsd_calc
896
897END MODULE xc_xbecke88
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public becke1988
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
objects that represent the structure of input sections and the data contained in an input section
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
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Module with functions to handle derivative descriptors. derivative description are strings have the f...
integer, parameter, public deriv_norm_drho
integer, parameter, public deriv_norm_drhoa
integer, parameter, public deriv_rhob
integer, parameter, public deriv_rhoa
integer, parameter, public deriv_rho
integer, parameter, public deriv_norm_drhob
represent a group ofunctional derivatives
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Provides types for the management of the xc-functionals and their derivatives.
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
contains the structure
contains the structure
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
calculates the Becke 88 exchange functional
Definition xc_xbecke88.F:14
subroutine, public xb88_lsd_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_xbecke88.F:90
subroutine, public xb88_lda_eval(rho_set, deriv_set, grad_deriv, xb88_params)
evaluates the becke 88 exchange functional for lda
subroutine, public xb88_lsd_eval(rho_set, deriv_set, grad_deriv, xb88_params)
evaluates the becke 88 exchange functional for lsd
subroutine, public xb88_lda_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_xbecke88.F:59
represent a pointer to a contiguous 3d array
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
represent a derivative of a functional
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
represent a density, with all the representation and data needed to perform a functional evaluation