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