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 2731 : SUBROUTINE xalpha_init(cutoff, xalpha)
62 :
63 : REAL(KIND=dp), INTENT(IN) :: cutoff
64 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: xalpha
65 :
66 2731 : eps_rho = cutoff
67 2731 : CALL set_util(cutoff)
68 2731 : IF (PRESENT(xalpha)) THEN
69 2731 : xparam = xalpha
70 : ELSE
71 0 : xparam = 2.0_dp/3.0_dp
72 : END IF
73 :
74 2731 : flda = -9.0_dp/8.0_dp*xparam*(3.0_dp/pi)**f13
75 2731 : flsd = flda*2.0_dp**f13
76 :
77 2731 : 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 3051 : 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 3051 : my_xparam = 2.0_dp/3.0_dp
100 3051 : IF (PRESENT(xa_parameter)) my_xparam = xa_parameter
101 3051 : my_scaling = 1.0_dp
102 3051 : IF (PRESENT(scaling)) my_scaling = scaling
103 :
104 3051 : IF (PRESENT(reference)) THEN
105 38 : 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 38 : "Dirac/Slater local exchange; parameter=", my_xparam
111 : END IF
112 38 : IF (.NOT. lsd) THEN
113 33 : IF (LEN_TRIM(reference) + 6 < LEN(reference)) THEN
114 33 : reference(LEN_TRIM(reference):LEN_TRIM(reference) + 6) = ' {LDA}'
115 : END IF
116 : END IF
117 : END IF
118 3051 : IF (PRESENT(shortform)) THEN
119 38 : 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 38 : WRITE (shortform, '(A,F8.4)') "Dirac/Slater exchange", my_xparam
123 : END IF
124 38 : IF (.NOT. lsd) THEN
125 33 : IF (LEN_TRIM(shortform) + 6 < LEN(shortform)) THEN
126 33 : shortform(LEN_TRIM(shortform):LEN_TRIM(shortform) + 6) = ' {LDA}'
127 : END IF
128 : END IF
129 : END IF
130 3051 : IF (PRESENT(needs)) THEN
131 3013 : IF (lsd) THEN
132 180 : needs%rho_spin = .TRUE.
133 180 : needs%rho_spin_1_3 = .TRUE.
134 : ELSE
135 2833 : needs%rho = .TRUE.
136 2833 : needs%rho_1_3 = .TRUE.
137 : END IF
138 : END IF
139 3051 : IF (PRESENT(max_deriv)) max_deriv = 3
140 :
141 3051 : 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 7683 : 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 2561 : POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho, &
165 2561 : r13, rho
166 : TYPE(xc_derivative_type), POINTER :: deriv
167 :
168 2561 : CALL timeset(routineN, handle)
169 :
170 2561 : 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 2561 : local_bounds=bo, rho_cutoff=epsilon_rho)
174 2561 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
175 2561 : CALL xalpha_init(epsilon_rho, xa_parameter)
176 :
177 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
178 2561 : allocate_deriv=.TRUE.)
179 2561 : CALL xc_derivative_get(deriv, deriv_data=e_0)
180 :
181 2561 : CALL xalpha_lda_0(npoints, rho, r13, e_0, sx)
182 :
183 2561 : IF (order >= 1) THEN
184 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
185 2545 : allocate_deriv=.TRUE.)
186 2545 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
187 :
188 2545 : CALL xalpha_lda_1(npoints, rho, r13, e_rho, sx)
189 : END IF
190 2561 : IF (order >= 2) THEN
191 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
192 290 : allocate_deriv=.TRUE.)
193 290 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
194 :
195 290 : CALL xalpha_lda_2(npoints, rho, r13, e_rho_rho, sx)
196 : END IF
197 2561 : IF (order >= 3) THEN
198 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
199 0 : allocate_deriv=.TRUE.)
200 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
201 :
202 0 : CALL xalpha_lda_3(npoints, rho, r13, e_rho_rho_rho, sx)
203 : END IF
204 2561 : IF (order > 3) THEN
205 0 : CPABORT("derivatives bigger than 3 not implemented")
206 : END IF
207 2561 : CALL timestop(handle)
208 :
209 2561 : END SUBROUTINE xalpha_lda_eval
210 :
211 : ! **************************************************************************************************
212 : !> \brief ...
213 : !> \param rho_set ...
214 : !> \param deriv_set ...
215 : !> \param order ...
216 : !> \param xa_params ...
217 : !> \param xa_parameter ...
218 : ! **************************************************************************************************
219 680 : SUBROUTINE xalpha_lsd_eval(rho_set, deriv_set, order, xa_params, xa_parameter)
220 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
221 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
222 : INTEGER, INTENT(in) :: order
223 : TYPE(section_vals_type), POINTER :: xa_params
224 : REAL(KIND=dp), INTENT(in), OPTIONAL :: xa_parameter
225 :
226 : CHARACTER(len=*), PARAMETER :: routineN = 'xalpha_lsd_eval'
227 : INTEGER, DIMENSION(2), PARAMETER :: rho_spin_name = [deriv_rhoa, deriv_rhob]
228 :
229 : INTEGER :: handle, i, ispin, npoints
230 : INTEGER, DIMENSION(2, 3) :: bo
231 : REAL(KIND=dp) :: epsilon_rho, sx
232 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
233 170 : POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
234 1020 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: rho, rho_1_3
235 : TYPE(xc_derivative_type), POINTER :: deriv
236 :
237 170 : CALL timeset(routineN, handle)
238 170 : NULLIFY (deriv)
239 510 : DO i = 1, 2
240 510 : NULLIFY (rho(i)%array, rho_1_3(i)%array)
241 : END DO
242 :
243 170 : CALL section_vals_val_get(xa_params, "scale_x", r_val=sx)
244 :
245 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
246 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
247 : rhob=rho(2)%array, rho_cutoff=epsilon_rho, &
248 170 : local_bounds=bo)
249 170 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
250 170 : CALL xalpha_init(epsilon_rho, xa_parameter)
251 :
252 510 : DO ispin = 1, 2
253 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
254 340 : allocate_deriv=.TRUE.)
255 340 : CALL xc_derivative_get(deriv, deriv_data=e_0)
256 :
257 : CALL xalpha_lsd_0(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
258 340 : e_0, sx)
259 340 : IF (order >= 1) THEN
260 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
261 680 : allocate_deriv=.TRUE.)
262 340 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
263 :
264 : CALL xalpha_lsd_1(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
265 340 : e_rho, sx)
266 : END IF
267 340 : IF (order >= 2) THEN
268 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
269 192 : rho_spin_name(ispin)], allocate_deriv=.TRUE.)
270 64 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
271 :
272 : CALL xalpha_lsd_2(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
273 64 : e_rho_rho, sx)
274 : END IF
275 340 : IF (order >= 3) THEN
276 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
277 : rho_spin_name(ispin), rho_spin_name(ispin)], &
278 0 : allocate_deriv=.TRUE.)
279 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
280 :
281 : CALL xalpha_lsd_3(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
282 0 : e_rho_rho_rho, sx)
283 : END IF
284 510 : IF (order > 3) THEN
285 0 : CPABORT("derivatives bigger than 3 not implemented")
286 : END IF
287 : END DO
288 170 : CALL timestop(handle)
289 170 : END SUBROUTINE xalpha_lsd_eval
290 :
291 : ! **************************************************************************************************
292 : !> \brief ...
293 : !> \param n ...
294 : !> \param rho ...
295 : !> \param r13 ...
296 : !> \param pot ...
297 : !> \param sx ...
298 : ! **************************************************************************************************
299 2561 : SUBROUTINE xalpha_lda_0(n, rho, r13, pot, sx)
300 :
301 : INTEGER, INTENT(IN) :: n
302 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
303 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
304 : REAL(KIND=dp), INTENT(IN) :: sx
305 :
306 : INTEGER :: ip
307 : REAL(KIND=dp) :: f
308 :
309 2561 : f = sx*flda
310 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE) &
311 2561 : !$OMP SHARED(n,rho,eps_rho,pot,f,r13)
312 :
313 : DO ip = 1, n
314 : IF (rho(ip) > eps_rho) THEN
315 : pot(ip) = pot(ip) + f*r13(ip)*rho(ip)
316 : END IF
317 : END DO
318 :
319 2561 : END SUBROUTINE xalpha_lda_0
320 :
321 : ! **************************************************************************************************
322 : !> \brief ...
323 : !> \param n ...
324 : !> \param rho ...
325 : !> \param r13 ...
326 : !> \param pot ...
327 : !> \param sx ...
328 : ! **************************************************************************************************
329 2545 : SUBROUTINE xalpha_lda_1(n, rho, r13, pot, sx)
330 :
331 : INTEGER, INTENT(IN) :: n
332 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
333 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
334 : REAL(KIND=dp) :: sx
335 :
336 : INTEGER :: ip
337 : REAL(KIND=dp) :: f
338 :
339 2545 : f = f43*flda*sx
340 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
341 2545 : !$OMP SHARED(n,rho,eps_rho,pot,f,r13)
342 : DO ip = 1, n
343 : IF (rho(ip) > eps_rho) THEN
344 : pot(ip) = pot(ip) + f*r13(ip)
345 : END IF
346 : END DO
347 :
348 2545 : END SUBROUTINE xalpha_lda_1
349 :
350 : ! **************************************************************************************************
351 : !> \brief ...
352 : !> \param n ...
353 : !> \param rho ...
354 : !> \param r13 ...
355 : !> \param pot ...
356 : !> \param sx ...
357 : ! **************************************************************************************************
358 290 : SUBROUTINE xalpha_lda_2(n, rho, r13, pot, sx)
359 :
360 : INTEGER, INTENT(IN) :: n
361 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
362 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
363 : REAL(KIND=dp) :: sx
364 :
365 : INTEGER :: ip
366 : REAL(KIND=dp) :: f
367 :
368 290 : f = f13*f43*flda*sx
369 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
370 290 : !$OMP SHARED(n,rho,eps_rho,pot,f,r13)
371 : DO ip = 1, n
372 : IF (rho(ip) > eps_rho) THEN
373 : pot(ip) = pot(ip) + f*r13(ip)/rho(ip)
374 : END IF
375 : END DO
376 :
377 290 : END SUBROUTINE xalpha_lda_2
378 :
379 : ! **************************************************************************************************
380 : !> \brief ...
381 : !> \param n ...
382 : !> \param rho ...
383 : !> \param r13 ...
384 : !> \param pot ...
385 : !> \param sx ...
386 : ! **************************************************************************************************
387 0 : SUBROUTINE xalpha_lda_3(n, rho, r13, pot, sx)
388 :
389 : INTEGER, INTENT(IN) :: n
390 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
391 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
392 : REAL(KIND=dp) :: sx
393 :
394 : INTEGER :: ip
395 : REAL(KIND=dp) :: f
396 :
397 0 : f = -f23*f13*f43*flda*sx
398 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
399 0 : !$OMP SHARED(n,rho,eps_rho,pot,f,r13)
400 : DO ip = 1, n
401 : IF (rho(ip) > eps_rho) THEN
402 : pot(ip) = pot(ip) + f*r13(ip)/(rho(ip)*rho(ip))
403 : END IF
404 : END DO
405 :
406 0 : END SUBROUTINE xalpha_lda_3
407 :
408 : ! **************************************************************************************************
409 : !> \brief ...
410 : !> \param n ...
411 : !> \param rhoa ...
412 : !> \param r13a ...
413 : !> \param pot ...
414 : !> \param sx ...
415 : ! **************************************************************************************************
416 340 : SUBROUTINE xalpha_lsd_0(n, rhoa, r13a, pot, sx)
417 :
418 : INTEGER, INTENT(IN) :: n
419 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
420 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
421 : REAL(KIND=dp) :: sx
422 :
423 : INTEGER :: ip
424 : REAL(KIND=dp) :: f
425 :
426 : ! number of points in array
427 :
428 340 : f = sx*flsd
429 :
430 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
431 340 : !$OMP SHARED(n,rhoa,eps_rho,pot,f,r13a)
432 : DO ip = 1, n
433 :
434 : IF (rhoa(ip) > eps_rho) THEN
435 : pot(ip) = pot(ip) + f*r13a(ip)*rhoa(ip)
436 : END IF
437 :
438 : END DO
439 :
440 340 : END SUBROUTINE xalpha_lsd_0
441 :
442 : ! **************************************************************************************************
443 : !> \brief ...
444 : !> \param n ...
445 : !> \param rhoa ...
446 : !> \param r13a ...
447 : !> \param pota ...
448 : !> \param sx ...
449 : ! **************************************************************************************************
450 340 : SUBROUTINE xalpha_lsd_1(n, rhoa, r13a, pota, sx)
451 :
452 : INTEGER, INTENT(IN) :: n
453 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
454 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pota
455 : REAL(KIND=dp) :: sx
456 :
457 : INTEGER :: ip
458 : REAL(KIND=dp) :: f
459 :
460 : ! number of points in array
461 :
462 340 : f = f43*flsd*sx
463 :
464 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
465 340 : !$OMP SHARED(n,rhoa,eps_rho,pota,f,r13a)
466 : DO ip = 1, n
467 :
468 : IF (rhoa(ip) > eps_rho) THEN
469 : pota(ip) = pota(ip) + f*r13a(ip)
470 : END IF
471 :
472 : END DO
473 :
474 340 : END SUBROUTINE xalpha_lsd_1
475 :
476 : ! **************************************************************************************************
477 : !> \brief ...
478 : !> \param n ...
479 : !> \param rhoa ...
480 : !> \param r13a ...
481 : !> \param potaa ...
482 : !> \param sx ...
483 : ! **************************************************************************************************
484 64 : SUBROUTINE xalpha_lsd_2(n, rhoa, r13a, potaa, sx)
485 :
486 : INTEGER, INTENT(IN) :: n
487 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
488 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: potaa
489 : REAL(KIND=dp) :: sx
490 :
491 : INTEGER :: ip
492 : REAL(KIND=dp) :: f
493 :
494 : ! number of points in array
495 :
496 64 : f = f13*f43*flsd*sx
497 :
498 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
499 64 : !$OMP SHARED(n,rhoa,eps_rho,potaa,f,r13a)
500 : DO ip = 1, n
501 :
502 : IF (rhoa(ip) > eps_rho) THEN
503 : potaa(ip) = potaa(ip) + f*r13a(ip)/rhoa(ip)
504 : END IF
505 :
506 : END DO
507 :
508 64 : END SUBROUTINE xalpha_lsd_2
509 :
510 : ! **************************************************************************************************
511 : !> \brief ...
512 : !> \param n ...
513 : !> \param rhoa ...
514 : !> \param r13a ...
515 : !> \param potaaa ...
516 : !> \param sx ...
517 : ! **************************************************************************************************
518 0 : SUBROUTINE xalpha_lsd_3(n, rhoa, r13a, potaaa, sx)
519 :
520 : INTEGER, INTENT(IN) :: n
521 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
522 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: potaaa
523 : REAL(KIND=dp) :: sx
524 :
525 : INTEGER :: ip
526 : REAL(KIND=dp) :: f
527 :
528 : ! number of points in array
529 :
530 0 : f = -f23*f13*f43*flsd*sx
531 :
532 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
533 0 : !$OMP SHARED(n,rhoa,eps_rho,potaaa,f,r13a)
534 : DO ip = 1, n
535 :
536 : IF (rhoa(ip) > eps_rho) THEN
537 : potaaa(ip) = potaaa(ip) + f*r13a(ip)/(rhoa(ip)*rhoa(ip))
538 : END IF
539 :
540 : END DO
541 :
542 0 : END SUBROUTINE xalpha_lsd_3
543 :
544 : ! **************************************************************************************************
545 : !> \brief ...
546 : !> \param rho_a ...
547 : !> \param rho_b ...
548 : !> \param fxc_aa ...
549 : !> \param fxc_bb ...
550 : !> \param scale_x ...
551 : !> \param eps_rho ...
552 : ! **************************************************************************************************
553 10 : SUBROUTINE xalpha_fxc_eval(rho_a, rho_b, fxc_aa, fxc_bb, scale_x, eps_rho)
554 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_a, rho_b
555 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: fxc_aa, fxc_bb
556 : REAL(KIND=dp), INTENT(IN) :: scale_x, eps_rho
557 :
558 : INTEGER :: i, j, k
559 : INTEGER, DIMENSION(2, 3) :: bo
560 : REAL(KIND=dp) :: eaa, ebb, f, flda, flsd, rhoa, rhob
561 :
562 10 : flda = -3.0_dp/4.0_dp*(3.0_dp/pi)**f13
563 10 : flsd = flda*2.0_dp**f13
564 10 : f = f13*f43*flsd*scale_x
565 100 : bo(1:2, 1:3) = rho_a%pw_grid%bounds_local(1:2, 1:3)
566 : !$OMP PARALLEL DO PRIVATE(i,j,k,rhoa,rhob,eaa,ebb) DEFAULT(NONE)&
567 10 : !$OMP SHARED(bo,rho_a,rho_b,fxc_aa,fxc_bb,f,eps_rho)
568 : DO k = bo(1, 3), bo(2, 3)
569 : DO j = bo(1, 2), bo(2, 2)
570 : DO i = bo(1, 1), bo(2, 1)
571 :
572 : rhoa = rho_a%array(i, j, k)
573 : IF (rhoa > eps_rho) THEN
574 : eaa = rhoa**(-f23)
575 : fxc_aa%array(i, j, k) = fxc_aa%array(i, j, k) + f*eaa
576 : END IF
577 : rhob = rho_b%array(i, j, k)
578 : IF (rhob > eps_rho) THEN
579 : ebb = rhob**(-f23)
580 : fxc_bb%array(i, j, k) = fxc_bb%array(i, j, k) + f*ebb
581 : END IF
582 :
583 : END DO
584 : END DO
585 : END DO
586 :
587 10 : END SUBROUTINE xalpha_fxc_eval
588 :
589 : END MODULE xc_xalpha
590 :
|
