Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate the Thomas-Fermi kinetic energy 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 : !> \author JGH (18.02.2002)
17 : ! **************************************************************************************************
18 : MODULE xc_thomas_fermi
19 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
20 : USE kinds, ONLY: dp
21 : USE xc_derivative_desc, ONLY: deriv_rho,&
22 : deriv_rhoa,&
23 : deriv_rhob
24 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
25 : xc_dset_get_derivative
26 : USE xc_derivative_types, ONLY: xc_derivative_get,&
27 : xc_derivative_type
28 : USE xc_functionals_utilities, ONLY: set_util
29 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
30 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
31 : xc_rho_set_type
32 : #include "../base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 :
36 : PRIVATE
37 :
38 : ! *** Global parameters ***
39 :
40 : REAL(KIND=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp
41 : REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
42 : f23 = 2.0_dp*f13, &
43 : f43 = 4.0_dp*f13, &
44 : f53 = 5.0_dp*f13
45 :
46 : PUBLIC :: thomas_fermi_info, thomas_fermi_lda_eval, thomas_fermi_lsd_eval
47 :
48 : REAL(KIND=dp) :: cf, flda, flsd
49 : REAL(KIND=dp) :: eps_rho
50 :
51 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_thomas_fermi'
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief ...
57 : !> \param cutoff ...
58 : ! **************************************************************************************************
59 216 : SUBROUTINE thomas_fermi_init(cutoff)
60 :
61 : REAL(KIND=dp), INTENT(IN) :: cutoff
62 :
63 216 : eps_rho = cutoff
64 216 : CALL set_util(cutoff)
65 :
66 216 : cf = 0.3_dp*(3.0_dp*pi*pi)**f23
67 216 : flda = cf
68 216 : flsd = flda*2.0_dp**f23
69 :
70 216 : END SUBROUTINE thomas_fermi_init
71 :
72 : ! **************************************************************************************************
73 : !> \brief ...
74 : !> \param lsd ...
75 : !> \param reference ...
76 : !> \param shortform ...
77 : !> \param needs ...
78 : !> \param max_deriv ...
79 : ! **************************************************************************************************
80 216 : SUBROUTINE thomas_fermi_info(lsd, reference, shortform, needs, max_deriv)
81 : LOGICAL, INTENT(in) :: lsd
82 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
83 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
84 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
85 :
86 216 : IF (PRESENT(reference)) THEN
87 0 : reference = "Thomas-Fermi kinetic energy functional: see Parr and Yang"
88 0 : IF (.NOT. lsd) THEN
89 0 : IF (LEN_TRIM(reference) + 6 < LEN(reference)) THEN
90 0 : reference(LEN_TRIM(reference):LEN_TRIM(reference) + 6) = ' {LDA}'
91 : END IF
92 : END IF
93 : END IF
94 216 : IF (PRESENT(shortform)) THEN
95 0 : shortform = "Thomas-Fermi kinetic energy functional"
96 0 : IF (.NOT. lsd) THEN
97 0 : IF (LEN_TRIM(shortform) + 6 < LEN(shortform)) THEN
98 0 : shortform(LEN_TRIM(shortform):LEN_TRIM(shortform) + 6) = ' {LDA}'
99 : END IF
100 : END IF
101 : END IF
102 216 : IF (PRESENT(needs)) THEN
103 216 : IF (lsd) THEN
104 0 : needs%rho_spin = .TRUE.
105 0 : needs%rho_spin_1_3 = .TRUE.
106 : ELSE
107 216 : needs%rho = .TRUE.
108 216 : needs%rho_1_3 = .TRUE.
109 : END IF
110 : END IF
111 216 : IF (PRESENT(max_deriv)) max_deriv = 3
112 :
113 216 : END SUBROUTINE thomas_fermi_info
114 :
115 : ! **************************************************************************************************
116 : !> \brief ...
117 : !> \param rho_set ...
118 : !> \param deriv_set ...
119 : !> \param order ...
120 : ! **************************************************************************************************
121 432 : SUBROUTINE thomas_fermi_lda_eval(rho_set, deriv_set, order)
122 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
123 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
124 : INTEGER, INTENT(in) :: order
125 :
126 : CHARACTER(len=*), PARAMETER :: routineN = 'thomas_fermi_lda_eval'
127 :
128 : INTEGER :: handle, npoints
129 : INTEGER, DIMENSION(2, 3) :: bo
130 : REAL(KIND=dp) :: epsilon_rho
131 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
132 216 : POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho, &
133 216 : r13, rho
134 : TYPE(xc_derivative_type), POINTER :: deriv
135 :
136 216 : CALL timeset(routineN, handle)
137 :
138 : CALL xc_rho_set_get(rho_set, rho_1_3=r13, rho=rho, &
139 216 : local_bounds=bo, rho_cutoff=epsilon_rho)
140 216 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
141 216 : CALL thomas_fermi_init(epsilon_rho)
142 :
143 216 : IF (order >= 0) THEN
144 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
145 216 : allocate_deriv=.TRUE.)
146 216 : CALL xc_derivative_get(deriv, deriv_data=e_0)
147 :
148 216 : CALL thomas_fermi_lda_0(rho, r13, e_0, npoints)
149 : END IF
150 216 : IF (order >= 1 .OR. order == -1) THEN
151 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
152 216 : allocate_deriv=.TRUE.)
153 216 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
154 :
155 216 : CALL thomas_fermi_lda_1(rho, r13, e_rho, npoints)
156 : END IF
157 216 : IF (order >= 2 .OR. order == -2) THEN
158 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
159 0 : allocate_deriv=.TRUE.)
160 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
161 :
162 0 : CALL thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
163 : END IF
164 216 : IF (order >= 3 .OR. order == -3) THEN
165 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
166 0 : allocate_deriv=.TRUE.)
167 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
168 :
169 0 : CALL thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
170 : END IF
171 216 : IF (order > 3 .OR. order < -3) THEN
172 0 : CPABORT("derivatives bigger than 3 not implemented")
173 : END IF
174 216 : CALL timestop(handle)
175 216 : END SUBROUTINE thomas_fermi_lda_eval
176 :
177 : ! **************************************************************************************************
178 : !> \brief ...
179 : !> \param rho_set ...
180 : !> \param deriv_set ...
181 : !> \param order ...
182 : ! **************************************************************************************************
183 0 : SUBROUTINE thomas_fermi_lsd_eval(rho_set, deriv_set, order)
184 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
185 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
186 : INTEGER, INTENT(in) :: order
187 :
188 : CHARACTER(len=*), PARAMETER :: routineN = 'thomas_fermi_lsd_eval'
189 : INTEGER, DIMENSION(2), PARAMETER :: rho_spin_name = [deriv_rhoa, deriv_rhob]
190 :
191 : INTEGER :: handle, i, ispin, npoints
192 : INTEGER, DIMENSION(2, 3) :: bo
193 : REAL(KIND=dp) :: epsilon_rho
194 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
195 0 : POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
196 0 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: rho, rho_1_3
197 : TYPE(xc_derivative_type), POINTER :: deriv
198 :
199 0 : CALL timeset(routineN, handle)
200 0 : NULLIFY (deriv)
201 0 : DO i = 1, 2
202 0 : NULLIFY (rho(i)%array, rho_1_3(i)%array)
203 : END DO
204 :
205 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
206 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
207 : rhob=rho(2)%array, &
208 : rho_cutoff=epsilon_rho, &
209 0 : local_bounds=bo)
210 0 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
211 0 : CALL thomas_fermi_init(epsilon_rho)
212 :
213 0 : DO ispin = 1, 2
214 0 : IF (order >= 0) THEN
215 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
216 0 : allocate_deriv=.TRUE.)
217 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
218 :
219 : CALL thomas_fermi_lsd_0(rho(ispin)%array, rho_1_3(ispin)%array, &
220 0 : e_0, npoints)
221 : END IF
222 0 : IF (order >= 1 .OR. order == -1) THEN
223 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
224 0 : allocate_deriv=.TRUE.)
225 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
226 :
227 : CALL thomas_fermi_lsd_1(rho(ispin)%array, rho_1_3(ispin)%array, &
228 0 : e_rho, npoints)
229 : END IF
230 0 : IF (order >= 2 .OR. order == -2) THEN
231 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
232 0 : rho_spin_name(ispin)], allocate_deriv=.TRUE.)
233 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
234 :
235 : CALL thomas_fermi_lsd_2(rho(ispin)%array, rho_1_3(ispin)%array, &
236 0 : e_rho_rho, npoints)
237 : END IF
238 0 : IF (order >= 3 .OR. order == -3) THEN
239 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
240 : rho_spin_name(ispin), rho_spin_name(ispin)], &
241 0 : allocate_deriv=.TRUE.)
242 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
243 :
244 : CALL thomas_fermi_lsd_3(rho(ispin)%array, rho_1_3(ispin)%array, &
245 0 : e_rho_rho_rho, npoints)
246 : END IF
247 0 : IF (order > 3 .OR. order < -3) THEN
248 0 : CPABORT("derivatives bigger than 3 not implemented")
249 : END IF
250 : END DO
251 0 : CALL timestop(handle)
252 0 : END SUBROUTINE thomas_fermi_lsd_eval
253 :
254 : ! **************************************************************************************************
255 : !> \brief ...
256 : !> \param rho ...
257 : !> \param r13 ...
258 : !> \param e_0 ...
259 : !> \param npoints ...
260 : ! **************************************************************************************************
261 216 : SUBROUTINE thomas_fermi_lda_0(rho, r13, e_0, npoints)
262 :
263 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
264 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
265 : INTEGER, INTENT(in) :: npoints
266 :
267 : INTEGER :: ip
268 :
269 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
270 216 : !$OMP SHARED(npoints,rho,eps_rho,e_0,flda,r13)
271 : DO ip = 1, npoints
272 :
273 : IF (rho(ip) > eps_rho) THEN
274 :
275 : e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip)
276 :
277 : END IF
278 :
279 : END DO
280 :
281 216 : END SUBROUTINE thomas_fermi_lda_0
282 :
283 : ! **************************************************************************************************
284 : !> \brief ...
285 : !> \param rho ...
286 : !> \param r13 ...
287 : !> \param e_rho ...
288 : !> \param npoints ...
289 : ! **************************************************************************************************
290 216 : SUBROUTINE thomas_fermi_lda_1(rho, r13, e_rho, npoints)
291 :
292 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
293 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho
294 : INTEGER, INTENT(in) :: npoints
295 :
296 : INTEGER :: ip
297 : REAL(KIND=dp) :: f
298 :
299 216 : f = f53*flda
300 :
301 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE) &
302 216 : !$OMP SHARED(npoints,rho,eps_rho,e_rho,f,r13)
303 : DO ip = 1, npoints
304 :
305 : IF (rho(ip) > eps_rho) THEN
306 :
307 : e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip)
308 :
309 : END IF
310 :
311 : END DO
312 :
313 216 : END SUBROUTINE thomas_fermi_lda_1
314 :
315 : ! **************************************************************************************************
316 : !> \brief ...
317 : !> \param rho ...
318 : !> \param r13 ...
319 : !> \param e_rho_rho ...
320 : !> \param npoints ...
321 : ! **************************************************************************************************
322 0 : SUBROUTINE thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
323 :
324 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
325 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho
326 : INTEGER, INTENT(in) :: npoints
327 :
328 : INTEGER :: ip
329 : REAL(KIND=dp) :: f
330 :
331 0 : f = f23*f53*flda
332 :
333 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
334 0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho,f,r13)
335 : DO ip = 1, npoints
336 :
337 : IF (rho(ip) > eps_rho) THEN
338 :
339 : e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip)
340 :
341 : END IF
342 :
343 : END DO
344 :
345 0 : END SUBROUTINE thomas_fermi_lda_2
346 :
347 : ! **************************************************************************************************
348 : !> \brief ...
349 : !> \param rho ...
350 : !> \param r13 ...
351 : !> \param e_rho_rho_rho ...
352 : !> \param npoints ...
353 : ! **************************************************************************************************
354 0 : SUBROUTINE thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
355 :
356 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
357 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho
358 : INTEGER, INTENT(in) :: npoints
359 :
360 : INTEGER :: ip
361 : REAL(KIND=dp) :: f
362 :
363 0 : f = -f13*f23*f53*flda
364 :
365 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
366 0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho_rho,f,r13)
367 : DO ip = 1, npoints
368 :
369 : IF (rho(ip) > eps_rho) THEN
370 :
371 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip))
372 :
373 : END IF
374 :
375 : END DO
376 :
377 0 : END SUBROUTINE thomas_fermi_lda_3
378 :
379 : ! **************************************************************************************************
380 : !> \brief ...
381 : !> \param rhoa ...
382 : !> \param r13a ...
383 : !> \param e_0 ...
384 : !> \param npoints ...
385 : ! **************************************************************************************************
386 0 : SUBROUTINE thomas_fermi_lsd_0(rhoa, r13a, e_0, npoints)
387 :
388 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
389 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
390 : INTEGER, INTENT(in) :: npoints
391 :
392 : INTEGER :: ip
393 :
394 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
395 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_0,flsd,r13a)
396 : DO ip = 1, npoints
397 :
398 : IF (rhoa(ip) > eps_rho) THEN
399 : e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip)
400 : END IF
401 :
402 : END DO
403 :
404 0 : END SUBROUTINE thomas_fermi_lsd_0
405 :
406 : ! **************************************************************************************************
407 : !> \brief ...
408 : !> \param rhoa ...
409 : !> \param r13a ...
410 : !> \param e_rho ...
411 : !> \param npoints ...
412 : ! **************************************************************************************************
413 0 : SUBROUTINE thomas_fermi_lsd_1(rhoa, r13a, e_rho, npoints)
414 :
415 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
416 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho
417 : INTEGER, INTENT(in) :: npoints
418 :
419 : INTEGER :: ip
420 : REAL(KIND=dp) :: f
421 :
422 0 : f = f53*flsd
423 :
424 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
425 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho,f,r13a)
426 : DO ip = 1, npoints
427 :
428 : IF (rhoa(ip) > eps_rho) THEN
429 : e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip)
430 : END IF
431 :
432 : END DO
433 :
434 0 : END SUBROUTINE thomas_fermi_lsd_1
435 :
436 : ! **************************************************************************************************
437 : !> \brief ...
438 : !> \param rhoa ...
439 : !> \param r13a ...
440 : !> \param e_rho_rho ...
441 : !> \param npoints ...
442 : ! **************************************************************************************************
443 0 : SUBROUTINE thomas_fermi_lsd_2(rhoa, r13a, e_rho_rho, npoints)
444 :
445 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
446 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho
447 : INTEGER, INTENT(in) :: npoints
448 :
449 : INTEGER :: ip
450 : REAL(KIND=dp) :: f
451 :
452 0 : f = f23*f53*flsd
453 :
454 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
455 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho,f,r13a)
456 :
457 : DO ip = 1, npoints
458 :
459 : IF (rhoa(ip) > eps_rho) THEN
460 : e_rho_rho(ip) = e_rho_rho(ip) + f/r13a(ip)
461 : END IF
462 :
463 : END DO
464 :
465 0 : END SUBROUTINE thomas_fermi_lsd_2
466 :
467 : ! **************************************************************************************************
468 : !> \brief ...
469 : !> \param rhoa ...
470 : !> \param r13a ...
471 : !> \param e_rho_rho_rho ...
472 : !> \param npoints ...
473 : ! **************************************************************************************************
474 0 : SUBROUTINE thomas_fermi_lsd_3(rhoa, r13a, e_rho_rho_rho, npoints)
475 :
476 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
477 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho
478 : INTEGER, INTENT(in) :: npoints
479 :
480 : INTEGER :: ip
481 : REAL(KIND=dp) :: f
482 :
483 0 : f = -f13*f23*f53*flsd
484 :
485 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
486 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho_rho,f,r13a)
487 : DO ip = 1, npoints
488 :
489 : IF (rhoa(ip) > eps_rho) THEN
490 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13a(ip)*rhoa(ip))
491 : END IF
492 :
493 : END DO
494 :
495 0 : END SUBROUTINE thomas_fermi_lsd_3
496 :
497 : END MODULE xc_thomas_fermi
498 :
|
