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