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 2 : PROGRAM xc_semilocal_unittest
9 2 : 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 2 : CALL cp2k_init()
62 2 : NULLIFY (specification, functional)
63 2 : CALL create_xc_fun_section(specification)
64 2 : CALL section_vals_create(functional, specification)
65 2 : CALL section_vals_val_set(functional, 'PBE%_SECTION_PARAMETERS_', l_val=.TRUE.)
66 2 : energy_error = 0.0_dp
67 2 : hessian_error = 0.0_dp
68 8 : DO variant = 1, SIZE(variants)
69 6 : CALL section_vals_val_set(functional, 'PBE%PARAMETRIZATION', i_val=variants(variant))
70 26 : DO configuration = 1, 3
71 60 : DO spin = 1, 2
72 36 : IF (configuration == 1 .AND. spin == 2) CYCLE
73 30 : x = [0.4_dp, 0.2_dp, 0.3_dp, 0.2_dp, 0.25_dp]
74 30 : IF (configuration == 2) x(spin) = -0.1_dp
75 30 : IF (configuration == 3) x(spin) = 1.e-13_dp
76 30 : CALL evaluate(x, energy, gradient, hessian)
77 198 : DO axis = 1, 5
78 150 : step = 1.e-5_dp
79 150 : IF (configuration == 3 .AND. axis == spin) step = 1.e-14_dp
80 150 : trial = x
81 150 : trial(axis) = x(axis) + step
82 150 : CALL evaluate(trial, plus, gp, scratch)
83 150 : trial(axis) = x(axis) - step
84 150 : CALL evaluate(trial, minus, gm, scratch)
85 150 : energy_error = MAX(energy_error, ABS((plus - minus)/(2.0_dp*step) - gradient(axis)))
86 936 : 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 2 : WRITE (*, '(A,2ES15.6)') 'PBE spin-floor first/second derivative errors:', energy_error, hessian_error
92 2 : IF (energy_error > 2.e-8_dp) ERROR STOP 'PBE potential is not the derivative of its clamped-density energy'
93 2 : IF (hessian_error > 2.e-7_dp) ERROR STOP 'PBE kernel is not the derivative of its clamped-density potential'
94 2 : CALL section_vals_val_set(functional, 'PBE%_SECTION_PARAMETERS_', l_val=.FALSE.)
95 2 : CALL section_vals_val_set(functional, 'TPSS%_SECTION_PARAMETERS_', l_val=.TRUE.)
96 2 : CALL check_tpss_zero_gradient()
97 2 : CALL section_vals_val_set(functional, 'TPSS%_SECTION_PARAMETERS_', l_val=.FALSE.)
98 2 : CALL section_vals_release(functional)
99 2 : CALL section_vals_create(functional, specification)
100 2 : CALL section_vals_val_set(functional, 'XWPBE%_SECTION_PARAMETERS_', l_val=.TRUE.)
101 2 : CALL check_xwpbe_limits()
102 2 : CALL section_vals_val_set(functional, 'XWPBE%_SECTION_PARAMETERS_', l_val=.FALSE.)
103 : #if defined (__LIBXC)
104 2 : CALL check_libxc_tau_projection()
105 2 : CALL check_libxc_c2_window()
106 2 : CALL check_libxc_c2_rounding()
107 2 : CALL check_libxc_individual_floor()
108 : #endif
109 2 : CALL section_vals_release(functional)
110 2 : CALL section_release(specification)
111 2 : CALL cp2k_finalize()
112 :
113 : CONTAINS
114 :
115 : #if defined (__LIBXC)
116 : ! **************************************************************************************************
117 : !> \brief Check meta-GGA potentials and kernels against independent differences of the projected energy.
118 : ! **************************************************************************************************
119 2 : 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 2 : REAL(dp), ALLOCATABLE :: gm(:), gp(:), gradient(:), &
129 2 : hessian(:, :), scratch(:, :), &
130 2 : trial(:), x(:)
131 :
132 2 : first_error = 0.0_dp
133 2 : second_error = 0.0_dp
134 6 : DO spin_case = 1, 2
135 4 : nvar = 3
136 4 : nconfig = 3
137 4 : IF (spin_case == 2) THEN
138 2 : nvar = 7
139 2 : nconfig = 6
140 : END IF
141 0 : ALLOCATE (x(nvar), trial(nvar), gradient(nvar), gp(nvar), gm(nvar), &
142 16 : hessian(nvar, nvar), scratch(nvar, nvar))
143 28 : DO family = 1, SIZE(names)
144 136 : DO configuration = 1, nconfig
145 108 : IF (spin_case == 1) THEN
146 144 : x(:) = [0.4_dp, 0.8_dp, 0.3_dp]
147 36 : IF (configuration == 2) x(3) = 0.02_dp
148 36 : IF (configuration == 3) x(3) = -0.02_dp
149 : ELSE
150 576 : x(:) = [0.3_dp, 0.2_dp, 0.6_dp, 0.4_dp, 0.3_dp, 0.12_dp, 0.10_dp]
151 96 : IF (configuration == 2) x(6:7) = [0.02_dp, 0.015_dp]
152 72 : IF (configuration == 3) x(6) = -0.02_dp
153 96 : IF (configuration == 4) x(6:7) = [-0.02_dp, -0.015_dp]
154 72 : IF (configuration == 5) x(1) = -0.1_dp
155 72 : IF (configuration == 6) x(2) = -0.1_dp
156 : END IF
157 108 : CALL evaluate_libxc_point(names(family), x, energy, gradient, hessian)
158 108 : IF (configuration >= 3) THEN
159 60 : IF (INDEX(names(family), 'MGGA_X_') == 1 .AND. ABS(energy) < 1.e-4_dp) THEN
160 0 : ERROR STOP 'Libxc discarded finite density at nonpositive tau'
161 : END IF
162 432 : trial(:) = x(:)
163 60 : trial(nvar - spin_case + 1) = 1.e-10_dp
164 60 : CALL evaluate_libxc_point(names(family), trial, plus, gp, scratch)
165 60 : trial(nvar - spin_case + 1) = -1.e-10_dp
166 60 : CALL evaluate_libxc_point(names(family), trial, minus, gm, scratch)
167 60 : IF (ABS(plus - minus) > 1.e-9_dp) ERROR STOP 'Libxc energy is discontinuous through zero tau'
168 : END IF
169 852 : DO axis = 1, nvar
170 4464 : trial(:) = x(:)
171 612 : trial(axis) = x(axis) + step
172 612 : CALL evaluate_libxc_point(names(family), trial, plus, gp, scratch)
173 612 : trial(axis) = x(axis) - step
174 612 : CALL evaluate_libxc_point(names(family), trial, minus, gm, scratch)
175 : first_error = MAX(first_error, &
176 612 : 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 4572 : (1.0_dp + ABS(hessian(:, axis)))))
180 : END DO
181 : END DO
182 : END DO
183 6 : DEALLOCATE (x, trial, gradient, gp, gm, hessian, scratch)
184 : END DO
185 2 : WRITE (*, '(A,2ES15.6)') 'Libxc projected-tau first/second derivative errors:', first_error, second_error
186 2 : IF (first_error > 3.e-8_dp) ERROR STOP 'Libxc projected-tau potential is not the energy derivative'
187 2 : IF (second_error > 2.e-6_dp) ERROR STOP 'Libxc projected-tau kernel is not the potential derivative'
188 4 : END SUBROUTINE check_libxc_tau_projection
189 :
190 : ! **************************************************************************************************
191 : !> \brief ...
192 : ! **************************************************************************************************
193 2 : 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 2 : REAL(dp), ALLOCATABLE :: gm(:), gmm(:), gp(:), gpp(:), gradient(:), hessian(:, :), hm(:, :), &
205 2 : hmm(:, :), hp(:, :), hpp(:, :), third(:, :, :), trial(:), x(:)
206 :
207 : ! Nonzero offsets keep the third-derivative differences well conditioned.
208 :
209 2 : first_error = 0.0_dp
210 2 : second_error = 0.0_dp
211 2 : third_error = 0.0_dp
212 6 : DO spin_case = 1, 2
213 4 : nvar = 3
214 4 : IF (spin_case == 2) nvar = 7
215 0 : ALLOCATE (x(nvar), trial(nvar), gradient(nvar), gm(nvar), gp(nvar), &
216 0 : gmm(nvar), gpp(nvar), hessian(nvar, nvar), &
217 24 : third(nvar, nvar, nvar), hm(nvar, nvar), hp(nvar, nvar), hmm(nvar, nvar), hpp(nvar, nvar))
218 20 : DO family = 1, SIZE(names)
219 68 : DO configuration = 1, SIZE(offsets)
220 48 : IF (spin_case == 1) THEN
221 96 : x(:) = [0.4_dp, 0.64_dp, 0.3_dp]
222 24 : wa = x(2)/(8.0_dp*x(1))
223 24 : widtha = relative_width*wa
224 24 : x(3) = wa + offsets(configuration)*widtha
225 : ELSE
226 192 : x(:) = [0.3_dp, 0.2_dp, 0.055_dp, 0.16_dp, 0.09_dp, 0.12_dp, 0.10_dp]
227 24 : wa = x(4)/(8.0_dp*x(1))
228 24 : wb = x(5)/(8.0_dp*x(2))
229 24 : widtha = relative_width*wa
230 24 : widthb = relative_width*wb
231 24 : x(6) = wa + offsets(configuration)*widtha
232 24 : x(7) = wb - offsets(configuration)*widthb
233 : END IF
234 : ! The installed r2SCAN functions expose only first and second derivatives.
235 48 : IF (family <= 2) THEN
236 24 : CALL evaluate_libxc_point(names(family), x, energy, gradient, hessian, third_all=third, contractions=.TRUE.)
237 : ELSE
238 24 : CALL evaluate_libxc_point(names(family), x, energy, gradient, hessian, contractions=.TRUE.)
239 : END IF
240 304 : DO axis = 1, nvar
241 240 : h = 0.01_dp*widtha
242 240 : IF (spin_case == 1) THEN
243 72 : IF (axis == 1) h = 0.01_dp*widtha/(wa/x(1))
244 72 : IF (axis == 2) h = 0.01_dp*widtha/(wa/x(2))
245 : ELSE
246 168 : IF (axis == 1) h = 0.01_dp*widtha/(wa/x(1))
247 168 : IF (axis == 2) h = 0.01_dp*widthb/(wb/x(2))
248 168 : IF (axis == 3) h = 1.e-5_dp
249 144 : IF (axis == 4) h = 0.01_dp*widtha/(wa/x(4))
250 168 : IF (axis == 5) h = 0.01_dp*widthb/(wb/x(5))
251 168 : IF (axis == 7) h = 0.01_dp*widthb
252 : END IF
253 1632 : trial(:) = x(:); trial(axis) = x(axis) + h
254 240 : CALL evaluate_libxc_point(names(family), trial, ep, gp, hp, contractions=.TRUE.)
255 240 : trial(axis) = x(axis) - h
256 240 : CALL evaluate_libxc_point(names(family), trial, em, gm, hm, contractions=.TRUE.)
257 240 : trial(axis) = x(axis) + 2.0_dp*h
258 240 : CALL evaluate_libxc_point(names(family), trial, epp, gpp, hpp, contractions=.TRUE.)
259 240 : trial(axis) = x(axis) - 2.0_dp*h
260 240 : CALL evaluate_libxc_point(names(family), trial, emm, gmm, hmm, contractions=.TRUE.)
261 240 : first_fd = (8.0_dp*(ep - em) - (epp - emm))/(12.0_dp*h)
262 240 : first_error = MAX(first_error, ABS(first_fd - gradient(axis))/(1.0_dp + ABS(gradient(axis))))
263 : second_error = MAX(second_error, &
264 1632 : MAXVAL(ABS((8.0_dp*(gp - gm) - (gpp - gmm))/(12.0_dp*h) - hessian(:, axis))/(1.0_dp + ABS(hessian(:, axis)))))
265 288 : 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 5256 : (1.0_dp + ABS(third(:, :, axis)))))
269 : END IF
270 : END DO
271 : END DO
272 : END DO
273 6 : DEALLOCATE (x, trial, gradient, gm, gp, gmm, gpp, hessian, third, hm, hp, hmm, hpp)
274 : END DO
275 2 : WRITE (*, '(A,3ES15.6)') 'Libxc C2-window first/second/third derivative errors:', first_error, second_error, third_error
276 2 : IF (first_error > 2.e-6_dp) ERROR STOP 'C2-window energy Jacobian failed'
277 2 : IF (second_error > 2.e-5_dp) ERROR STOP 'C2-window potential Hessian failed'
278 2 : IF (third_error > 2.e-5_dp) ERROR STOP 'C2-window kernel third derivatives failed'
279 2 : END SUBROUTINE check_libxc_c2_window
280 :
281 : ! **************************************************************************************************
282 : !> \brief ...
283 : ! **************************************************************************************************
284 2 : 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 2 : REAL(dp), ALLOCATABLE :: gp(:), gradient(:), hessian(:, :), &
294 2 : hp(:, :), probe(:), x(:)
295 :
296 2 : error = 0.0_dp
297 6 : DO spin_case = 1, 2
298 4 : nvar = 3
299 4 : IF (spin_case == 2) nvar = 7
300 16 : ALLOCATE (x(nvar), probe(nvar), gradient(nvar), gp(nvar), hessian(nvar, nvar), hp(nvar, nvar))
301 20 : DO family = 1, SIZE(names)
302 16 : IF (spin_case == 1) THEN
303 32 : x(:) = [0.4_dp, 0.8_dp, 0.3_dp]
304 8 : w(1) = x(2)**2/(8.0_dp*x(1))
305 : ELSE
306 64 : x(:) = [0.3_dp, 0.2_dp, 0.6_dp, 0.4_dp, 0.3_dp, 0.12_dp, 0.10_dp]
307 8 : w(1) = x(4)**2/(8.0_dp*x(1))
308 8 : w(2) = x(5)**2/(8.0_dp*x(2))
309 : END IF
310 96 : probe(:) = x(:)
311 40 : DO channel = 1, spin_case
312 24 : tau_axis = nvar - spin_case + channel
313 24 : width(channel) = relative_width*w(channel)
314 24 : x(tau_axis) = w(channel) + 0.99999_dp*width(channel)
315 : ! The nearby unprojected point supplies raw Libxc derivatives.
316 40 : probe(tau_axis) = w(channel) + 1.00001_dp*width(channel)
317 : END DO
318 16 : CALL evaluate_libxc_point(names(family), x, energy, gradient, hessian)
319 16 : CALL evaluate_libxc_point(names(family), probe, energy, gp, hp)
320 44 : DO channel = 1, spin_case
321 24 : tau_axis = nvar - spin_case + channel
322 24 : u = (x(tau_axis) - w(channel))/width(channel)
323 24 : slope = (1.0_dp + u)**2*(2.0_dp - u)/4.0_dp
324 24 : curve = 0.75_dp*(1.0_dp - u*u)/width(channel)
325 24 : expected = hp(tau_axis, tau_axis)*slope*slope + gp(tau_axis)*curve
326 40 : error = MAX(error, ABS(hessian(tau_axis, tau_axis) - expected)/(1.0_dp + ABS(expected)))
327 : END DO
328 : END DO
329 6 : DEALLOCATE (x, probe, gradient, gp, hessian, hp)
330 : END DO
331 2 : WRITE (*, '(A,ES15.6)') 'Libxc rounded-window curvature error:', error
332 2 : IF (error > 2.e-6_dp) ERROR STOP 'Rounded energy value discarded C2 projection curvature'
333 2 : END SUBROUTINE check_libxc_c2_rounding
334 :
335 : ! **************************************************************************************************
336 : !> \brief Check the original-density energy normalization below Libxc's individual spin floors.
337 : ! **************************************************************************************************
338 2 : 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 2 : error_first = 0.0_dp
347 2 : error_second = 0.0_dp
348 2 : error_third = 0.0_dp
349 14 : DO f = 1, SIZE(names)
350 86 : DO configuration = 1, 6
351 72 : x = [0.6_dp, 0.01_dp, 0.21_dp, 0.2_dp, 0.02_dp, 0.3_dp, 0.2_dp]
352 12 : SELECT CASE (configuration)
353 : CASE (2)
354 12 : x(2) = -0.01_dp
355 : CASE (3)
356 36 : x(1:2) = 0.03_dp
357 : CASE (4)
358 36 : x(1:2) = [0.01_dp, 0.6_dp]
359 : CASE (5)
360 12 : x(7) = -0.02_dp ! Compose the density floor with the kinetic-density projection.
361 : CASE (6)
362 108 : x(1:2) = [-0.01_dp, 0.6_dp]
363 : END SELECT
364 72 : IF (f <= 2) THEN
365 24 : CALL evaluate_libxc_point(TRIM(names(f)), x, energy, gradient, hessian, 0.05_dp, third)
366 : ELSE
367 48 : CALL evaluate_libxc_point(TRIM(names(f)), x, energy, gradient, hessian, 0.05_dp)
368 : END IF
369 588 : DO axis = 1, SIZE(x)
370 504 : step = 1.e-6_dp
371 504 : trial = x
372 504 : trial(axis) = x(axis) + step
373 504 : CALL evaluate_libxc_point(TRIM(names(f)), trial, plus, gp, hp, 0.05_dp)
374 504 : trial(axis) = x(axis) - step
375 504 : CALL evaluate_libxc_point(TRIM(names(f)), trial, minus, gm, hm, 0.05_dp)
376 504 : error_first = MAX(error_first, ABS((plus - minus)/(2.0_dp*step) - gradient(axis)))
377 4032 : error_second = MAX(error_second, MAXVAL(ABS((gp - gm)/(2.0_dp*step) - hessian(:, axis))))
378 576 : 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 336 : (2.0_dp*step) - third(:, :, axis))))
381 : END IF
382 : END DO
383 : END DO
384 : END DO
385 2 : WRITE (*, '(A,3ES15.6)') 'Libxc individual spin-floor first/second/third derivative errors:', &
386 4 : error_first, error_second, error_third
387 2 : IF (error_first > 2.e-7_dp) ERROR STOP 'Libxc individual spin-floor potential is not an energy derivative'
388 2 : IF (error_second > 2.e-6_dp) ERROR STOP 'Libxc individual spin-floor kernel is not a potential derivative'
389 2 : IF (error_third > 2.e-5_dp) ERROR STOP 'Libxc individual spin-floor third derivatives are inconsistent'
390 2 : 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 3596 : 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 3572 : 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 35720 : bounds = 1
425 3572 : use_contractions = .FALSE.
426 3572 : IF (PRESENT(contractions)) use_contractions = contractions
427 3572 : labels = 0
428 3572 : needs%norm_drho = .TRUE.
429 3572 : IF (SIZE(x) == 3) THEN
430 604 : nspins = 1
431 2416 : labels(1:3) = [deriv_rho, deriv_norm_drho, deriv_tau]
432 604 : needs%rho = .TRUE.
433 604 : needs%tau = .TRUE.
434 604 : needs%laplace_rho = .TRUE.
435 : ELSE
436 2968 : nspins = 2
437 2968 : labels = [deriv_rhoa, deriv_rhob, deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob, deriv_tau_a, deriv_tau_b]
438 2968 : needs%rho_spin = .TRUE.
439 2968 : needs%norm_drho_spin = .TRUE.
440 2968 : needs%tau_spin = .TRUE.
441 2968 : needs%laplace_rho_spin = .TRUE.
442 : END IF
443 3572 : IF (use_contractions) THEN
444 1008 : IF (nspins == 1) THEN
445 312 : labels(2) = deriv_gamma
446 : ELSE
447 2784 : labels(3:5) = [deriv_gamma_ab, deriv_gamma_aa, deriv_gamma_bb]
448 : END IF
449 : END IF
450 3572 : cutoff = 1.e-20_dp
451 3572 : IF (PRESENT(density_cutoff)) cutoff = density_cutoff
452 3572 : CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=cutoff, tau_cutoff=1.e-20_dp)
453 3572 : CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
454 3572 : IF (nspins == 1) THEN
455 2416 : rho_set%rho = x(1)
456 2416 : rho_set%norm_drho = x(2)
457 2416 : rho_set%tau = x(3)
458 2416 : rho_set%laplace_rho = 0.0_dp
459 1540 : IF (use_contractions) rho_set%norm_drho = SQRT(x(2))
460 : ELSE
461 11872 : rho_set%rhoa = x(1)
462 11872 : rho_set%rhob = x(2)
463 11872 : rho_set%norm_drho = x(3)
464 11872 : rho_set%norm_drhoa = x(4)
465 11872 : rho_set%norm_drhob = x(5)
466 11872 : rho_set%tau_a = x(6)
467 11872 : rho_set%tau_b = x(7)
468 11872 : rho_set%laplace_rhoa = 0.0_dp
469 11872 : rho_set%laplace_rhob = 0.0_dp
470 2968 : IF (use_contractions) THEN
471 2784 : rho_set%norm_drho = SQRT(2.0_dp*x(3) + x(4) + x(5))
472 2784 : rho_set%norm_drhoa = SQRT(x(4))
473 2784 : rho_set%norm_drhob = SQRT(x(5))
474 : END IF
475 : END IF
476 3572 : rho_set%has = needs
477 3572 : CALL xc_dset_create(derivatives, local_bounds=bounds)
478 26160 : DO i = 1, SIZE(x)
479 45176 : derivative => xc_dset_get_derivative(derivatives, [labels(i)], allocate_deriv=.TRUE.)
480 22588 : CALL xc_derivative_get(derivative, deriv_data=values)
481 90352 : values = 0.125_dp
482 112888 : DO j = i, SIZE(x)
483 260184 : derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j)], allocate_deriv=.TRUE.)
484 86728 : CALL xc_derivative_get(derivative, deriv_data=values)
485 369500 : values = 0.125_dp
486 : END DO
487 : END DO
488 3572 : order = 2
489 3572 : IF (PRESENT(third_rho) .OR. PRESENT(third_all)) THEN
490 48 : IF (PRESENT(third_rho) .AND. nspins /= 2) ERROR STOP 'Spin third-derivative oracle requires two spins'
491 48 : order = 3
492 336 : DO i = 1, SIZE(x)
493 2208 : DO j = 1, SIZE(x)
494 14832 : DO k = 1, SIZE(x)
495 50688 : derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j), labels(k)], allocate_deriv=.TRUE.)
496 12672 : CALL xc_derivative_get(derivative, deriv_data=values)
497 52560 : values = 0.125_dp
498 : END DO
499 : END DO
500 : END DO
501 : END IF
502 3572 : NULLIFY (parameters)
503 3572 : IF (nspins == 1) THEN
504 604 : CALL libxc_spin_unpolarized_eval(rho_set, derivatives, order, parameters, func_name_override=name)
505 : ELSE
506 2968 : CALL libxc_spin_polarized_eval(rho_set, derivatives, order, parameters, func_name_override=name)
507 : END IF
508 3572 : derivative => xc_dset_get_derivative(derivatives, [INTEGER::])
509 3572 : CALL xc_derivative_get(derivative, deriv_data=values)
510 3572 : energy = values(1, 1, 1)
511 26160 : DO i = 1, SIZE(x)
512 45176 : derivative => xc_dset_get_derivative(derivatives, [labels(i)])
513 22588 : CALL xc_derivative_get(derivative, deriv_data=values)
514 22588 : gradient(i) = values(1, 1, 1) - 0.125_dp
515 112888 : DO j = i, SIZE(x)
516 260184 : derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j)])
517 86728 : CALL xc_derivative_get(derivative, deriv_data=values)
518 86728 : hessian(i, j) = values(1, 1, 1) - 0.125_dp
519 109316 : hessian(j, i) = hessian(i, j)
520 : END DO
521 : END DO
522 3572 : IF (PRESENT(third_rho)) THEN
523 72 : DO i = 1, 2
524 168 : DO j = 1, 2
525 336 : DO k = 1, 2
526 768 : derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j), labels(k)])
527 192 : CALL xc_derivative_get(derivative, deriv_data=values)
528 288 : 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 3572 : IF (PRESENT(third_all)) THEN
534 144 : DO k = 1, SIZE(x)
535 840 : DO j = 1, SIZE(x)
536 5256 : DO i = 1, SIZE(x)
537 17760 : derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j), labels(k)])
538 4440 : CALL xc_derivative_get(derivative, deriv_data=values)
539 5136 : third_all(i, j, k) = values(1, 1, 1) - 0.125_dp
540 : END DO
541 : END DO
542 : END DO
543 5280 : IF (.NOT. ALL(ieee_is_finite(third_all))) ERROR STOP 'Nonfinite Libxc third derivatives'
544 : END IF
545 203188 : IF (.NOT. ieee_is_finite(energy) .OR. .NOT. ALL(ieee_is_finite(gradient)) .OR. &
546 0 : .NOT. ALL(ieee_is_finite(hessian))) ERROR STOP 'Nonfinite Libxc projected-tau result'
547 3572 : CALL xc_dset_release(derivatives)
548 3572 : CALL xc_rho_set_release(rho_set)
549 67892 : 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 892 : 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 446 : 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 446 : order = 2
573 116 : IF (PRESENT(derivative_order)) order = derivative_order
574 4460 : bounds = 1
575 446 : needs = xc_functionals_get_needs(functional, lsd=.TRUE., calc_potential=.TRUE.)
576 446 : needs%norm_drho = .TRUE.
577 446 : CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=1.e-20_dp)
578 446 : CALL xc_rho_set_atom_update(rho_set, needs, 2, bounds)
579 1784 : rho_set%rhoa = x(1)
580 1784 : rho_set%rhob = x(2)
581 1784 : rho_set%norm_drho = x(3)
582 1784 : rho_set%norm_drhoa = x(4)
583 1784 : rho_set%norm_drhob = x(5)
584 446 : rho_set%has%rho_spin = .TRUE.
585 446 : rho_set%has%norm_drho = .TRUE.
586 446 : rho_set%has%norm_drho_spin = .TRUE.
587 446 : CALL xc_dset_create(derivatives, local_bounds=bounds)
588 : ! Seed every first/second derivative to catch overwriting earlier functional contributions.
589 2676 : DO i = 1, 5
590 4460 : derivative => xc_dset_get_derivative(derivatives, [labels(i)], allocate_deriv=.TRUE.)
591 2230 : CALL xc_derivative_get(derivative, deriv_data=values)
592 8920 : values = 0.125_dp
593 9366 : DO j = i, 5
594 20070 : derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j)], allocate_deriv=.TRUE.)
595 6690 : CALL xc_derivative_get(derivative, deriv_data=values)
596 28990 : values = 0.125_dp
597 : END DO
598 : END DO
599 446 : CALL xc_functionals_eval(functional, lsd=.TRUE., rho_set=rho_set, deriv_set=derivatives, deriv_order=order)
600 446 : derivative => xc_dset_get_derivative(derivatives, [INTEGER::])
601 446 : CALL xc_derivative_get(derivative, deriv_data=values)
602 446 : energy = values(1, 1, 1)
603 2676 : DO i = 1, 5
604 4460 : derivative => xc_dset_get_derivative(derivatives, [labels(i)])
605 2230 : CALL xc_derivative_get(derivative, deriv_data=values)
606 2230 : gradient(i) = values(1, 1, 1) - 0.125_dp
607 9366 : DO j = i, 5
608 20070 : derivative => xc_dset_get_derivative(derivatives, [labels(i), labels(j)])
609 6690 : CALL xc_derivative_get(derivative, deriv_data=values)
610 6690 : hessian(i, j) = values(1, 1, 1) - 0.125_dp
611 6690 : IF (order == 1) hessian(i, j) = 0.0_dp
612 8920 : hessian(j, i) = hessian(i, j)
613 : END DO
614 : END DO
615 446 : CALL xc_dset_release(derivatives)
616 446 : CALL xc_rho_set_release(rho_set)
617 9812 : END SUBROUTINE evaluate
618 :
619 : ! **************************************************************************************************
620 : !> \brief Check native TPSS derivatives at stationary density points and below the von Weizsaecker bound.
621 : ! **************************************************************************************************
622 2 : 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 2 : error = 0.0_dp
630 2 : x(1) = 0.4_dp
631 2 : 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 8 : DO configuration = 1, SIZE(tau_scales)
633 6 : x(2) = 0.0_dp
634 6 : x(3) = tau_scales(configuration)*tau_unif
635 6 : CALL evaluate_tpss_point(x, energy, gradient)
636 6 : IF (ABS(gradient(2)) > 1.e-13_dp) ERROR STOP 'TPSS has a nonzero derivative at zero density gradient'
637 12 : DO axis = 1, 3, 2
638 12 : trial(:) = x(:)
639 12 : trial(axis) = x(axis) + step
640 12 : CALL evaluate_tpss_point(trial, plus, scratch)
641 12 : trial(axis) = x(axis) - step
642 12 : CALL evaluate_tpss_point(trial, minus, scratch)
643 12 : error = MAX(error, ABS((plus - minus)/(2.0_dp*step) - gradient(axis)))
644 : END DO
645 6 : trial = x
646 6 : trial(2) = 1.e-100_dp
647 6 : CALL evaluate_tpss_point(trial, plus, scratch)
648 6 : IF (ABS(plus - energy) > 1.e-13_dp .OR. ABS(scratch(2)) > 1.e-20_dp) THEN
649 0 : ERROR STOP 'TPSS does not approach its zero-gradient limit continuously'
650 : END IF
651 6 : x(2) = 1.e-4_dp
652 6 : CALL evaluate_tpss_point(x, energy, gradient)
653 6 : trial = x
654 6 : trial(2) = x(2) + step
655 6 : CALL evaluate_tpss_point(trial, plus, scratch)
656 6 : trial(2) = x(2) - step
657 6 : CALL evaluate_tpss_point(trial, minus, scratch)
658 20 : 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 2 : x = [0.4_dp, 1.2_dp, 0.05_dp]
662 2 : CALL evaluate_tpss_point(x, energy, gradient)
663 8 : DO axis = 1, 3
664 6 : trial = x
665 6 : trial(axis) = x(axis) + step
666 6 : CALL evaluate_tpss_point(trial, plus, scratch)
667 6 : trial(axis) = x(axis) - step
668 6 : CALL evaluate_tpss_point(trial, minus, scratch)
669 8 : 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 6 : DO configuration = 1, 2
673 4 : x = [0.4_dp, 1.2_dp, -0.05_dp]
674 4 : IF (configuration == 2) x(2) = 0.0_dp
675 4 : CALL evaluate_tpss_point(x, energy, gradient)
676 4 : IF (ABS(energy) < 1.e-4_dp .OR. gradient(3) /= 0.0_dp) THEN
677 0 : ERROR STOP 'TPSS negative-tau continuation is not finite and constant in tau'
678 : END IF
679 22 : DO axis = 1, 3
680 12 : IF (axis == 2 .AND. configuration == 2) CYCLE
681 10 : trial = x
682 10 : trial(axis) = x(axis) + step
683 10 : CALL evaluate_tpss_point(trial, plus, scratch)
684 10 : trial(axis) = x(axis) - step
685 10 : CALL evaluate_tpss_point(trial, minus, scratch)
686 16 : error = MAX(error, ABS((plus - minus)/(2.0_dp*step) - gradient(axis)))
687 : END DO
688 : END DO
689 2 : x = [0.4_dp, 1.2_dp, 0.0_dp]
690 2 : CALL evaluate_tpss_point(x, energy, gradient)
691 2 : trial = x
692 2 : trial(3) = 1.e-10_dp
693 2 : CALL evaluate_tpss_point(trial, plus, scratch)
694 2 : trial(3) = -1.e-10_dp
695 2 : CALL evaluate_tpss_point(trial, minus, scratch)
696 2 : IF (MAX(ABS(plus - energy), ABS(minus - energy)) > 1.e-9_dp) THEN
697 0 : ERROR STOP 'TPSS energy is discontinuous through zero kinetic density'
698 : END IF
699 2 : WRITE (*, '(A,ES15.6)') 'TPSS limiting-density first-derivative error:', error
700 2 : IF (error > 3.e-8_dp) ERROR STOP 'TPSS limiting-density potential is not the energy derivative'
701 4 : END SUBROUTINE check_tpss_zero_gradient
702 :
703 : ! **************************************************************************************************
704 : !> \brief Check screened exchange limits and gradient scaling against energy differences.
705 : ! **************************************************************************************************
706 2 : 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 2 : error = 0.0_dp
716 2 : CALL section_vals_val_set(functional, 'XWPBE%SCALE_X', r_val=-0.25_dp)
717 2 : CALL section_vals_val_set(functional, 'XWPBE%SCALE_X0', r_val=1.0_dp)
718 6 : DO range_case = 1, 2
719 4 : CALL section_vals_val_set(functional, 'XWPBE%OMEGA', r_val=0.11_dp*REAL(range_case - 1, dp))
720 18 : DO configuration = 1, SIZE(reduced_gradients)
721 12 : xr = [0.4_dp, 0.0_dp, 0.5_dp]
722 12 : 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 12 : CALL evaluate_tpss_point(xr, energy, gr)
724 36 : DO axis = 1, 2
725 24 : IF (axis == 2 .AND. configuration == 1) CYCLE
726 20 : trialr = xr; trialr(axis) = xr(axis) + step
727 20 : CALL evaluate_tpss_point(trialr, plus, scratchr)
728 20 : trialr(axis) = xr(axis) - step
729 20 : CALL evaluate_tpss_point(trialr, minus, scratchr)
730 36 : error = MAX(error, ABS((plus - minus)/(2.0_dp*step) - gr(axis))/(1.0_dp + ABS(gr(axis))))
731 : END DO
732 12 : IF (configuration == 1) THEN
733 4 : IF (ABS(gr(2)) > 1.e-13_dp) ERROR STOP 'XWPBE zero-gradient potential is not stationary'
734 4 : trialr = xr; trialr(2) = 1.e-8_dp
735 4 : CALL evaluate_tpss_point(trialr, plus, scratchr)
736 4 : IF (ABS(plus - energy) > 1.e-12_dp) ERROR STOP 'XWPBE zero-gradient energy is discontinuous'
737 : END IF
738 12 : xs = [0.3_dp, 0.2_dp, 0.3_dp, 0.0_dp, 0.0_dp]
739 12 : 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 12 : 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 12 : CALL evaluate(xs, energy, gs, hessian, derivative_order=1)
742 72 : DO axis = 1, 5
743 60 : IF (axis >= 4 .AND. configuration == 1) CYCLE
744 52 : trials = xs; trials(axis) = xs(axis) + step
745 52 : CALL evaluate(trials, plus, scratchs, hessian, derivative_order=1)
746 52 : trials(axis) = xs(axis) - step
747 52 : CALL evaluate(trials, minus, scratchs, hessian, derivative_order=1)
748 72 : error = MAX(error, ABS((plus - minus)/(2.0_dp*step) - gs(axis))/(1.0_dp + ABS(gs(axis))))
749 : END DO
750 52 : 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 2 : WRITE (*, '(A,ES15.6)') 'XWPBE limiting/large-gradient first-derivative error:', error
754 2 : IF (error > 3.e-8_dp) ERROR STOP 'XWPBE potential is not the derivative of its scaled-gradient energy'
755 2 : 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 308 : 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 154 : 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 1540 : bounds = 1
777 154 : needs = xc_functionals_get_needs(functional, lsd=.FALSE., calc_potential=.TRUE.)
778 154 : needs%tau = .TRUE.
779 154 : CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=1.e-20_dp, tau_cutoff=1.e-20_dp)
780 154 : CALL xc_rho_set_atom_update(rho_set, needs, 1, bounds)
781 616 : rho_set%rho = x(1)
782 616 : rho_set%norm_drho = x(2)
783 616 : rho_set%tau = x(3)
784 154 : rho_set%has%rho = .TRUE.
785 154 : rho_set%has%norm_drho = .TRUE.
786 154 : rho_set%has%tau = .TRUE.
787 154 : CALL xc_dset_create(derivatives, local_bounds=bounds)
788 616 : DO i = 1, SIZE(labels)
789 924 : derivative => xc_dset_get_derivative(derivatives, [labels(i)], allocate_deriv=.TRUE.)
790 462 : CALL xc_derivative_get(derivative, deriv_data=values)
791 2002 : values = 0.125_dp
792 : END DO
793 154 : CALL xc_functionals_eval(functional, lsd=.FALSE., rho_set=rho_set, deriv_set=derivatives, deriv_order=1)
794 154 : derivative => xc_dset_get_derivative(derivatives, [INTEGER::])
795 154 : CALL xc_derivative_get(derivative, deriv_data=values)
796 154 : energy = values(1, 1, 1)
797 616 : DO i = 1, SIZE(labels)
798 924 : derivative => xc_dset_get_derivative(derivatives, [labels(i)])
799 462 : CALL xc_derivative_get(derivative, deriv_data=values)
800 616 : gradient(i) = values(1, 1, 1) - 0.125_dp
801 : END DO
802 616 : IF (.NOT. ieee_is_finite(energy) .OR. .NOT. ALL(ieee_is_finite(gradient))) THEN
803 0 : ERROR STOP 'TPSS energy or first derivative is nonfinite at a stationary density point'
804 : END IF
805 154 : CALL xc_dset_release(derivatives)
806 154 : CALL xc_rho_set_release(rho_set)
807 3388 : END SUBROUTINE evaluate_tpss_point
808 :
809 : END PROGRAM xc_semilocal_unittest
|
