(git:08825b9)
Loading...
Searching...
No Matches
xc_semilocal_unittest.F
Go to the documentation of this file.
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
8PROGRAM xc_semilocal_unittest
9 USE ieee_arithmetic, ONLY: ieee_is_finite
10 USE input_cp2k_xc, ONLY: create_xc_fun_section
11 USE input_section_types, ONLY: section_release, &
12 section_type, &
13 section_vals_create, &
14 section_vals_release, &
15 section_vals_type, &
16 section_vals_val_set
17 USE kinds, ONLY: dp
18 USE libcp2k, ONLY: cp2k_finalize, &
19 cp2k_init
20 USE mathconstants, ONLY: pi
21 USE xc_atom, ONLY: xc_rho_set_atom_update
22 USE xc_derivative_desc, ONLY: deriv_norm_drho, &
23 deriv_norm_drhoa, &
24 deriv_norm_drhob, &
25 deriv_rho, &
26 deriv_rhoa, &
27 deriv_rhob, &
28 deriv_tau
29 USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
30 xc_dset_create, &
31 xc_dset_get_derivative, &
32 xc_dset_release
33 USE xc_derivative_types, ONLY: xc_derivative_get, &
34 xc_derivative_type
35 USE xc_derivatives, ONLY: xc_functionals_eval, &
36 xc_functionals_get_needs
37 USE xc_input_constants, ONLY: xc_pbe_orig, &
38 xc_pbe_rev, &
39 xc_pbe_sol
40 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
41 USE xc_rho_set_types, ONLY: xc_rho_set_create, &
42 xc_rho_set_release, &
43 xc_rho_set_type
44#if defined (__LIBXC)
45 USE xc_derivative_desc, ONLY: deriv_gamma, deriv_gamma_aa, deriv_gamma_ab, deriv_gamma_bb, &
46 deriv_tau_a, &
47 deriv_tau_b
48 USE xc_libxc, ONLY: libxc_spin_polarized_eval, &
49 libxc_spin_unpolarized_eval
50#endif
51
52 IMPLICIT NONE
53 INTEGER, PARAMETER :: labels(5) = [deriv_rhoa, deriv_rhob, deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob]
54 INTEGER, PARAMETER :: variants(3) = [xc_pbe_orig, xc_pbe_rev, xc_pbe_sol]
55 TYPE(section_type), POINTER :: specification
56 TYPE(section_vals_type), POINTER :: functional
57 INTEGER :: axis, configuration, spin, variant
58 REAL(dp) :: energy, gradient(5), hessian(5, 5), minus, plus, gm(5), gp(5), scratch(5, 5), &
59 step, trial(5), x(5), energy_error, hessian_error
60
61 CALL cp2k_init()
62 NULLIFY (specification, functional)
63 CALL create_xc_fun_section(specification)
64 CALL section_vals_create(functional, specification)
65 CALL section_vals_val_set(functional, 'PBE%_SECTION_PARAMETERS_', l_val=.true.)
66 energy_error = 0.0_dp
67 hessian_error = 0.0_dp
68 DO variant = 1, SIZE(variants)
69 CALL section_vals_val_set(functional, 'PBE%PARAMETRIZATION', i_val=variants(variant))
70 DO configuration = 1, 3
71 DO spin = 1, 2
72 IF (configuration == 1 .AND. spin == 2) cycle
73 x = [0.4_dp, 0.2_dp, 0.3_dp, 0.2_dp, 0.25_dp]
74 IF (configuration == 2) x(spin) = -0.1_dp
75 IF (configuration == 3) x(spin) = 1.e-13_dp
76 CALL evaluate(x, energy, gradient, hessian)
77 DO axis = 1, 5
78 step = 1.e-5_dp
79 IF (configuration == 3 .AND. axis == spin) step = 1.e-14_dp
80 trial = x
81 trial(axis) = x(axis) + step
82 CALL evaluate(trial, plus, gp, scratch)
83 trial(axis) = x(axis) - step
84 CALL evaluate(trial, minus, gm, scratch)
85 energy_error = max(energy_error, abs((plus - minus)/(2.0_dp*step) - gradient(axis)))
86 hessian_error = max(hessian_error, maxval(abs((gp - gm)/(2.0_dp*step) - hessian(:, axis))))
87 END DO
88 END DO
89 END DO
90 END DO
91 WRITE (*, '(A,2ES15.6)') 'PBE spin-floor first/second derivative errors:', energy_error, hessian_error
92 IF (energy_error > 2.e-8_dp) error stop 'PBE potential is not the derivative of its clamped-density energy'
93 IF (hessian_error > 2.e-7_dp) error stop 'PBE kernel is not the derivative of its clamped-density potential'
94 CALL section_vals_val_set(functional, 'PBE%_SECTION_PARAMETERS_', l_val=.false.)
95 CALL section_vals_val_set(functional, 'TPSS%_SECTION_PARAMETERS_', l_val=.true.)
96 CALL check_tpss_zero_gradient()
97 CALL section_vals_val_set(functional, 'TPSS%_SECTION_PARAMETERS_', l_val=.false.)
98 CALL section_vals_release(functional)
99 CALL section_vals_create(functional, specification)
100 CALL section_vals_val_set(functional, 'XWPBE%_SECTION_PARAMETERS_', l_val=.true.)
101 CALL check_xwpbe_limits()
102 CALL section_vals_val_set(functional, 'XWPBE%_SECTION_PARAMETERS_', l_val=.false.)
103#if defined (__LIBXC)
104 CALL check_libxc_tau_projection()
105 CALL check_libxc_c2_window()
106 CALL check_libxc_c2_rounding()
107 CALL check_libxc_individual_floor()
108#endif
109 CALL section_vals_release(functional)
110 CALL section_release(specification)
111 CALL cp2k_finalize()
112
113CONTAINS
114
115#if defined (__LIBXC)
116! **************************************************************************************************
117!> \brief Check meta-GGA potentials and kernels against independent differences of the projected energy.
118! **************************************************************************************************
119 SUBROUTINE check_libxc_tau_projection()
120 CHARACTER(len=32), PARAMETER :: names(6) = [character(len=32):: 'MGGA_X_TPSS', 'MGGA_C_TPSS',&
121 'MGGA_X_R2SCAN', 'MGGA_C_R2SCAN', 'LDA_C_PW', 'GGA_C_PBE']
122 REAL(dp), PARAMETER :: step = 1.e-5_dp
123
124 INTEGER :: axis, configuration, family, nconfig, &
125 nvar, spin_case
126 REAL(dp) :: energy, first_error, minus, plus, &
127 second_error
128 REAL(dp), ALLOCATABLE :: gm(:), gp(:), gradient(:), &
129 hessian(:, :), scratch(:, :), &
130 trial(:), x(:)
131
132 first_error = 0.0_dp
133 second_error = 0.0_dp
134 DO spin_case = 1, 2
135 nvar = 3
136 nconfig = 3
137 IF (spin_case == 2) THEN
138 nvar = 7
139 nconfig = 6
140 END IF
141 ALLOCATE (x(nvar), trial(nvar), gradient(nvar), gp(nvar), gm(nvar), &
142 hessian(nvar, nvar), scratch(nvar, nvar))
143 DO family = 1, SIZE(names)
144 DO configuration = 1, nconfig
145 IF (spin_case == 1) THEN
146 x(:) = [0.4_dp, 0.8_dp, 0.3_dp]
147 IF (configuration == 2) x(3) = 0.02_dp
148 IF (configuration == 3) x(3) = -0.02_dp
149 ELSE
150 x(:) = [0.3_dp, 0.2_dp, 0.6_dp, 0.4_dp, 0.3_dp, 0.12_dp, 0.10_dp]
151 IF (configuration == 2) x(6:7) = [0.02_dp, 0.015_dp]
152 IF (configuration == 3) x(6) = -0.02_dp
153 IF (configuration == 4) x(6:7) = [-0.02_dp, -0.015_dp]
154 IF (configuration == 5) x(1) = -0.1_dp
155 IF (configuration == 6) x(2) = -0.1_dp
156 END IF
157 CALL evaluate_libxc_point(names(family), x, energy, gradient, hessian)
158 IF (configuration >= 3) THEN
159 IF (index(names(family), 'MGGA_X_') == 1 .AND. abs(energy) < 1.e-4_dp) THEN
160 error stop 'Libxc discarded finite density at nonpositive tau'
161 END IF
162 trial(:) = x(:)
163 trial(nvar - spin_case + 1) = 1.e-10_dp
164 CALL evaluate_libxc_point(names(family), trial, plus, gp, scratch)
165 trial(nvar - spin_case + 1) = -1.e-10_dp
166 CALL evaluate_libxc_point(names(family), trial, minus, gm, scratch)
167 IF (abs(plus - minus) > 1.e-9_dp) error stop 'Libxc energy is discontinuous through zero tau'
168 END IF
169 DO axis = 1, nvar
170 trial(:) = x(:)
171 trial(axis) = x(axis) + step
172 CALL evaluate_libxc_point(names(family), trial, plus, gp, scratch)
173 trial(axis) = x(axis) - step
174 CALL evaluate_libxc_point(names(family), trial, minus, gm, scratch)
175 first_error = max(first_error, &
176 abs((plus - minus)/(2.0_dp*step) - gradient(axis))/(1.0_dp + abs(gradient(axis))))
177 second_error = max(second_error, &
178 maxval(abs((gp - gm)/(2.0_dp*step) - hessian(:, axis))/ &
179 (1.0_dp + abs(hessian(:, axis)))))
180 END DO
181 END DO
182 END DO
183 DEALLOCATE (x, trial, gradient, gp, gm, hessian, scratch)
184 END DO
185 WRITE (*, '(A,2ES15.6)') 'Libxc projected-tau first/second derivative errors:', first_error, second_error
186 IF (first_error > 3.e-8_dp) error stop 'Libxc projected-tau potential is not the energy derivative'
187 IF (second_error > 2.e-6_dp) error stop 'Libxc projected-tau kernel is not the potential derivative'
188 END SUBROUTINE check_libxc_tau_projection
189
190! **************************************************************************************************
191!> \brief ...
192! **************************************************************************************************
193 SUBROUTINE check_libxc_c2_window()
194 CHARACTER(len=32), PARAMETER :: names(4) = [character(len=32):: 'MGGA_X_TPSS', 'MGGA_C_TPSS',&
195 'MGGA_X_R2SCAN', 'MGGA_C_R2SCAN']
196 REAL(dp), PARAMETER :: offsets(3) = [-0.6_dp, 0.2_dp, 0.6_dp], &
197 relative_width = 64.0_dp*sqrt(epsilon(0.0_dp))
198
199 INTEGER :: axis, configuration, family, nvar, &
200 spin_case
201 REAL(dp) :: em, emm, energy, ep, epp, first_error, &
202 first_fd, h, second_error, &
203 third_error, wa, wb, widtha, widthb
204 REAL(dp), ALLOCATABLE :: gm(:), gmm(:), gp(:), gpp(:), gradient(:), hessian(:, :), hm(:, :), &
205 hmm(:, :), hp(:, :), hpp(:, :), third(:, :, :), trial(:), x(:)
206
207! Nonzero offsets keep the third-derivative differences well conditioned.
208
209 first_error = 0.0_dp
210 second_error = 0.0_dp
211 third_error = 0.0_dp
212 DO spin_case = 1, 2
213 nvar = 3
214 IF (spin_case == 2) nvar = 7
215 ALLOCATE (x(nvar), trial(nvar), gradient(nvar), gm(nvar), gp(nvar), &
216 gmm(nvar), gpp(nvar), hessian(nvar, nvar), &
217 third(nvar, nvar, nvar), hm(nvar, nvar), hp(nvar, nvar), hmm(nvar, nvar), hpp(nvar, nvar))
218 DO family = 1, SIZE(names)
219 DO configuration = 1, SIZE(offsets)
220 IF (spin_case == 1) THEN
221 x(:) = [0.4_dp, 0.64_dp, 0.3_dp]
222 wa = x(2)/(8.0_dp*x(1))
223 widtha = relative_width*wa
224 x(3) = wa + offsets(configuration)*widtha
225 ELSE
226 x(:) = [0.3_dp, 0.2_dp, 0.055_dp, 0.16_dp, 0.09_dp, 0.12_dp, 0.10_dp]
227 wa = x(4)/(8.0_dp*x(1))
228 wb = x(5)/(8.0_dp*x(2))
229 widtha = relative_width*wa
230 widthb = relative_width*wb
231 x(6) = wa + offsets(configuration)*widtha
232 x(7) = wb - offsets(configuration)*widthb
233 END IF
234 ! The installed r2SCAN functions expose only first and second derivatives.
235 IF (family <= 2) THEN
236 CALL evaluate_libxc_point(names(family), x, energy, gradient, hessian, third_all=third, contractions=.true.)
237 ELSE
238 CALL evaluate_libxc_point(names(family), x, energy, gradient, hessian, contractions=.true.)
239 END IF
240 DO axis = 1, nvar
241 h = 0.01_dp*widtha
242 IF (spin_case == 1) THEN
243 IF (axis == 1) h = 0.01_dp*widtha/(wa/x(1))
244 IF (axis == 2) h = 0.01_dp*widtha/(wa/x(2))
245 ELSE
246 IF (axis == 1) h = 0.01_dp*widtha/(wa/x(1))
247 IF (axis == 2) h = 0.01_dp*widthb/(wb/x(2))
248 IF (axis == 3) h = 1.e-5_dp
249 IF (axis == 4) h = 0.01_dp*widtha/(wa/x(4))
250 IF (axis == 5) h = 0.01_dp*widthb/(wb/x(5))
251 IF (axis == 7) h = 0.01_dp*widthb
252 END IF
253 trial(:) = x(:); trial(axis) = x(axis) + h
254 CALL evaluate_libxc_point(names(family), trial, ep, gp, hp, contractions=.true.)
255 trial(axis) = x(axis) - h
256 CALL evaluate_libxc_point(names(family), trial, em, gm, hm, contractions=.true.)
257 trial(axis) = x(axis) + 2.0_dp*h
258 CALL evaluate_libxc_point(names(family), trial, epp, gpp, hpp, contractions=.true.)
259 trial(axis) = x(axis) - 2.0_dp*h
260 CALL evaluate_libxc_point(names(family), trial, emm, gmm, hmm, contractions=.true.)
261 first_fd = (8.0_dp*(ep - em) - (epp - emm))/(12.0_dp*h)
262 first_error = max(first_error, abs(first_fd - gradient(axis))/(1.0_dp + abs(gradient(axis))))
263 second_error = max(second_error, &
264 maxval(abs((8.0_dp*(gp - gm) - (gpp - gmm))/(12.0_dp*h) - hessian(:, axis))/(1.0_dp + abs(hessian(:, axis)))))
265 IF (family <= 2) THEN
266 third_error = max(third_error, &
267 maxval(abs((8.0_dp*(hp - hm) - (hpp - hmm))/(12.0_dp*h) - third(:, :, axis))/ &
268 (1.0_dp + abs(third(:, :, axis)))))
269 END IF
270 END DO
271 END DO
272 END DO
273 DEALLOCATE (x, trial, gradient, gm, gp, gmm, gpp, hessian, third, hm, hp, hmm, hpp)
274 END DO
275 WRITE (*, '(A,3ES15.6)') 'Libxc C2-window first/second/third derivative errors:', first_error, second_error, third_error
276 IF (first_error > 2.e-6_dp) error stop 'C2-window energy Jacobian failed'
277 IF (second_error > 2.e-5_dp) error stop 'C2-window potential Hessian failed'
278 IF (third_error > 2.e-5_dp) error stop 'C2-window kernel third derivatives failed'
279 END SUBROUTINE check_libxc_c2_window
280
281! **************************************************************************************************
282!> \brief ...
283! **************************************************************************************************
284 SUBROUTINE check_libxc_c2_rounding()
285 CHARACTER(len=32), PARAMETER :: names(4) = [character(len=32):: 'MGGA_X_TPSS', 'MGGA_C_TPSS',&
286 'MGGA_X_R2SCAN', 'MGGA_C_R2SCAN']
287 REAL(dp), PARAMETER :: relative_width = 64.0_dp*sqrt(epsilon(0.0_dp))
288
289 INTEGER :: channel, family, nvar, spin_case, &
290 tau_axis
291 REAL(dp) :: curve, energy, error, expected, slope, &
292 u, w(2), width(2)
293 REAL(dp), ALLOCATABLE :: gp(:), gradient(:), hessian(:, :), &
294 hp(:, :), probe(:), x(:)
295
296 error = 0.0_dp
297 DO spin_case = 1, 2
298 nvar = 3
299 IF (spin_case == 2) nvar = 7
300 ALLOCATE (x(nvar), probe(nvar), gradient(nvar), gp(nvar), hessian(nvar, nvar), hp(nvar, nvar))
301 DO family = 1, SIZE(names)
302 IF (spin_case == 1) THEN
303 x(:) = [0.4_dp, 0.8_dp, 0.3_dp]
304 w(1) = x(2)**2/(8.0_dp*x(1))
305 ELSE
306 x(:) = [0.3_dp, 0.2_dp, 0.6_dp, 0.4_dp, 0.3_dp, 0.12_dp, 0.10_dp]
307 w(1) = x(4)**2/(8.0_dp*x(1))
308 w(2) = x(5)**2/(8.0_dp*x(2))
309 END IF
310 probe(:) = x(:)
311 DO channel = 1, spin_case
312 tau_axis = nvar - spin_case + channel
313 width(channel) = relative_width*w(channel)
314 x(tau_axis) = w(channel) + 0.99999_dp*width(channel)
315 ! The nearby unprojected point supplies raw Libxc derivatives.
316 probe(tau_axis) = w(channel) + 1.00001_dp*width(channel)
317 END DO
318 CALL evaluate_libxc_point(names(family), x, energy, gradient, hessian)
319 CALL evaluate_libxc_point(names(family), probe, energy, gp, hp)
320 DO channel = 1, spin_case
321 tau_axis = nvar - spin_case + channel
322 u = (x(tau_axis) - w(channel))/width(channel)
323 slope = (1.0_dp + u)**2*(2.0_dp - u)/4.0_dp
324 curve = 0.75_dp*(1.0_dp - u*u)/width(channel)
325 expected = hp(tau_axis, tau_axis)*slope*slope + gp(tau_axis)*curve
326 error = max(error, abs(hessian(tau_axis, tau_axis) - expected)/(1.0_dp + abs(expected)))
327 END DO
328 END DO
329 DEALLOCATE (x, probe, gradient, gp, hessian, hp)
330 END DO
331 WRITE (*, '(A,ES15.6)') 'Libxc rounded-window curvature error:', error
332 IF (error > 2.e-6_dp) error stop 'Rounded energy value discarded C2 projection curvature'
333 END SUBROUTINE check_libxc_c2_rounding
334
335! **************************************************************************************************
336!> \brief Check the original-density energy normalization below Libxc's individual spin floors.
337! **************************************************************************************************
338 SUBROUTINE check_libxc_individual_floor()
339 CHARACTER(len=24), PARAMETER :: names(6) = [character(len=24) :: 'LDA_C_PW', 'GGA_C_PBE', &
340 'MGGA_X_TPSS', 'MGGA_C_TPSS', 'MGGA_X_R2SCAN', 'MGGA_C_R2SCAN']
341
342 INTEGER :: axis, configuration, f
343 REAL(dp) :: energy, error_first, error_second, error_third, gm(7), gp(7), gradient(7), &
344 hessian(7, 7), hm(7, 7), hp(7, 7), minus, plus, step, third(2, 2, 2), trial(7), x(7)
345
346 error_first = 0.0_dp
347 error_second = 0.0_dp
348 error_third = 0.0_dp
349 DO f = 1, SIZE(names)
350 DO configuration = 1, 6
351 x = [0.6_dp, 0.01_dp, 0.21_dp, 0.2_dp, 0.02_dp, 0.3_dp, 0.2_dp]
352 SELECT CASE (configuration)
353 CASE (2)
354 x(2) = -0.01_dp
355 CASE (3)
356 x(1:2) = 0.03_dp
357 CASE (4)
358 x(1:2) = [0.01_dp, 0.6_dp]
359 CASE (5)
360 x(7) = -0.02_dp ! Compose the density floor with the kinetic-density projection.
361 CASE (6)
362 x(1:2) = [-0.01_dp, 0.6_dp]
363 END SELECT
364 IF (f <= 2) THEN
365 CALL evaluate_libxc_point(trim(names(f)), x, energy, gradient, hessian, 0.05_dp, third)
366 ELSE
367 CALL evaluate_libxc_point(trim(names(f)), x, energy, gradient, hessian, 0.05_dp)
368 END IF
369 DO axis = 1, SIZE(x)
370 step = 1.e-6_dp
371 trial = x
372 trial(axis) = x(axis) + step
373 CALL evaluate_libxc_point(trim(names(f)), trial, plus, gp, hp, 0.05_dp)
374 trial(axis) = x(axis) - step
375 CALL evaluate_libxc_point(trim(names(f)), trial, minus, gm, hm, 0.05_dp)
376 error_first = max(error_first, abs((plus - minus)/(2.0_dp*step) - gradient(axis)))
377 error_second = max(error_second, maxval(abs((gp - gm)/(2.0_dp*step) - hessian(:, axis))))
378 IF (f <= 2 .AND. axis <= 2) THEN
379 error_third = max(error_third, maxval(abs((hp(1:2, 1:2) - hm(1:2, 1:2))/ &
380 (2.0_dp*step) - third(:, :, axis))))
381 END IF
382 END DO
383 END DO
384 END DO
385 WRITE (*, '(A,3ES15.6)') 'Libxc individual spin-floor first/second/third derivative errors:', &
386 error_first, error_second, error_third
387 IF (error_first > 2.e-7_dp) error stop 'Libxc individual spin-floor potential is not an energy derivative'
388 IF (error_second > 2.e-6_dp) error stop 'Libxc individual spin-floor kernel is not a potential derivative'
389 IF (error_third > 2.e-5_dp) error stop 'Libxc individual spin-floor third derivatives are inconsistent'
390 END SUBROUTINE check_libxc_individual_floor
391
392! **************************************************************************************************
393!> \brief Evaluate Libxc with independently varied density, gradient-norm and kinetic-density inputs.
394!> \param name Functional name.
395!> \param x Three restricted variables or seven spin variables.
396!> \param energy Energy density.
397!> \param gradient First derivatives with a preceding functional contribution removed.
398!> \param hessian Second derivatives in the same independent variables.
399!> \param density_cutoff Optional cutoff for the individual-spin-floor oracle.
400!> \param third_rho Optional third derivatives in the two spin densities.
401!> \param third_all Optional third derivatives in every independent variable.
402!> \param contractions Use Libxc gradient contractions instead of CP2K gradient norms.
403! **************************************************************************************************
404 SUBROUTINE evaluate_libxc_point(name, x, energy, gradient, hessian, density_cutoff, third_rho, third_all, contractions)
405 CHARACTER(len=*), INTENT(IN) :: name
406 REAL(dp), INTENT(IN) :: x(:)
407 REAL(dp), INTENT(OUT) :: energy, gradient(size(x)), &
408 hessian(size(x), size(x))
409 REAL(dp), INTENT(IN), OPTIONAL :: density_cutoff
410 REAL(dp), INTENT(OUT), OPTIONAL :: third_rho(2, 2, 2), third_all(size(x), size(x), size(x))
411 LOGICAL, INTENT(IN), OPTIONAL :: contractions
412
413 INTEGER :: bounds(2, 3), i, j, k, labels(7), &
414 nspins, order
415 LOGICAL :: use_contractions
416 REAL(dp) :: cutoff
417 REAL(dp), POINTER :: values(:, :, :)
418 TYPE(section_vals_type), POINTER :: parameters
419 TYPE(xc_derivative_set_type) :: derivatives
420 TYPE(xc_derivative_type), POINTER :: derivative
421 TYPE(xc_rho_cflags_type) :: needs
422 TYPE(xc_rho_set_type) :: rho_set
423
424 bounds = 1
425 use_contractions = .false.
426 IF (PRESENT(contractions)) use_contractions = contractions
427 labels = 0
428 needs%norm_drho = .true.
429 IF (SIZE(x) == 3) THEN
430 nspins = 1
431 labels(1:3) = [deriv_rho, deriv_norm_drho, deriv_tau]
432 needs%rho = .true.
433 needs%tau = .true.
434 needs%laplace_rho = .true.
435 ELSE
436 nspins = 2
437 labels = [deriv_rhoa, deriv_rhob, deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob, deriv_tau_a, deriv_tau_b]
438 needs%rho_spin = .true.
439 needs%norm_drho_spin = .true.
440 needs%tau_spin = .true.
441 needs%laplace_rho_spin = .true.
442 END IF
443 IF (use_contractions) THEN
444 IF (nspins == 1) THEN
445 labels(2) = deriv_gamma
446 ELSE
447 labels(3:5) = [deriv_gamma_ab, deriv_gamma_aa, deriv_gamma_bb]
448 END IF
449 END IF
450 cutoff = 1.e-20_dp
451 IF (PRESENT(density_cutoff)) cutoff = density_cutoff
452 CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=cutoff, tau_cutoff=1.e-20_dp)
453 CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
454 IF (nspins == 1) THEN
455 rho_set%rho = x(1)
456 rho_set%norm_drho = x(2)
457 rho_set%tau = x(3)
458 rho_set%laplace_rho = 0.0_dp
459 IF (use_contractions) rho_set%norm_drho = sqrt(x(2))
460 ELSE
461 rho_set%rhoa = x(1)
462 rho_set%rhob = x(2)
463 rho_set%norm_drho = x(3)
464 rho_set%norm_drhoa = x(4)
465 rho_set%norm_drhob = x(5)
466 rho_set%tau_a = x(6)
467 rho_set%tau_b = x(7)
468 rho_set%laplace_rhoa = 0.0_dp
469 rho_set%laplace_rhob = 0.0_dp
470 IF (use_contractions) THEN
471 rho_set%norm_drho = sqrt(2.0_dp*x(3) + x(4) + x(5))
472 rho_set%norm_drhoa = sqrt(x(4))
473 rho_set%norm_drhob = sqrt(x(5))
474 END IF
475 END IF
476 rho_set%has = needs
477 CALL xc_dset_create(derivatives, local_bounds=bounds)
478 DO i = 1, SIZE(x)
479 derivative => xc_dset_get_derivative(derivatives, [labels(i)], allocate_deriv=.true.)
480 CALL xc_derivative_get(derivative, deriv_data=values)
481 values = 0.125_dp
482 DO j = i, SIZE(x)
483 derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j)], allocate_deriv=.true.)
484 CALL xc_derivative_get(derivative, deriv_data=values)
485 values = 0.125_dp
486 END DO
487 END DO
488 order = 2
489 IF (PRESENT(third_rho) .OR. PRESENT(third_all)) THEN
490 IF (PRESENT(third_rho) .AND. nspins /= 2) error stop 'Spin third-derivative oracle requires two spins'
491 order = 3
492 DO i = 1, SIZE(x)
493 DO j = 1, SIZE(x)
494 DO k = 1, SIZE(x)
495 derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j), labels(k)], allocate_deriv=.true.)
496 CALL xc_derivative_get(derivative, deriv_data=values)
497 values = 0.125_dp
498 END DO
499 END DO
500 END DO
501 END IF
502 NULLIFY (parameters)
503 IF (nspins == 1) THEN
504 CALL libxc_spin_unpolarized_eval(rho_set, derivatives, order, parameters, func_name_override=name)
505 ELSE
506 CALL libxc_spin_polarized_eval(rho_set, derivatives, order, parameters, func_name_override=name)
507 END IF
508 derivative => xc_dset_get_derivative(derivatives, [INTEGER::])
509 CALL xc_derivative_get(derivative, deriv_data=values)
510 energy = values(1, 1, 1)
511 DO i = 1, SIZE(x)
512 derivative => xc_dset_get_derivative(derivatives, [labels(i)])
513 CALL xc_derivative_get(derivative, deriv_data=values)
514 gradient(i) = values(1, 1, 1) - 0.125_dp
515 DO j = i, SIZE(x)
516 derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j)])
517 CALL xc_derivative_get(derivative, deriv_data=values)
518 hessian(i, j) = values(1, 1, 1) - 0.125_dp
519 hessian(j, i) = hessian(i, j)
520 END DO
521 END DO
522 IF (PRESENT(third_rho)) THEN
523 DO i = 1, 2
524 DO j = 1, 2
525 DO k = 1, 2
526 derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j), labels(k)])
527 CALL xc_derivative_get(derivative, deriv_data=values)
528 third_rho(i, j, k) = values(1, 1, 1) - 0.125_dp
529 END DO
530 END DO
531 END DO
532 END IF
533 IF (PRESENT(third_all)) THEN
534 DO k = 1, SIZE(x)
535 DO j = 1, SIZE(x)
536 DO i = 1, SIZE(x)
537 derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j), labels(k)])
538 CALL xc_derivative_get(derivative, deriv_data=values)
539 third_all(i, j, k) = values(1, 1, 1) - 0.125_dp
540 END DO
541 END DO
542 END DO
543 IF (.NOT. all(ieee_is_finite(third_all))) error stop 'Nonfinite Libxc third derivatives'
544 END IF
545 IF (.NOT. ieee_is_finite(energy) .OR. .NOT. all(ieee_is_finite(gradient)) .OR. &
546 .NOT. all(ieee_is_finite(hessian))) error stop 'Nonfinite Libxc projected-tau result'
547 CALL xc_dset_release(derivatives)
548 CALL xc_rho_set_release(rho_set)
549 END SUBROUTINE evaluate_libxc_point
550#endif
551
552! **************************************************************************************************
553!> \brief Evaluate native PBE in its independent spin-density and gradient-norm variables.
554!> \param x Two spin densities followed by total, alpha and beta gradient norms.
555!> \param energy Energy density.
556!> \param gradient First derivatives, with the accumulation of another functional preserved.
557!> \param hessian Second derivatives in the same independent variables.
558!> \param derivative_order ...
559! **************************************************************************************************
560 SUBROUTINE evaluate(x, energy, gradient, hessian, derivative_order)
561 REAL(dp), INTENT(IN) :: x(5)
562 REAL(dp), INTENT(OUT) :: energy, gradient(5), hessian(5, 5)
563 INTEGER, INTENT(IN), OPTIONAL :: derivative_order
564
565 INTEGER :: bounds(2, 3), i, j, order
566 REAL(dp), POINTER :: values(:, :, :)
567 TYPE(xc_derivative_set_type) :: derivatives
568 TYPE(xc_derivative_type), POINTER :: derivative
569 TYPE(xc_rho_cflags_type) :: needs
570 TYPE(xc_rho_set_type) :: rho_set
571
572 order = 2
573 IF (PRESENT(derivative_order)) order = derivative_order
574 bounds = 1
575 needs = xc_functionals_get_needs(functional, lsd=.true., calc_potential=.true.)
576 needs%norm_drho = .true.
577 CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=1.e-20_dp)
578 CALL xc_rho_set_atom_update(rho_set, needs, 2, bounds)
579 rho_set%rhoa = x(1)
580 rho_set%rhob = x(2)
581 rho_set%norm_drho = x(3)
582 rho_set%norm_drhoa = x(4)
583 rho_set%norm_drhob = x(5)
584 rho_set%has%rho_spin = .true.
585 rho_set%has%norm_drho = .true.
586 rho_set%has%norm_drho_spin = .true.
587 CALL xc_dset_create(derivatives, local_bounds=bounds)
588 ! Seed every first/second derivative to catch overwriting earlier functional contributions.
589 DO i = 1, 5
590 derivative => xc_dset_get_derivative(derivatives, [labels(i)], allocate_deriv=.true.)
591 CALL xc_derivative_get(derivative, deriv_data=values)
592 values = 0.125_dp
593 DO j = i, 5
594 derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j)], allocate_deriv=.true.)
595 CALL xc_derivative_get(derivative, deriv_data=values)
596 values = 0.125_dp
597 END DO
598 END DO
599 CALL xc_functionals_eval(functional, lsd=.true., rho_set=rho_set, deriv_set=derivatives, deriv_order=order)
600 derivative => xc_dset_get_derivative(derivatives, [INTEGER::])
601 CALL xc_derivative_get(derivative, deriv_data=values)
602 energy = values(1, 1, 1)
603 DO i = 1, 5
604 derivative => xc_dset_get_derivative(derivatives, [labels(i)])
605 CALL xc_derivative_get(derivative, deriv_data=values)
606 gradient(i) = values(1, 1, 1) - 0.125_dp
607 DO j = i, 5
608 derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j)])
609 CALL xc_derivative_get(derivative, deriv_data=values)
610 hessian(i, j) = values(1, 1, 1) - 0.125_dp
611 IF (order == 1) hessian(i, j) = 0.0_dp
612 hessian(j, i) = hessian(i, j)
613 END DO
614 END DO
615 CALL xc_dset_release(derivatives)
616 CALL xc_rho_set_release(rho_set)
617 END SUBROUTINE evaluate
618
619! **************************************************************************************************
620!> \brief Check native TPSS derivatives at stationary density points and below the von Weizsaecker bound.
621! **************************************************************************************************
622 SUBROUTINE check_tpss_zero_gradient()
623 REAL(dp), PARAMETER :: step = 1.e-5_dp, tau_scales(3) = [0.2_dp, 1.0_dp, 2.0_dp]
624
625 INTEGER :: axis, configuration
626 REAL(dp) :: energy, error, gradient(3), minus, plus, &
627 scratch(3), tau_unif, trial(3), x(3)
628
629 error = 0.0_dp
630 x(1) = 0.4_dp
631 tau_unif = 0.3_dp*(3.0_dp*pi**2)**(2.0_dp/3.0_dp)*x(1)**(5.0_dp/3.0_dp)
632 DO configuration = 1, SIZE(tau_scales)
633 x(2) = 0.0_dp
634 x(3) = tau_scales(configuration)*tau_unif
635 CALL evaluate_tpss_point(x, energy, gradient)
636 IF (abs(gradient(2)) > 1.e-13_dp) error stop 'TPSS has a nonzero derivative at zero density gradient'
637 DO axis = 1, 3, 2
638 trial(:) = x(:)
639 trial(axis) = x(axis) + step
640 CALL evaluate_tpss_point(trial, plus, scratch)
641 trial(axis) = x(axis) - step
642 CALL evaluate_tpss_point(trial, minus, scratch)
643 error = max(error, abs((plus - minus)/(2.0_dp*step) - gradient(axis)))
644 END DO
645 trial = x
646 trial(2) = 1.e-100_dp
647 CALL evaluate_tpss_point(trial, plus, scratch)
648 IF (abs(plus - energy) > 1.e-13_dp .OR. abs(scratch(2)) > 1.e-20_dp) THEN
649 error stop 'TPSS does not approach its zero-gradient limit continuously'
650 END IF
651 x(2) = 1.e-4_dp
652 CALL evaluate_tpss_point(x, energy, gradient)
653 trial = x
654 trial(2) = x(2) + step
655 CALL evaluate_tpss_point(trial, plus, scratch)
656 trial(2) = x(2) - step
657 CALL evaluate_tpss_point(trial, minus, scratch)
658 error = max(error, abs((plus - minus)/(2.0_dp*step) - gradient(2)))
659 END DO
660 ! NLCC can give tau below |grad rho|^2/(8*rho). Check the projection Jacobian.
661 x = [0.4_dp, 1.2_dp, 0.05_dp]
662 CALL evaluate_tpss_point(x, energy, gradient)
663 DO axis = 1, 3
664 trial = x
665 trial(axis) = x(axis) + step
666 CALL evaluate_tpss_point(trial, plus, scratch)
667 trial(axis) = x(axis) - step
668 CALL evaluate_tpss_point(trial, minus, scratch)
669 error = max(error, abs((plus - minus)/(2.0_dp*step) - gradient(axis)))
670 END DO
671 ! Negative Fourier-interpolated tau must not discard a finite-density energy.
672 DO configuration = 1, 2
673 x = [0.4_dp, 1.2_dp, -0.05_dp]
674 IF (configuration == 2) x(2) = 0.0_dp
675 CALL evaluate_tpss_point(x, energy, gradient)
676 IF (abs(energy) < 1.e-4_dp .OR. gradient(3) /= 0.0_dp) THEN
677 error stop 'TPSS negative-tau continuation is not finite and constant in tau'
678 END IF
679 DO axis = 1, 3
680 IF (axis == 2 .AND. configuration == 2) cycle
681 trial = x
682 trial(axis) = x(axis) + step
683 CALL evaluate_tpss_point(trial, plus, scratch)
684 trial(axis) = x(axis) - step
685 CALL evaluate_tpss_point(trial, minus, scratch)
686 error = max(error, abs((plus - minus)/(2.0_dp*step) - gradient(axis)))
687 END DO
688 END DO
689 x = [0.4_dp, 1.2_dp, 0.0_dp]
690 CALL evaluate_tpss_point(x, energy, gradient)
691 trial = x
692 trial(3) = 1.e-10_dp
693 CALL evaluate_tpss_point(trial, plus, scratch)
694 trial(3) = -1.e-10_dp
695 CALL evaluate_tpss_point(trial, minus, scratch)
696 IF (max(abs(plus - energy), abs(minus - energy)) > 1.e-9_dp) THEN
697 error stop 'TPSS energy is discontinuous through zero kinetic density'
698 END IF
699 WRITE (*, '(A,ES15.6)') 'TPSS limiting-density first-derivative error:', error
700 IF (error > 3.e-8_dp) error stop 'TPSS limiting-density potential is not the energy derivative'
701 END SUBROUTINE check_tpss_zero_gradient
702
703! **************************************************************************************************
704!> \brief Check screened exchange limits and gradient scaling against energy differences.
705! **************************************************************************************************
706 SUBROUTINE check_xwpbe_limits()
707 REAL(dp), PARAMETER :: reduced_gradients(3) = [0.0_dp, 0.2_dp, 12.0_dp], step = 1.e-5_dp
708
709 INTEGER :: axis, configuration, range_case
710 REAL(dp) :: energy, error, gr(3), gs(5), &
711 hessian(5, 5), minus, plus, &
712 scratchr(3), scratchs(5), trialr(3), &
713 trials(5), xr(3), xs(5)
714
715 error = 0.0_dp
716 CALL section_vals_val_set(functional, 'XWPBE%SCALE_X', r_val=-0.25_dp)
717 CALL section_vals_val_set(functional, 'XWPBE%SCALE_X0', r_val=1.0_dp)
718 DO range_case = 1, 2
719 CALL section_vals_val_set(functional, 'XWPBE%OMEGA', r_val=0.11_dp*real(range_case - 1, dp))
720 DO configuration = 1, SIZE(reduced_gradients)
721 xr = [0.4_dp, 0.0_dp, 0.5_dp]
722 xr(2) = reduced_gradients(configuration)*2.0_dp*(3.0_dp*pi**2)**(1.0_dp/3.0_dp)*xr(1)**(4.0_dp/3.0_dp)
723 CALL evaluate_tpss_point(xr, energy, gr)
724 DO axis = 1, 2
725 IF (axis == 2 .AND. configuration == 1) cycle
726 trialr = xr; trialr(axis) = xr(axis) + step
727 CALL evaluate_tpss_point(trialr, plus, scratchr)
728 trialr(axis) = xr(axis) - step
729 CALL evaluate_tpss_point(trialr, minus, scratchr)
730 error = max(error, abs((plus - minus)/(2.0_dp*step) - gr(axis))/(1.0_dp + abs(gr(axis))))
731 END DO
732 IF (configuration == 1) THEN
733 IF (abs(gr(2)) > 1.e-13_dp) error stop 'XWPBE zero-gradient potential is not stationary'
734 trialr = xr; trialr(2) = 1.e-8_dp
735 CALL evaluate_tpss_point(trialr, plus, scratchr)
736 IF (abs(plus - energy) > 1.e-12_dp) error stop 'XWPBE zero-gradient energy is discontinuous'
737 END IF
738 xs = [0.3_dp, 0.2_dp, 0.3_dp, 0.0_dp, 0.0_dp]
739 xs(4) = reduced_gradients(configuration)*2.0_dp*(3.0_dp*pi**2)**(1.0_dp/3.0_dp)*(2.0_dp*xs(1))**(4.0_dp/3.0_dp)/2.0_dp
740 xs(5) = reduced_gradients(configuration)*2.0_dp*(3.0_dp*pi**2)**(1.0_dp/3.0_dp)*(2.0_dp*xs(2))**(4.0_dp/3.0_dp)/2.0_dp
741 CALL evaluate(xs, energy, gs, hessian, derivative_order=1)
742 DO axis = 1, 5
743 IF (axis >= 4 .AND. configuration == 1) cycle
744 trials = xs; trials(axis) = xs(axis) + step
745 CALL evaluate(trials, plus, scratchs, hessian, derivative_order=1)
746 trials(axis) = xs(axis) - step
747 CALL evaluate(trials, minus, scratchs, hessian, derivative_order=1)
748 error = max(error, abs((plus - minus)/(2.0_dp*step) - gs(axis))/(1.0_dp + abs(gs(axis))))
749 END DO
750 IF (configuration == 1 .AND. maxval(abs(gs(4:5))) > 1.e-13_dp) error stop 'XWPBE spin zero-gradient potential is not stationary'
751 END DO
752 END DO
753 WRITE (*, '(A,ES15.6)') 'XWPBE limiting/large-gradient first-derivative error:', error
754 IF (error > 3.e-8_dp) error stop 'XWPBE potential is not the derivative of its scaled-gradient energy'
755 END SUBROUTINE check_xwpbe_limits
756
757! **************************************************************************************************
758!> \brief Evaluate an unpolarized semilocal functional in density, gradient-norm and kinetic-density variables.
759!> \param x ...
760!> \param energy ...
761!> \param gradient ...
762! **************************************************************************************************
763 SUBROUTINE evaluate_tpss_point(x, energy, gradient)
764 REAL(dp), INTENT(IN) :: x(3)
765 REAL(dp), INTENT(OUT) :: energy, gradient(3)
766
767 INTEGER, PARAMETER :: labels(3) = [deriv_rho, deriv_norm_drho, deriv_tau]
768
769 INTEGER :: bounds(2, 3), i
770 REAL(dp), POINTER :: values(:, :, :)
771 TYPE(xc_derivative_set_type) :: derivatives
772 TYPE(xc_derivative_type), POINTER :: derivative
773 TYPE(xc_rho_cflags_type) :: needs
774 TYPE(xc_rho_set_type) :: rho_set
775
776 bounds = 1
777 needs = xc_functionals_get_needs(functional, lsd=.false., calc_potential=.true.)
778 needs%tau = .true.
779 CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=1.e-20_dp, tau_cutoff=1.e-20_dp)
780 CALL xc_rho_set_atom_update(rho_set, needs, 1, bounds)
781 rho_set%rho = x(1)
782 rho_set%norm_drho = x(2)
783 rho_set%tau = x(3)
784 rho_set%has%rho = .true.
785 rho_set%has%norm_drho = .true.
786 rho_set%has%tau = .true.
787 CALL xc_dset_create(derivatives, local_bounds=bounds)
788 DO i = 1, SIZE(labels)
789 derivative => xc_dset_get_derivative(derivatives, [labels(i)], allocate_deriv=.true.)
790 CALL xc_derivative_get(derivative, deriv_data=values)
791 values = 0.125_dp
792 END DO
793 CALL xc_functionals_eval(functional, lsd=.false., rho_set=rho_set, deriv_set=derivatives, deriv_order=1)
794 derivative => xc_dset_get_derivative(derivatives, [INTEGER::])
795 CALL xc_derivative_get(derivative, deriv_data=values)
796 energy = values(1, 1, 1)
797 DO i = 1, SIZE(labels)
798 derivative => xc_dset_get_derivative(derivatives, [labels(i)])
799 CALL xc_derivative_get(derivative, deriv_data=values)
800 gradient(i) = values(1, 1, 1) - 0.125_dp
801 END DO
802 IF (.NOT. ieee_is_finite(energy) .OR. .NOT. all(ieee_is_finite(gradient))) THEN
803 error stop 'TPSS energy or first derivative is nonfinite at a stationary density point'
804 END IF
805 CALL xc_dset_release(derivatives)
806 CALL xc_rho_set_release(rho_set)
807 END SUBROUTINE evaluate_tpss_point
808
809END PROGRAM xc_semilocal_unittest
input_cp2k_xc
function that build the xc section of the input
Definition input_cp2k_xc.F:14
input_cp2k_xc::create_xc_fun_section
subroutine, public create_xc_fun_section(section)
creates the structure of the section needed to select the xc functional
Definition input_cp2k_xc.F:77
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_create
recursive subroutine, public section_vals_create(section_vals, section)
creates a object where to store the values of a section
Definition input_section_types.F:612
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition input_section_types.F:1233
input_section_types::section_release
recursive subroutine, public section_release(section)
releases the given keyword list (see doc/ReferenceCounting.html)
Definition input_section_types.F:238
input_section_types::section_vals_release
recursive subroutine, public section_vals_release(section_vals)
releases the given object
Definition input_section_types.F:657
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
libcp2k
CP2K C/C++ interface.
Definition libcp2k.F:20
libcp2k::cp2k_finalize
subroutine, public cp2k_finalize()
...
Definition libcp2k.F:139
libcp2k::cp2k_init
subroutine, public cp2k_init()
...
Definition libcp2k.F:104
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
xc_atom
Definition xc_atom.F:9
xc_atom::xc_rho_set_atom_update
subroutine, public xc_rho_set_atom_update(rho_set, needs, nspins, bo)
...
Definition xc_atom.F:523
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition xc_derivative_desc.F:26
xc_derivative_desc::deriv_norm_drho
integer, parameter, public deriv_norm_drho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_gamma
integer, parameter, public deriv_gamma
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_gamma_aa
integer, parameter, public deriv_gamma_aa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_tau
integer, parameter, public deriv_tau
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_tau_b
integer, parameter, public deriv_tau_b
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_tau_a
integer, parameter, public deriv_tau_a
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_gamma_bb
integer, parameter, public deriv_gamma_bb
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_gamma_ab
integer, parameter, public deriv_gamma_ab
Definition xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition xc_derivative_set_types.F:65
xc_derivative_set_types::xc_dset_release
subroutine, public xc_dset_release(derivative_set)
releases a derivative set
Definition xc_derivative_set_types.F:143
xc_derivative_set_types::xc_dset_create
subroutine, public xc_dset_create(derivative_set, pw_pool, local_bounds)
creates a derivative set object
Definition xc_derivative_set_types.F:111
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition xc_derivative_types.F:103
xc_derivatives
Definition xc_derivatives.F:9
xc_derivatives::xc_functionals_get_needs
type(xc_rho_cflags_type) function, public xc_functionals_get_needs(functionals, lsd, calc_potential)
...
Definition xc_derivatives.F:774
xc_derivatives::xc_functionals_eval
subroutine, public xc_functionals_eval(functionals, lsd, rho_set, deriv_set, deriv_order)
...
Definition xc_derivatives.F:739
xc_input_constants
input constants for xc
Definition xc_input_constants.F:10
xc_input_constants::xc_pbe_orig
integer, parameter, public xc_pbe_orig
Definition xc_input_constants.F:41
xc_input_constants::xc_pbe_rev
integer, parameter, public xc_pbe_rev
Definition xc_input_constants.F:41
xc_input_constants::xc_pbe_sol
integer, parameter, public xc_pbe_sol
Definition xc_input_constants.F:41
xc_libxc
calculates a functional from libxc and its derivatives
Definition xc_libxc.F:29
xc_libxc::libxc_spin_unpolarized_eval
subroutine, public libxc_spin_unpolarized_eval(rho_set, deriv_set, grad_deriv, libxc_params, func_name_override)
evaluates the functional from libxc
Definition xc_libxc.F:887
xc_libxc::libxc_spin_polarized_eval
subroutine, public libxc_spin_polarized_eval(rho_set, deriv_set, grad_deriv, libxc_params, func_name_override)
evaluates the functional from libxc
Definition xc_libxc.F:1444
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_create
subroutine, public xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, tau_cutoff)
allocates and does (minimal) initialization of a rho_set
Definition xc_rho_set_types.F:111
xc_rho_set_types::xc_rho_set_release
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
Definition xc_rho_set_types.F:129
input_section_types::section_type
represent a section of the input file
Definition input_section_types.F:102
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78
evaluate
subroutine evaluate(x, energy, gradient, hessian, derivative_order)
Evaluate native PBE in its independent spin-density and gradient-norm variables.
Definition xc_semilocal_unittest.F:561
xc_semilocal_unittest
program xc_semilocal_unittest
Definition xc_semilocal_unittest.F:8