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 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
143 216 : allocate_deriv=.TRUE.)
144 216 : CALL xc_derivative_get(deriv, deriv_data=e_0)
145 :
146 216 : CALL thomas_fermi_lda_0(rho, r13, e_0, npoints)
147 216 : IF (order >= 1) THEN
148 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
149 216 : allocate_deriv=.TRUE.)
150 216 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
151 :
152 216 : CALL thomas_fermi_lda_1(rho, r13, e_rho, npoints)
153 : END IF
154 216 : IF (order >= 2) THEN
155 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
156 0 : allocate_deriv=.TRUE.)
157 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
158 :
159 0 : CALL thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
160 : END IF
161 216 : IF (order >= 3) THEN
162 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
163 0 : allocate_deriv=.TRUE.)
164 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
165 :
166 0 : CALL thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
167 : END IF
168 216 : IF (order > 3) THEN
169 0 : CPABORT("derivatives bigger than 3 not implemented")
170 : END IF
171 216 : CALL timestop(handle)
172 216 : END SUBROUTINE thomas_fermi_lda_eval
173 :
174 : ! **************************************************************************************************
175 : !> \brief ...
176 : !> \param rho_set ...
177 : !> \param deriv_set ...
178 : !> \param order ...
179 : ! **************************************************************************************************
180 0 : SUBROUTINE thomas_fermi_lsd_eval(rho_set, deriv_set, order)
181 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
182 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
183 : INTEGER, INTENT(in) :: order
184 :
185 : CHARACTER(len=*), PARAMETER :: routineN = 'thomas_fermi_lsd_eval'
186 : INTEGER, DIMENSION(2), PARAMETER :: rho_spin_name = [deriv_rhoa, deriv_rhob]
187 :
188 : INTEGER :: handle, i, ispin, npoints
189 : INTEGER, DIMENSION(2, 3) :: bo
190 : REAL(KIND=dp) :: epsilon_rho
191 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
192 0 : POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
193 0 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: rho, rho_1_3
194 : TYPE(xc_derivative_type), POINTER :: deriv
195 :
196 0 : CALL timeset(routineN, handle)
197 0 : NULLIFY (deriv)
198 0 : DO i = 1, 2
199 0 : NULLIFY (rho(i)%array, rho_1_3(i)%array)
200 : END DO
201 :
202 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
203 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
204 : rhob=rho(2)%array, &
205 : rho_cutoff=epsilon_rho, &
206 0 : local_bounds=bo)
207 0 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
208 0 : CALL thomas_fermi_init(epsilon_rho)
209 :
210 0 : DO ispin = 1, 2
211 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
212 0 : allocate_deriv=.TRUE.)
213 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
214 :
215 : CALL thomas_fermi_lsd_0(rho(ispin)%array, rho_1_3(ispin)%array, &
216 0 : e_0, npoints)
217 0 : IF (order >= 1) THEN
218 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
219 0 : allocate_deriv=.TRUE.)
220 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
221 :
222 : CALL thomas_fermi_lsd_1(rho(ispin)%array, rho_1_3(ispin)%array, &
223 0 : e_rho, npoints)
224 : END IF
225 0 : IF (order >= 2) THEN
226 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
227 0 : rho_spin_name(ispin)], allocate_deriv=.TRUE.)
228 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
229 :
230 : CALL thomas_fermi_lsd_2(rho(ispin)%array, rho_1_3(ispin)%array, &
231 0 : e_rho_rho, npoints)
232 : END IF
233 0 : IF (order >= 3) THEN
234 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
235 : rho_spin_name(ispin), rho_spin_name(ispin)], &
236 0 : allocate_deriv=.TRUE.)
237 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
238 :
239 : CALL thomas_fermi_lsd_3(rho(ispin)%array, rho_1_3(ispin)%array, &
240 0 : e_rho_rho_rho, npoints)
241 : END IF
242 0 : IF (order > 3) THEN
243 0 : CPABORT("derivatives bigger than 3 not implemented")
244 : END IF
245 : END DO
246 0 : CALL timestop(handle)
247 0 : END SUBROUTINE thomas_fermi_lsd_eval
248 :
249 : ! **************************************************************************************************
250 : !> \brief ...
251 : !> \param rho ...
252 : !> \param r13 ...
253 : !> \param e_0 ...
254 : !> \param npoints ...
255 : ! **************************************************************************************************
256 216 : SUBROUTINE thomas_fermi_lda_0(rho, r13, e_0, npoints)
257 :
258 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
259 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
260 : INTEGER, INTENT(in) :: npoints
261 :
262 : INTEGER :: ip
263 :
264 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
265 216 : !$OMP SHARED(npoints,rho,eps_rho,e_0,flda,r13)
266 : DO ip = 1, npoints
267 :
268 : IF (rho(ip) > eps_rho) THEN
269 :
270 : e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip)
271 :
272 : END IF
273 :
274 : END DO
275 :
276 216 : END SUBROUTINE thomas_fermi_lda_0
277 :
278 : ! **************************************************************************************************
279 : !> \brief ...
280 : !> \param rho ...
281 : !> \param r13 ...
282 : !> \param e_rho ...
283 : !> \param npoints ...
284 : ! **************************************************************************************************
285 216 : SUBROUTINE thomas_fermi_lda_1(rho, r13, e_rho, npoints)
286 :
287 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
288 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho
289 : INTEGER, INTENT(in) :: npoints
290 :
291 : INTEGER :: ip
292 : REAL(KIND=dp) :: f
293 :
294 216 : f = f53*flda
295 :
296 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE) &
297 216 : !$OMP SHARED(npoints,rho,eps_rho,e_rho,f,r13)
298 : DO ip = 1, npoints
299 :
300 : IF (rho(ip) > eps_rho) THEN
301 :
302 : e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip)
303 :
304 : END IF
305 :
306 : END DO
307 :
308 216 : END SUBROUTINE thomas_fermi_lda_1
309 :
310 : ! **************************************************************************************************
311 : !> \brief ...
312 : !> \param rho ...
313 : !> \param r13 ...
314 : !> \param e_rho_rho ...
315 : !> \param npoints ...
316 : ! **************************************************************************************************
317 0 : SUBROUTINE thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
318 :
319 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
320 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho
321 : INTEGER, INTENT(in) :: npoints
322 :
323 : INTEGER :: ip
324 : REAL(KIND=dp) :: f
325 :
326 0 : f = f23*f53*flda
327 :
328 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
329 0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho,f,r13)
330 : DO ip = 1, npoints
331 :
332 : IF (rho(ip) > eps_rho) THEN
333 :
334 : e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip)
335 :
336 : END IF
337 :
338 : END DO
339 :
340 0 : END SUBROUTINE thomas_fermi_lda_2
341 :
342 : ! **************************************************************************************************
343 : !> \brief ...
344 : !> \param rho ...
345 : !> \param r13 ...
346 : !> \param e_rho_rho_rho ...
347 : !> \param npoints ...
348 : ! **************************************************************************************************
349 0 : SUBROUTINE thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
350 :
351 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
352 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho
353 : INTEGER, INTENT(in) :: npoints
354 :
355 : INTEGER :: ip
356 : REAL(KIND=dp) :: f
357 :
358 0 : f = -f13*f23*f53*flda
359 :
360 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
361 0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho_rho,f,r13)
362 : DO ip = 1, npoints
363 :
364 : IF (rho(ip) > eps_rho) THEN
365 :
366 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip))
367 :
368 : END IF
369 :
370 : END DO
371 :
372 0 : END SUBROUTINE thomas_fermi_lda_3
373 :
374 : ! **************************************************************************************************
375 : !> \brief ...
376 : !> \param rhoa ...
377 : !> \param r13a ...
378 : !> \param e_0 ...
379 : !> \param npoints ...
380 : ! **************************************************************************************************
381 0 : SUBROUTINE thomas_fermi_lsd_0(rhoa, r13a, e_0, npoints)
382 :
383 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
384 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
385 : INTEGER, INTENT(in) :: npoints
386 :
387 : INTEGER :: ip
388 :
389 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
390 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_0,flsd,r13a)
391 : DO ip = 1, npoints
392 :
393 : IF (rhoa(ip) > eps_rho) THEN
394 : e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip)
395 : END IF
396 :
397 : END DO
398 :
399 0 : END SUBROUTINE thomas_fermi_lsd_0
400 :
401 : ! **************************************************************************************************
402 : !> \brief ...
403 : !> \param rhoa ...
404 : !> \param r13a ...
405 : !> \param e_rho ...
406 : !> \param npoints ...
407 : ! **************************************************************************************************
408 0 : SUBROUTINE thomas_fermi_lsd_1(rhoa, r13a, e_rho, npoints)
409 :
410 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
411 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho
412 : INTEGER, INTENT(in) :: npoints
413 :
414 : INTEGER :: ip
415 : REAL(KIND=dp) :: f
416 :
417 0 : f = f53*flsd
418 :
419 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
420 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho,f,r13a)
421 : DO ip = 1, npoints
422 :
423 : IF (rhoa(ip) > eps_rho) THEN
424 : e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip)
425 : END IF
426 :
427 : END DO
428 :
429 0 : END SUBROUTINE thomas_fermi_lsd_1
430 :
431 : ! **************************************************************************************************
432 : !> \brief ...
433 : !> \param rhoa ...
434 : !> \param r13a ...
435 : !> \param e_rho_rho ...
436 : !> \param npoints ...
437 : ! **************************************************************************************************
438 0 : SUBROUTINE thomas_fermi_lsd_2(rhoa, r13a, e_rho_rho, npoints)
439 :
440 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
441 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho
442 : INTEGER, INTENT(in) :: npoints
443 :
444 : INTEGER :: ip
445 : REAL(KIND=dp) :: f
446 :
447 0 : f = f23*f53*flsd
448 :
449 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
450 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho,f,r13a)
451 :
452 : DO ip = 1, npoints
453 :
454 : IF (rhoa(ip) > eps_rho) THEN
455 : e_rho_rho(ip) = e_rho_rho(ip) + f/r13a(ip)
456 : END IF
457 :
458 : END DO
459 :
460 0 : END SUBROUTINE thomas_fermi_lsd_2
461 :
462 : ! **************************************************************************************************
463 : !> \brief ...
464 : !> \param rhoa ...
465 : !> \param r13a ...
466 : !> \param e_rho_rho_rho ...
467 : !> \param npoints ...
468 : ! **************************************************************************************************
469 0 : SUBROUTINE thomas_fermi_lsd_3(rhoa, r13a, e_rho_rho_rho, npoints)
470 :
471 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
472 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho
473 : INTEGER, INTENT(in) :: npoints
474 :
475 : INTEGER :: ip
476 : REAL(KIND=dp) :: f
477 :
478 0 : f = -f13*f23*f53*flsd
479 :
480 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
481 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho_rho,f,r13a)
482 : DO ip = 1, npoints
483 :
484 : IF (rhoa(ip) > eps_rho) THEN
485 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13a(ip)*rhoa(ip))
486 : END IF
487 :
488 : END DO
489 :
490 0 : END SUBROUTINE thomas_fermi_lsd_3
491 :
492 : END MODULE xc_thomas_fermi
493 :
|
