Line data Source code
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 Calculate the local exchange functional
10 : !> \note
11 : !> Order of derivatives is: LDA 0; 1; 2; 3;
12 : !> LSD 0; a b; aa bb; aaa bbb;
13 : !> \par History
14 : !> JGH (26.02.2003) : OpenMP enabled
15 : !> fawzi (04.2004) : adapted to the new xc interface
16 : !> MG (01.2007) : added scaling
17 : !> \author JGH (17.02.2002)
18 : ! **************************************************************************************************
19 : MODULE xc_xalpha
20 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
21 : USE input_section_types, ONLY: section_vals_type,&
22 : section_vals_val_get
23 : USE kinds, ONLY: dp
24 : USE mathconstants, ONLY: pi
25 : USE pw_types, ONLY: pw_r3d_rs_type
26 : USE xc_derivative_desc, ONLY: deriv_rho,&
27 : deriv_rhoa,&
28 : deriv_rhob
29 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
30 : xc_dset_get_derivative
31 : USE xc_derivative_types, ONLY: xc_derivative_get,&
32 : xc_derivative_type
33 : USE xc_functionals_utilities, ONLY: set_util
34 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
35 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
36 : xc_rho_set_type
37 : #include "../base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 :
41 : PRIVATE
42 :
43 : ! *** Global parameters ***
44 : REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
45 : f23 = 2.0_dp*f13, &
46 : f43 = 4.0_dp*f13
47 :
48 : PUBLIC :: xalpha_info, xalpha_lda_eval, xalpha_lsd_eval, xalpha_fxc_eval
49 :
50 : REAL(KIND=dp) :: xparam, flda, flsd
51 : REAL(KIND=dp) :: eps_rho
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xalpha'
53 :
54 : CONTAINS
55 :
56 : ! **************************************************************************************************
57 : !> \brief ...
58 : !> \param cutoff ...
59 : !> \param xalpha ...
60 : ! **************************************************************************************************
61 2799 : SUBROUTINE xalpha_init(cutoff, xalpha)
62 :
63 : REAL(KIND=dp), INTENT(IN) :: cutoff
64 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: xalpha
65 :
66 2799 : eps_rho = cutoff
67 2799 : CALL set_util(cutoff)
68 2799 : IF (PRESENT(xalpha)) THEN
69 2799 : xparam = xalpha
70 : ELSE
71 0 : xparam = 2.0_dp/3.0_dp
72 : END IF
73 :
74 2799 : flda = -9.0_dp/8.0_dp*xparam*(3.0_dp/pi)**f13
75 2799 : flsd = flda*2.0_dp**f13
76 :
77 2799 : END SUBROUTINE xalpha_init
78 :
79 : ! **************************************************************************************************
80 : !> \brief ...
81 : !> \param lsd ...
82 : !> \param reference ...
83 : !> \param shortform ...
84 : !> \param needs ...
85 : !> \param max_deriv ...
86 : !> \param xa_parameter ...
87 : !> \param scaling ...
88 : ! **************************************************************************************************
89 2760 : SUBROUTINE xalpha_info(lsd, reference, shortform, needs, max_deriv, &
90 : xa_parameter, scaling)
91 : LOGICAL, INTENT(in) :: lsd
92 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
93 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
94 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
95 : REAL(KIND=dp), INTENT(in), OPTIONAL :: xa_parameter, scaling
96 :
97 : REAL(KIND=dp) :: my_scaling, my_xparam
98 :
99 2760 : my_xparam = 2.0_dp/3.0_dp
100 2760 : IF (PRESENT(xa_parameter)) my_xparam = xa_parameter
101 2760 : my_scaling = 1.0_dp
102 2760 : IF (PRESENT(scaling)) my_scaling = scaling
103 :
104 2760 : IF (PRESENT(reference)) THEN
105 39 : IF (my_scaling /= 1._dp) THEN
106 : WRITE (reference, '(A,F8.4,A,F8.4)') &
107 0 : "Dirac/Slater local exchange; parameter=", my_xparam, " scaling=", my_scaling
108 : ELSE
109 : WRITE (reference, '(A,F8.4)') &
110 39 : "Dirac/Slater local exchange; parameter=", my_xparam
111 : END IF
112 39 : IF (.NOT. lsd) THEN
113 34 : IF (LEN_TRIM(reference) + 6 < LEN(reference)) THEN
114 34 : reference(LEN_TRIM(reference):LEN_TRIM(reference) + 6) = ' {LDA}'
115 : END IF
116 : END IF
117 : END IF
118 2760 : IF (PRESENT(shortform)) THEN
119 39 : IF (my_scaling /= 1._dp) THEN
120 0 : WRITE (shortform, '(A,F8.4,F8.4)') "Dirac/Slater exchange", my_xparam, my_scaling
121 : ELSE
122 39 : WRITE (shortform, '(A,F8.4)') "Dirac/Slater exchange", my_xparam
123 : END IF
124 39 : IF (.NOT. lsd) THEN
125 34 : IF (LEN_TRIM(shortform) + 6 < LEN(shortform)) THEN
126 34 : shortform(LEN_TRIM(shortform):LEN_TRIM(shortform) + 6) = ' {LDA}'
127 : END IF
128 : END IF
129 : END IF
130 2760 : IF (PRESENT(needs)) THEN
131 2721 : IF (lsd) THEN
132 130 : needs%rho_spin = .TRUE.
133 130 : needs%rho_spin_1_3 = .TRUE.
134 : ELSE
135 2591 : needs%rho = .TRUE.
136 2591 : needs%rho_1_3 = .TRUE.
137 : END IF
138 : END IF
139 2760 : IF (PRESENT(max_deriv)) max_deriv = 3
140 :
141 2760 : END SUBROUTINE xalpha_info
142 :
143 : ! **************************************************************************************************
144 : !> \brief ...
145 : !> \param rho_set ...
146 : !> \param deriv_set ...
147 : !> \param order ...
148 : !> \param xa_params ...
149 : !> \param xa_parameter ...
150 : ! **************************************************************************************************
151 7899 : SUBROUTINE xalpha_lda_eval(rho_set, deriv_set, order, xa_params, xa_parameter)
152 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
153 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
154 : INTEGER, INTENT(in) :: order
155 : TYPE(section_vals_type), POINTER :: xa_params
156 : REAL(KIND=dp), INTENT(in), OPTIONAL :: xa_parameter
157 :
158 : CHARACTER(len=*), PARAMETER :: routineN = 'xalpha_lda_eval'
159 :
160 : INTEGER :: handle, npoints
161 : INTEGER, DIMENSION(2, 3) :: bo
162 : REAL(KIND=dp) :: epsilon_rho, sx
163 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
164 2633 : POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho, &
165 2633 : r13, rho
166 : TYPE(xc_derivative_type), POINTER :: deriv
167 :
168 2633 : CALL timeset(routineN, handle)
169 :
170 2633 : CALL section_vals_val_get(xa_params, "scale_x", r_val=sx)
171 :
172 : CALL xc_rho_set_get(rho_set, rho_1_3=r13, rho=rho, &
173 2633 : local_bounds=bo, rho_cutoff=epsilon_rho)
174 2633 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
175 2633 : CALL xalpha_init(epsilon_rho, xa_parameter)
176 :
177 2633 : IF (order >= 0) THEN
178 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
179 2633 : allocate_deriv=.TRUE.)
180 2633 : CALL xc_derivative_get(deriv, deriv_data=e_0)
181 :
182 2633 : CALL xalpha_lda_0(npoints, rho, r13, e_0, sx)
183 :
184 : END IF
185 2633 : IF (order >= 1 .OR. order == -1) THEN
186 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
187 2623 : allocate_deriv=.TRUE.)
188 2623 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
189 :
190 2623 : CALL xalpha_lda_1(npoints, rho, r13, e_rho, sx)
191 : END IF
192 2633 : IF (order >= 2 .OR. order == -2) THEN
193 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
194 380 : allocate_deriv=.TRUE.)
195 380 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
196 :
197 380 : CALL xalpha_lda_2(npoints, rho, r13, e_rho_rho, sx)
198 : END IF
199 2633 : IF (order >= 3 .OR. order == -3) THEN
200 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
201 0 : allocate_deriv=.TRUE.)
202 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
203 :
204 0 : CALL xalpha_lda_3(npoints, rho, r13, e_rho_rho_rho, sx)
205 : END IF
206 2633 : IF (order > 3 .OR. order < -3) THEN
207 0 : CPABORT("derivatives bigger than 3 not implemented")
208 : END IF
209 2633 : CALL timestop(handle)
210 :
211 2633 : END SUBROUTINE xalpha_lda_eval
212 :
213 : ! **************************************************************************************************
214 : !> \brief ...
215 : !> \param rho_set ...
216 : !> \param deriv_set ...
217 : !> \param order ...
218 : !> \param xa_params ...
219 : !> \param xa_parameter ...
220 : ! **************************************************************************************************
221 664 : SUBROUTINE xalpha_lsd_eval(rho_set, deriv_set, order, xa_params, xa_parameter)
222 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
223 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
224 : INTEGER, INTENT(in) :: order
225 : TYPE(section_vals_type), POINTER :: xa_params
226 : REAL(KIND=dp), INTENT(in), OPTIONAL :: xa_parameter
227 :
228 : CHARACTER(len=*), PARAMETER :: routineN = 'xalpha_lsd_eval'
229 : INTEGER, DIMENSION(2), PARAMETER :: rho_spin_name = [deriv_rhoa, deriv_rhob]
230 :
231 : INTEGER :: handle, i, ispin, npoints
232 : INTEGER, DIMENSION(2, 3) :: bo
233 : REAL(KIND=dp) :: epsilon_rho, sx
234 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
235 166 : POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
236 996 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: rho, rho_1_3
237 : TYPE(xc_derivative_type), POINTER :: deriv
238 :
239 166 : CALL timeset(routineN, handle)
240 166 : NULLIFY (deriv)
241 498 : DO i = 1, 2
242 498 : NULLIFY (rho(i)%array, rho_1_3(i)%array)
243 : END DO
244 :
245 166 : CALL section_vals_val_get(xa_params, "scale_x", r_val=sx)
246 :
247 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
248 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
249 : rhob=rho(2)%array, rho_cutoff=epsilon_rho, &
250 166 : local_bounds=bo)
251 166 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
252 166 : CALL xalpha_init(epsilon_rho, xa_parameter)
253 :
254 498 : DO ispin = 1, 2
255 332 : IF (order >= 0) THEN
256 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
257 332 : allocate_deriv=.TRUE.)
258 332 : CALL xc_derivative_get(deriv, deriv_data=e_0)
259 :
260 : CALL xalpha_lsd_0(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
261 332 : e_0, sx)
262 : END IF
263 332 : IF (order >= 1 .OR. order == -1) THEN
264 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
265 664 : allocate_deriv=.TRUE.)
266 332 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
267 :
268 : CALL xalpha_lsd_1(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
269 332 : e_rho, sx)
270 : END IF
271 332 : IF (order >= 2 .OR. order == -2) THEN
272 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
273 192 : rho_spin_name(ispin)], allocate_deriv=.TRUE.)
274 64 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
275 :
276 : CALL xalpha_lsd_2(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
277 64 : e_rho_rho, sx)
278 : END IF
279 332 : IF (order >= 3 .OR. order == -3) THEN
280 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
281 : rho_spin_name(ispin), rho_spin_name(ispin)], &
282 0 : allocate_deriv=.TRUE.)
283 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
284 :
285 : CALL xalpha_lsd_3(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
286 0 : e_rho_rho_rho, sx)
287 : END IF
288 498 : IF (order > 3 .OR. order < -3) THEN
289 0 : CPABORT("derivatives bigger than 3 not implemented")
290 : END IF
291 : END DO
292 166 : CALL timestop(handle)
293 166 : END SUBROUTINE xalpha_lsd_eval
294 :
295 : ! **************************************************************************************************
296 : !> \brief ...
297 : !> \param n ...
298 : !> \param rho ...
299 : !> \param r13 ...
300 : !> \param pot ...
301 : !> \param sx ...
302 : ! **************************************************************************************************
303 2633 : SUBROUTINE xalpha_lda_0(n, rho, r13, pot, sx)
304 :
305 : INTEGER, INTENT(IN) :: n
306 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
307 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
308 : REAL(KIND=dp), INTENT(IN) :: sx
309 :
310 : INTEGER :: ip
311 : REAL(KIND=dp) :: f
312 :
313 2633 : f = sx*flda
314 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE) &
315 2633 : !$OMP SHARED(n,rho,eps_rho,pot,f,r13)
316 :
317 : DO ip = 1, n
318 : IF (rho(ip) > eps_rho) THEN
319 : pot(ip) = pot(ip) + f*r13(ip)*rho(ip)
320 : END IF
321 : END DO
322 :
323 2633 : END SUBROUTINE xalpha_lda_0
324 :
325 : ! **************************************************************************************************
326 : !> \brief ...
327 : !> \param n ...
328 : !> \param rho ...
329 : !> \param r13 ...
330 : !> \param pot ...
331 : !> \param sx ...
332 : ! **************************************************************************************************
333 2623 : SUBROUTINE xalpha_lda_1(n, rho, r13, pot, sx)
334 :
335 : INTEGER, INTENT(IN) :: n
336 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
337 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
338 : REAL(KIND=dp) :: sx
339 :
340 : INTEGER :: ip
341 : REAL(KIND=dp) :: f
342 :
343 2623 : f = f43*flda*sx
344 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
345 2623 : !$OMP SHARED(n,rho,eps_rho,pot,f,r13)
346 : DO ip = 1, n
347 : IF (rho(ip) > eps_rho) THEN
348 : pot(ip) = pot(ip) + f*r13(ip)
349 : END IF
350 : END DO
351 :
352 2623 : END SUBROUTINE xalpha_lda_1
353 :
354 : ! **************************************************************************************************
355 : !> \brief ...
356 : !> \param n ...
357 : !> \param rho ...
358 : !> \param r13 ...
359 : !> \param pot ...
360 : !> \param sx ...
361 : ! **************************************************************************************************
362 380 : SUBROUTINE xalpha_lda_2(n, rho, r13, pot, sx)
363 :
364 : INTEGER, INTENT(IN) :: n
365 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
366 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
367 : REAL(KIND=dp) :: sx
368 :
369 : INTEGER :: ip
370 : REAL(KIND=dp) :: f
371 :
372 380 : f = f13*f43*flda*sx
373 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
374 380 : !$OMP SHARED(n,rho,eps_rho,pot,f,r13)
375 : DO ip = 1, n
376 : IF (rho(ip) > eps_rho) THEN
377 : pot(ip) = pot(ip) + f*r13(ip)/rho(ip)
378 : END IF
379 : END DO
380 :
381 380 : END SUBROUTINE xalpha_lda_2
382 :
383 : ! **************************************************************************************************
384 : !> \brief ...
385 : !> \param n ...
386 : !> \param rho ...
387 : !> \param r13 ...
388 : !> \param pot ...
389 : !> \param sx ...
390 : ! **************************************************************************************************
391 0 : SUBROUTINE xalpha_lda_3(n, rho, r13, pot, sx)
392 :
393 : INTEGER, INTENT(IN) :: n
394 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
395 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
396 : REAL(KIND=dp) :: sx
397 :
398 : INTEGER :: ip
399 : REAL(KIND=dp) :: f
400 :
401 0 : f = -f23*f13*f43*flda*sx
402 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
403 0 : !$OMP SHARED(n,rho,eps_rho,pot,f,r13)
404 : DO ip = 1, n
405 : IF (rho(ip) > eps_rho) THEN
406 : pot(ip) = pot(ip) + f*r13(ip)/(rho(ip)*rho(ip))
407 : END IF
408 : END DO
409 :
410 0 : END SUBROUTINE xalpha_lda_3
411 :
412 : ! **************************************************************************************************
413 : !> \brief ...
414 : !> \param n ...
415 : !> \param rhoa ...
416 : !> \param r13a ...
417 : !> \param pot ...
418 : !> \param sx ...
419 : ! **************************************************************************************************
420 332 : SUBROUTINE xalpha_lsd_0(n, rhoa, r13a, pot, sx)
421 :
422 : INTEGER, INTENT(IN) :: n
423 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
424 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
425 : REAL(KIND=dp) :: sx
426 :
427 : INTEGER :: ip
428 : REAL(KIND=dp) :: f
429 :
430 : ! number of points in array
431 :
432 332 : f = sx*flsd
433 :
434 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
435 332 : !$OMP SHARED(n,rhoa,eps_rho,pot,f,r13a)
436 : DO ip = 1, n
437 :
438 : IF (rhoa(ip) > eps_rho) THEN
439 : pot(ip) = pot(ip) + f*r13a(ip)*rhoa(ip)
440 : END IF
441 :
442 : END DO
443 :
444 332 : END SUBROUTINE xalpha_lsd_0
445 :
446 : ! **************************************************************************************************
447 : !> \brief ...
448 : !> \param n ...
449 : !> \param rhoa ...
450 : !> \param r13a ...
451 : !> \param pota ...
452 : !> \param sx ...
453 : ! **************************************************************************************************
454 332 : SUBROUTINE xalpha_lsd_1(n, rhoa, r13a, pota, sx)
455 :
456 : INTEGER, INTENT(IN) :: n
457 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
458 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pota
459 : REAL(KIND=dp) :: sx
460 :
461 : INTEGER :: ip
462 : REAL(KIND=dp) :: f
463 :
464 : ! number of points in array
465 :
466 332 : f = f43*flsd*sx
467 :
468 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
469 332 : !$OMP SHARED(n,rhoa,eps_rho,pota,f,r13a)
470 : DO ip = 1, n
471 :
472 : IF (rhoa(ip) > eps_rho) THEN
473 : pota(ip) = pota(ip) + f*r13a(ip)
474 : END IF
475 :
476 : END DO
477 :
478 332 : END SUBROUTINE xalpha_lsd_1
479 :
480 : ! **************************************************************************************************
481 : !> \brief ...
482 : !> \param n ...
483 : !> \param rhoa ...
484 : !> \param r13a ...
485 : !> \param potaa ...
486 : !> \param sx ...
487 : ! **************************************************************************************************
488 64 : SUBROUTINE xalpha_lsd_2(n, rhoa, r13a, potaa, sx)
489 :
490 : INTEGER, INTENT(IN) :: n
491 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
492 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: potaa
493 : REAL(KIND=dp) :: sx
494 :
495 : INTEGER :: ip
496 : REAL(KIND=dp) :: f
497 :
498 : ! number of points in array
499 :
500 64 : f = f13*f43*flsd*sx
501 :
502 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
503 64 : !$OMP SHARED(n,rhoa,eps_rho,potaa,f,r13a)
504 : DO ip = 1, n
505 :
506 : IF (rhoa(ip) > eps_rho) THEN
507 : potaa(ip) = potaa(ip) + f*r13a(ip)/rhoa(ip)
508 : END IF
509 :
510 : END DO
511 :
512 64 : END SUBROUTINE xalpha_lsd_2
513 :
514 : ! **************************************************************************************************
515 : !> \brief ...
516 : !> \param n ...
517 : !> \param rhoa ...
518 : !> \param r13a ...
519 : !> \param potaaa ...
520 : !> \param sx ...
521 : ! **************************************************************************************************
522 0 : SUBROUTINE xalpha_lsd_3(n, rhoa, r13a, potaaa, sx)
523 :
524 : INTEGER, INTENT(IN) :: n
525 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
526 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: potaaa
527 : REAL(KIND=dp) :: sx
528 :
529 : INTEGER :: ip
530 : REAL(KIND=dp) :: f
531 :
532 : ! number of points in array
533 :
534 0 : f = -f23*f13*f43*flsd*sx
535 :
536 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
537 0 : !$OMP SHARED(n,rhoa,eps_rho,potaaa,f,r13a)
538 : DO ip = 1, n
539 :
540 : IF (rhoa(ip) > eps_rho) THEN
541 : potaaa(ip) = potaaa(ip) + f*r13a(ip)/(rhoa(ip)*rhoa(ip))
542 : END IF
543 :
544 : END DO
545 :
546 0 : END SUBROUTINE xalpha_lsd_3
547 :
548 : ! **************************************************************************************************
549 : !> \brief ...
550 : !> \param rho_a ...
551 : !> \param rho_b ...
552 : !> \param fxc_aa ...
553 : !> \param fxc_bb ...
554 : !> \param scale_x ...
555 : !> \param eps_rho ...
556 : ! **************************************************************************************************
557 10 : SUBROUTINE xalpha_fxc_eval(rho_a, rho_b, fxc_aa, fxc_bb, scale_x, eps_rho)
558 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_a, rho_b
559 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: fxc_aa, fxc_bb
560 : REAL(KIND=dp), INTENT(IN) :: scale_x, eps_rho
561 :
562 : INTEGER :: i, j, k
563 : INTEGER, DIMENSION(2, 3) :: bo
564 : REAL(KIND=dp) :: eaa, ebb, f, flda, flsd, rhoa, rhob
565 :
566 10 : flda = -3.0_dp/4.0_dp*(3.0_dp/pi)**f13
567 10 : flsd = flda*2.0_dp**f13
568 10 : f = f13*f43*flsd*scale_x
569 100 : bo(1:2, 1:3) = rho_a%pw_grid%bounds_local(1:2, 1:3)
570 : !$OMP PARALLEL DO PRIVATE(i,j,k,rhoa,rhob,eaa,ebb) DEFAULT(NONE)&
571 10 : !$OMP SHARED(bo,rho_a,rho_b,fxc_aa,fxc_bb,f,eps_rho)
572 : DO k = bo(1, 3), bo(2, 3)
573 : DO j = bo(1, 2), bo(2, 2)
574 : DO i = bo(1, 1), bo(2, 1)
575 :
576 : rhoa = rho_a%array(i, j, k)
577 : IF (rhoa > eps_rho) THEN
578 : eaa = rhoa**(-f23)
579 : fxc_aa%array(i, j, k) = fxc_aa%array(i, j, k) + f*eaa
580 : END IF
581 : rhob = rho_b%array(i, j, k)
582 : IF (rhob > eps_rho) THEN
583 : ebb = rhob**(-f23)
584 : fxc_bb%array(i, j, k) = fxc_bb%array(i, j, k) + f*ebb
585 : END IF
586 :
587 : END DO
588 : END DO
589 : END DO
590 :
591 10 : END SUBROUTINE xalpha_fxc_eval
592 :
593 : END MODULE xc_xalpha
594 :
|
