(git:6344619)
Loading...
Searching...
No Matches
spectral_localizer_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
8MODULE spectral_localizer_test_cases
9 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
10
11 IMPLICIT NONE
12 PRIVATE
13 PUBLIC :: check_localizer_2d, check_localizer_3d, check_localizer_flattening
14CONTAINS
15
16! **************************************************************************************************
17!> \brief Independent finite-model, metric, orientation and Pfaffian tests.
18!> \author Thomas D. Kuehne, tkuehne@cp2k.org
19! **************************************************************************************************
20 SUBROUTINE check_localizer_2d()
21#if defined(__HAS_IEEE_EXCEPTIONS)
22 USE ieee_exceptions, ONLY: ieee_all, ieee_get_halting_mode
23#endif
24 USE kinds, ONLY: dp
25 USE spectral_localizer, ONLY: localizer_aii_skew, &
26 localizer_evaluate, &
27 localizer_inertia, &
28 localizer_periodic_evaluate, &
29 localizer_pfaffian_sign
30
31 IMPLICIT NONE
32#if defined(__HAS_IEEE_EXCEPTIONS)
33 LOGICAL :: original_halt(size(ieee_all)), current_halt(size(ieee_all))
34#endif
35 COMPLEX(KIND=dp), ALLOCATABLE :: h(:, :), s(:, :), r(:, :, :), transform(:, :), &
36 hc(:, :), sc(:, :), rc(:, :, :), hs(:, :), ss(:, :), rs(:, :, :)
37 INTEGER :: axis, i, index, reference, status, n, pf, permutation(4), io
38 REAL(kind=dp) :: gap, gap_reference, ratio, skew(4, 4), exact
39
40 io = cp_logger_get_default_io_unit()
41#if defined(__HAS_IEEE_EXCEPTIONS)
42 CALL ieee_get_halting_mode(ieee_all, original_halt)
43#endif
44 CALL inertia_tests()
45 CALL real_inertia_tests()
46 CALL gap_tests()
47 CALL pfaffian_orientation_tests()
48 CALL aii_transform_tests()
49 CALL aii_roundoff_tests()
50 CALL periodic_tests()
51 ! Independent Pfaffian polynomial for size four, including odd permutations.
52 skew = 0.0_dp
53 skew(1, 2) = 0.3_dp
54 skew(1, 3) = 2.0_dp
55 skew(1, 4) = -1.0_dp
56 skew(2, 3) = 0.7_dp
57 skew(2, 4) = 3.0_dp
58 skew(3, 4) = 0.8_dp
59 skew = skew - transpose(skew)
60 exact = skew(1, 2)*skew(3, 4) - skew(1, 3)*skew(2, 4) + skew(1, 4)*skew(2, 3)
61 CALL localizer_pfaffian_sign(skew, pf, status)
62 IF (status /= 0 .OR. pf /= int(sign(1.0_dp, exact))) error stop 'Pfaffian polynomial failed'
63 permutation = [2, 1, 3, 4]
64 CALL localizer_pfaffian_sign(skew(permutation, permutation), pf, status)
65 IF (status /= 0 .OR. pf /= -int(sign(1.0_dp, exact))) error stop 'Pfaffian orientation failed'
66 skew = 0.0_dp
67 CALL localizer_pfaffian_sign(skew, pf, status)
68 IF (status == 0 .OR. pf /= 0) error stop 'Singular Pfaffian accepted'
69
70 ALLOCATE (h(1, 1), s(1, 1), r(1, 1, 2))
71 h = 4.0_dp
72 s = 2.0_dp
73 r(1, 1, :) = [6.0_dp, 8.0_dp]
74 CALL evaluate(h, s, r, [0.0_dp, 0.0_dp], 1.0_dp, 0.2_dp, .false., index, gap, status)
75 IF (status /= 0 .OR. index /= 0 .OR. abs(gap - sqrt(2.0_dp)) > 1.e-12_dp) THEN
76 error stop 'Analytic generalized localizer failed'
77 END IF
78 CALL evaluate(h, s, r, [3.0_dp, 4.0_dp], 2.0_dp, 0.2_dp, .false., index, gap, status)
79 IF (status /= 1) error stop 'Zero gap was assigned an index'
80 DEALLOCATE (h, s, r)
81
82 ALLOCATE (h(2, 2), s(2, 2), r(2, 2, 2))
83 h = 0.0_dp
84 s = 0.0_dp
85 r = 0.0_dp
86 s(1, 1) = 1.e-9_dp
87 s(2, 2) = 1.0_dp
88 h(1, 1) = 1.e-16_dp
89 h(2, 2) = 1.0_dp
90 CALL evaluate(h, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.2_dp, .false., index, gap, status)
91 IF (status /= 1) error stop 'Metric-amplified roundoff ignored'
92 DEALLOCATE (h, s, r)
93
94 CALL qwz_model(7, -1.0_dp, h, s, r)
95 n = SIZE(h, 1)
96 CALL evaluate(h, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .false., reference, gap_reference, status)
97 IF (status /= 0 .OR. reference /= -1 .OR. abs(gap_reference - 0.740327333453389_dp) > 1.e-11_dp) THEN
98 error stop 'Nontrivial open QWZ model failed'
99 END IF
100 CALL evaluate(h, s, r(:, :, [2, 1]), [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .false., index, gap, status)
101 IF (status /= 0 .OR. index /= -reference .OR. abs(gap - gap_reference) > 1.e-11_dp) THEN
102 error stop 'Oriented-plane reversal failed'
103 END IF
104 CALL evaluate(conjg(h), s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .false., index, gap, status)
105 IF (status /= 0 .OR. index /= -reference) error stop 'Time reversed Chern sign failed'
106 CALL evaluate(h, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.15_dp, .false., index, gap, status)
107 IF (status /= 0 .OR. index /= reference) error stop 'Kappa plateau failed'
108 CALL evaluate(h, s, r, [20.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .false., index, gap, status)
109 IF (status /= 0 .OR. index /= 0) error stop 'Outside-vacuum limit failed'
110
111 ALLOCATE (transform(n, n), hc(n, n), sc(n, n), rc(n, n, 2))
112 transform = 0.0_dp
113 DO i = 1, n
114 transform(i, i) = (1.0_dp + real(mod(i, 5), dp)/5.0_dp)*exp(cmplx(0.0_dp, 0.17_dp*i, dp))
115 IF (i < n) transform(i, i + 1) = cmplx(0.05_dp, 0.03_dp, dp)
116 END DO
117 hc(:, :) = matmul(conjg(transpose(transform)), matmul(h, transform))
118 sc(:, :) = matmul(conjg(transpose(transform)), matmul(s, transform))
119 DO axis = 1, 2
120 rc(:, :, axis) = matmul(conjg(transpose(transform)), matmul(r(:, :, axis), transform))
121 END DO
122 CALL evaluate(hc, sc, rc, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .false., index, gap, status)
123 IF (status /= 0 .OR. index /= reference .OR. abs(gap - gap_reference) > 1.e-11_dp) THEN
124 error stop 'Complex nonorthogonal basis invariance failed'
125 END IF
126 rc(:, :, 1) = rc(:, :, 1) + 13.0_dp*sc
127 rc(:, :, 2) = rc(:, :, 2) - 7.0_dp*sc
128 hc(:, :) = hc + 11.0_dp*sc
129 CALL evaluate(hc, sc, rc, [13.0_dp, -7.0_dp], 11.0_dp, 0.25_dp, .false., index, gap, status)
130 IF (status /= 0 .OR. index /= reference .OR. abs(gap - gap_reference) > 1.e-11_dp) THEN
131 error stop 'Origin/energy shift covariance failed'
132 END IF
133 sc = 0.0_dp
134 DO i = 1, n
135 sc(i, i) = 1.0_dp
136 END DO
137 sc(1, 1) = 1.e-14_dp
138 CALL evaluate(hc, sc, rc, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .false., index, gap, status)
139 IF (status /= -2) error stop 'Ill-conditioned metric accepted'
140 sc(1, 1) = -1.0_dp
141 CALL evaluate(hc, sc, rc, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .false., index, gap, status)
142 IF (status /= -2) error stop 'Indefinite metric accepted'
143 hc(1, 2) = hc(1, 2) + 1.0_dp
144 CALL evaluate(hc, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .false., index, gap, status)
145 IF (status /= -1) error stop 'Non-Hermitian Hamiltonian accepted'
146
147 ! Two time-reversed Chern copies, with TR-preserving, spin-nonconserving mixing.
148 ALLOCATE (hs(2*n, 2*n), ss(2*n, 2*n), rs(2*n, 2*n, 2))
149 hs = 0.0_dp
150 ss = 0.0_dp
151 rs = 0.0_dp
152 hs(:n, :n) = h
153 hs(n + 1:, n + 1:) = conjg(h)
154 ss(:n, :n) = s
155 ss(n + 1:, n + 1:) = s
156 rs(:n, :n, :) = r
157 rs(n + 1:, n + 1:, :) = r
158 DO i = 1, n, 2
159 hs(i, n + i + 1) = 0.1_dp
160 hs(i + 1, n + i) = -0.1_dp
161 END DO
162 hs(n + 1:, :n) = conjg(transpose(hs(:n, n + 1:)))
163 CALL evaluate(hs, ss, rs, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .true., index, gap, status)
164 IF (status /= 0 .OR. index /= 1) error stop 'Spin-mixed nontrivial AII model failed'
165 gap_reference = gap
166 ! Real, nonorthogonal AO change applied to both spin components preserves fixed Theta.
167 transform = real(transform, dp)
168 DO i = 1, n
169 transform(i, i) = 1.0_dp + real(mod(i, 5), dp)/5.0_dp
170 END DO
171 DO axis = 1, 2
172 DO i = 1, 2
173 hs((i - 1)*n + 1:i*n, (axis - 1)*n + 1:axis*n) = &
174 matmul(transpose(transform), matmul(hs((i - 1)*n + 1:i*n, (axis - 1)*n + 1:axis*n), transform))
175 END DO
176 END DO
177 sc(:, :) = matmul(transpose(transform), matmul(s, transform))
178 ss(:n, :n) = sc
179 ss(n + 1:, n + 1:) = sc
180 DO axis = 1, 2
181 rc(:, :, axis) = matmul(transpose(transform), matmul(r(:, :, axis), transform))
182 rs(:n, :n, axis) = rc(:, :, axis)
183 rs(n + 1:, n + 1:, axis) = rc(:, :, axis)
184 END DO
185 CALL evaluate(hs, ss, rs, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .true., index, gap, status)
186 IF (status /= 0 .OR. index /= 1 .OR. abs(gap - gap_reference) > 1.e-11_dp) THEN
187 error stop 'AII metric covariance failed'
188 END IF
189 hs(1, 1) = hs(1, 1) + 1.e-12_dp
190 CALL evaluate(hs, ss, rs, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .true., index, gap, status)
191 IF (status /= 0 .OR. index /= 1) error stop 'Roundoff-level AII residual failed'
192 hs(1, 1) = hs(1, 1) + 1.0_dp
193 CALL evaluate(hs, ss, rs, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .true., index, gap, status)
194 IF (status /= -3) error stop 'Broken physical time reversal accepted'
195 hs(:, :) = ss
196 CALL evaluate(hs, ss, rs, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .true., index, gap, status)
197 IF (status /= 0 .OR. index /= 0) error stop 'Atomic AII calibration failed'
198
199 DEALLOCATE (h, s, r)
200 CALL qwz_model(7, -3.0_dp, h, s, r)
201 CALL evaluate(h, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .false., index, gap, status)
202 IF (status /= 0 .OR. index /= 0) error stop 'Trivial QWZ model failed'
203#if defined(__HAS_IEEE_EXCEPTIONS)
204 CALL ieee_get_halting_mode(ieee_all, current_halt)
205 IF (any(current_halt .NEQV. original_halt)) error stop 'Caller floating-point traps changed'
206#endif
207
208 CONTAINS
209
210! **************************************************************************************************
211!> \brief Signed-permutation AII conversion against the explicit dense unitary transform.
212! **************************************************************************************************
213 SUBROUTINE aii_transform_tests()
214 COMPLEX(dp), ALLOCATABLE :: matrix(:, :), q(:, :), reference(:, :)
215 INTEGER :: i, j, m, n, status
216 REAL(dp) :: error, largest, residual
217 REAL(dp), ALLOCATABLE :: actual(:, :), expected(:, :)
218
219 largest = 0.0_dp
220 DO n = 1, 13
221 m = 4*n
222 ALLOCATE (q(m, m), matrix(m, m), reference(m, m), expected(m, m), actual(m, m))
223 q = 0.0_dp
224 DO i = 1, m
225 q(i, i) = 1.0_dp
226 DO j = 1, m
227 expected(i, j) = sin(real(i + 2*j, dp)) + cos(real(3*i - j, dp))
228 END DO
229 END DO
230 expected(:, :) = expected - transpose(expected)
231 DO i = 1, n
232 q(i, 3*n + i) = cmplx(0.0_dp, -1.0_dp, dp)
233 q(n + i, 2*n + i) = cmplx(0.0_dp, 1.0_dp, dp)
234 q(2*n + i, n + i) = cmplx(0.0_dp, 1.0_dp, dp)
235 q(3*n + i, i) = cmplx(0.0_dp, -1.0_dp, dp)
236 END DO
237 q = q/sqrt(2.0_dp)
238 matrix(:, :) = cmplx(0.0_dp, -1.0_dp, dp)*matmul(q, matmul(expected, conjg(transpose(q))))
239 reference(:, :) = matmul(conjg(transpose(q)), matmul(matrix, q))
240 CALL localizer_aii_skew(matrix, actual, 1.e-10_dp, residual, status)
241 error = max(maxval(abs(actual - expected)), maxval(abs(actual + aimag(reference))))
242 IF (status /= 0 .OR. error > 1.e-13_dp) error stop 'AII sparse-unitary transform failed'
243 largest = max(largest, error)
244 DO i = 1, m
245 matrix(i, i) = matrix(i, i) - cmplx(0.0_dp, 5.e-12_dp, dp)
246 END DO
247 CALL localizer_aii_skew(matrix, actual, 1.e-10_dp, residual, status)
248 IF (status /= 0 .OR. maxval(abs(actual - expected)) > 1.e-13_dp .OR. residual < 4.9e-12_dp) THEN
249 error stop 'AII skew projection or discarded-error bound failed'
250 END IF
251 IF (maxval(abs(actual + transpose(actual))) /= 0.0_dp) error stop 'AII projection is not exactly skew'
252 DO i = 1, m
253 matrix(i, i) = matrix(i, i) + cmplx(0.0_dp, 5.e-12_dp, dp)
254 END DO
255 matrix(1, 1) = matrix(1, 1) + 0.1_dp
256 CALL localizer_aii_skew(matrix, actual, 1.e-10_dp, residual, status)
257 IF (status == 0) error stop 'Broken AII transform accepted'
258 CALL localizer_aii_skew(matrix(:m - 1, :m - 1), actual, 1.e-10_dp, residual, status)
259 IF (status == 0) error stop 'Invalid AII dimensions accepted'
260 DEALLOCATE (q, matrix, reference, expected, actual)
261 END DO
262 IF (io > 0) WRITE (io, '(A,ES14.5)') 'AII signed-permutation maximum error: ', largest
263 END SUBROUTINE aii_transform_tests
264
265! **************************************************************************************************
266!> \brief Accepted Hermiticity roundoff must not fail the stricter Pfaffian skew check.
267! **************************************************************************************************
268 SUBROUTINE aii_roundoff_tests()
269 COMPLEX(dp) :: h0(4, 4), r0(4, 4, 2), s0(4, 4)
270 INTEGER :: i, index, status
271 REAL(dp) :: gap
272
273 h0 = 0.0_dp
274 s0 = 0.0_dp
275 r0 = 0.0_dp
276 DO i = 1, 4
277 s0(i, i) = 1.0_dp
278 h0(i, i) = real(3*mod(i + 1, 2) - 1, dp)
279 END DO
280 h0(1, 2) = 5.e-12_dp
281 h0(3, 4) = 5.e-12_dp
282 CALL evaluate(h0, s0, r0, [0.0_dp, 0.0_dp], 0.0_dp, 0.2_dp, .true., index, gap, status)
283 IF (status /= 0 .OR. index /= 0 .OR. abs(gap - 1.0_dp) > 1.e-12_dp) error stop 'AII roundoff rejected'
284 h0(1, 1) = -1.e-10_dp
285 h0(3, 3) = -1.e-10_dp
286 CALL evaluate(h0, s0, r0, [0.0_dp, 0.0_dp], 0.0_dp, 0.2_dp, .true., index, gap, status)
287 IF (status /= 1) error stop 'AII unresolved gap assigned an index'
288 h0(1, 2) = 1.e-5_dp
289 CALL evaluate(h0, s0, r0, [0.0_dp, 0.0_dp], 0.0_dp, 0.2_dp, .true., index, gap, status)
290 IF (status /= -1) error stop 'AII non-Hermitian input projected away'
291 END SUBROUTINE aii_roundoff_tests
292
293! **************************************************************************************************
294!> \brief Torus models with known Chern/Z2 phases, origin periodicity and metric covariance.
295! **************************************************************************************************
296 SUBROUTINE periodic_tests()
297 INTEGER, PARAMETER :: length = 8
298
299 COMPLEX(dp), ALLOCATABLE :: co(:, :, :), ct(:, :, :), h0(:, :), ht(:, :), pos(:, :, :), &
300 s0(:, :), si(:, :, :), sn(:, :, :), spin_c(:, :, :), spin_h(:, :), spin_s(:, :), &
301 spin_sin(:, :, :), st(:, :), t(:, :)
302 INTEGER :: col, expected, i, idx, info, j, kx, ky, &
303 n, orbital, phase, row, x, y
304 REAL(dp) :: cond, g, g0, mass, two_pi, wave(2)
305
306 two_pi = 8.0_dp*atan(1.0_dp)
307 DO phase = 1, 3
308 mass = -1.0_dp
309 expected = -1
310 IF (phase == 2) THEN
311 mass = 1.0_dp
312 expected = 1
313 ELSE IF (phase == 3) THEN
314 mass = -3.0_dp
315 expected = 0
316 END IF
317 CALL qwz_model(length, mass, h0, s0, pos, torus=.true.)
318 n = SIZE(h0, 1)
319 ALLOCATE (co(n, n, 2), si(n, n, 2))
320 co = 0.0_dp
321 si = 0.0_dp
322 DO j = 1, 2
323 DO i = 1, n
324 co(i, i, j) = cos(two_pi*real(pos(i, i, j), dp)/length)
325 si(i, i, j) = sin(two_pi*real(pos(i, i, j), dp)/length)
326 END DO
327 END DO
328 CALL localizer_periodic_evaluate(h0, s0, co, si, [0.0_dp, 0.0_dp], 0.0_dp, 1.0_dp, .false., &
329 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g0, cond, info)
330 IF (info /= 0 .OR. idx /= expected) error stop 'Periodic QWZ Chern phase failed'
331 CALL localizer_periodic_evaluate(h0, s0, co, si, [two_pi, -2*two_pi], 0.0_dp, 1.0_dp, .false., &
332 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
333 IF (info /= 0 .OR. idx /= expected .OR. abs(g - g0) > 1.e-11_dp) THEN
334 error stop 'Torus origin periodicity failed'
335 END IF
336 CALL localizer_periodic_evaluate(h0, s0, co, si, two_pi*[1.0_dp, 2.0_dp]/length, &
337 0.0_dp, 1.0_dp, .false., &
338 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
339 IF (info /= 0 .OR. idx /= expected .OR. abs(g - g0) > 1.e-11_dp) THEN
340 error stop 'Periodic lattice translation covariance failed'
341 END IF
342 CALL localizer_periodic_evaluate(h0, s0, co(:, :, [2, 1]), si(:, :, [2, 1]), &
343 [0.0_dp, 0.0_dp], 0.0_dp, 1.0_dp, .false., &
344 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
345 IF (info /= 0 .OR. idx /= -expected .OR. abs(g - g0) > 1.e-11_dp) THEN
346 error stop 'Periodic orientation reversal failed'
347 END IF
348
349 ALLOCATE (spin_h(2*n, 2*n), spin_s(2*n, 2*n), spin_c(2*n, 2*n, 2), spin_sin(2*n, 2*n, 2))
350 spin_h = 0.0_dp
351 spin_s = 0.0_dp
352 spin_c = 0.0_dp
353 spin_sin = 0.0_dp
354 spin_h(:n, :n) = h0
355 spin_h(n + 1:, n + 1:) = conjg(h0)
356 DO i = 1, n, 2
357 spin_h(i, n + i + 1) = 0.1_dp
358 spin_h(i + 1, n + i) = -0.1_dp
359 END DO
360 spin_h(n + 1:, :n) = conjg(transpose(spin_h(:n, n + 1:)))
361 spin_s(:n, :n) = s0
362 spin_s(n + 1:, n + 1:) = s0
363 spin_c(:n, :n, :) = co
364 spin_c(n + 1:, n + 1:, :) = co
365 spin_sin(:n, :n, :) = si
366 spin_sin(n + 1:, n + 1:, :) = si
367 CALL localizer_periodic_evaluate(spin_h, spin_s, spin_c, spin_sin, [0.0_dp, 0.0_dp], &
368 0.0_dp, 1.0_dp, .true., &
369 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
370 IF (info /= 0 .OR. idx /= abs(expected)) error stop 'Periodic spin-mixed Z2 phase failed'
371 spin_h(:, :) = spin_s
372 CALL localizer_periodic_evaluate(spin_h, spin_s, spin_c, spin_sin, [0.0_dp, 0.0_dp], &
373 0.0_dp, 1.0_dp, .true., &
374 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
375 IF (info /= 0 .OR. idx /= 0) error stop 'Periodic atomic Pfaffian calibration failed'
376 DEALLOCATE (spin_h, spin_s, spin_c, spin_sin)
377
378 IF (phase == 1) THEN
379 ALLOCATE (t(n, n), ht(n, n), st(n, n), ct(n, n, 2), sn(n, n, 2))
380 ! Full Fourier transform: H becomes block diagonal, but torus coordinates do not.
381 t = 0.0_dp
382 DO ky = 0, length - 1
383 DO kx = 0, length - 1
384 col = 2*(ky*length + kx) + 1
385 DO y = 0, length - 1
386 DO x = 0, length - 1
387 row = 2*(y*length + x) + 1
388 DO orbital = 0, 1
389 t(row + orbital, col + orbital) = &
390 exp(cmplx(0.0_dp, two_pi*(x*kx + y*ky)/length, dp))/length
391 END DO
392 END DO
393 END DO
394 END DO
395 END DO
396 ht(:, :) = matmul(conjg(transpose(t)), matmul(h0, t))
397 st(:, :) = s0
398 DO ky = 0, length - 1
399 DO kx = 0, length - 1
400 col = 2*(ky*length + kx) + 1
401 wave = two_pi*[real(kx, dp), real(ky, dp)]/length
402 IF (abs(ht(col, col) - (mass + sum(cos(wave)))) > 1.e-12_dp) THEN
403 error stop 'Periodic real-space/k-space Hamiltonian mismatch'
404 END IF
405 IF (abs(ht(col, col + 1) - cmplx(sin(wave(1)), -sin(wave(2)), dp)) > 1.e-12_dp) THEN
406 error stop 'Periodic Bloch Hamiltonian offdiagonal mismatch'
407 END IF
408 END DO
409 END DO
410 DO j = 1, 2
411 ct(:, :, j) = matmul(conjg(transpose(t)), matmul(co(:, :, j), t))
412 sn(:, :, j) = matmul(conjg(transpose(t)), matmul(si(:, :, j), t))
413 END DO
414 CALL localizer_periodic_evaluate(ht, st, ct, sn, [0.0_dp, 0.0_dp], &
415 0.0_dp, 1.0_dp, .false., &
416 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
417 IF (info /= 0 .OR. idx /= expected .OR. abs(g - g0) > 1.e-11_dp) THEN
418 error stop 'Full k-mesh/real-space periodic localizer mismatch'
419 END IF
420 t = 0.0_dp
421 DO i = 1, n
422 t(i, i) = (1.0_dp + mod(i, 5)/5.0_dp)*exp(cmplx(0.0_dp, 0.17_dp*i, dp))
423 IF (i < n) t(i, i + 1) = cmplx(0.07_dp, -0.04_dp, dp)
424 END DO
425 ht(:, :) = matmul(conjg(transpose(t)), matmul(h0, t))
426 st(:, :) = matmul(conjg(transpose(t)), t)
427 DO j = 1, 2
428 ct(:, :, j) = matmul(conjg(transpose(t)), matmul(co(:, :, j), t))
429 sn(:, :, j) = matmul(conjg(transpose(t)), matmul(si(:, :, j), t))
430 END DO
431 CALL localizer_periodic_evaluate(ht + 7.0_dp*st, st, ct, sn, [0.0_dp, 0.0_dp], &
432 7.0_dp, 1.0_dp, .false., &
433 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
434 IF (info /= 0 .OR. idx /= expected .OR. abs(g - g0) > 1.e-11_dp) THEN
435 error stop 'Periodic metric/energy covariance failed'
436 END IF
437 CALL localizer_periodic_evaluate(2*ht, st, ct, sn, [0.0_dp, 0.0_dp], &
438 0.0_dp, 2.0_dp, .false., &
439 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
440 IF (info /= 0 .OR. idx /= expected .OR. abs(g - 2*g0) > 1.e-11_dp) THEN
441 error stop 'Periodic energy-scale covariance failed'
442 END IF
443 DEALLOCATE (t, ht, st, ct, sn)
444 END IF
445 CALL localizer_periodic_evaluate(h0, s0, co, si, [0.0_dp, 0.0_dp], 0.0_dp, 0.0_dp, .false., &
446 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
447 IF (info /= -1) error stop 'Zero torus energy scale accepted'
448 co(1, 2, 1) = co(1, 2, 1) + 0.1_dp
449 CALL localizer_periodic_evaluate(h0, s0, co, si, [0.0_dp, 0.0_dp], 0.0_dp, 1.0_dp, .false., &
450 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
451 IF (info /= -1) error stop 'Non-Hermitian torus cosine accepted'
452 h0 = 0.0_dp
453 si = 0.0_dp
454 co(:, :, 1) = s0
455 co(:, :, 2) = s0
456 CALL localizer_periodic_evaluate(h0, s0, co, si, [0.0_dp, 0.0_dp], 0.0_dp, 1.0_dp, .false., &
457 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
458 IF (info /= 1) error stop 'Closed periodic localizer assigned an index'
459 DEALLOCATE (h0, s0, pos, co, si)
460 END DO
461 END SUBROUTINE periodic_tests
462
463! **************************************************************************************************
464!> \brief Exercise 2x2 pivots, complex congruences, scale changes and singularity.
465! **************************************************************************************************
466 SUBROUTINE inertia_tests()
467 COMPLEX(KIND=dp) :: change(4, 4), matrix(4, 4), &
468 transformed(4, 4)
469 INTEGER :: exponent, i, j, m, p, RESULT
470
471 matrix = 0.0_dp
472 matrix(1, 2) = cmplx(0.0_dp, 2.0_dp, dp)
473 matrix(2, 1) = conjg(matrix(1, 2))
474 matrix(3, 3) = 1.e-12_dp
475 matrix(4, 4) = 3.0_dp
476 change = 0.0_dp
477 DO i = 1, 4
478 change(i, i) = cmplx(1.0_dp + 0.2_dp*i, 0.1_dp*i, dp)
479 DO j = i + 1, 4
480 change(i, j) = cmplx(0.03_dp*j, -0.07_dp*i, dp)
481 END DO
482 END DO
483 transformed = matmul(conjg(transpose(change)), matmul(matrix, change))
484 DO exponent = -150, 150, 150
485 CALL localizer_inertia((10.0_dp**exponent)*transformed, p, m, result)
486 IF (result /= 0 .OR. p /= 3 .OR. m /= 1) error stop 'LDL block inertia or congruence failed'
487 END DO
488 matrix(3, 3) = -1.e-12_dp
489 CALL localizer_inertia(matrix, p, m, result)
490 IF (result /= 0 .OR. p /= 2 .OR. m /= 2) error stop 'LDL gap crossing failed'
491 matrix(3, 3) = 0.0_dp
492 CALL localizer_inertia(matrix, p, m, result)
493 IF (result /= 1) error stop 'Singular LDL accepted'
494 matrix(2, 1) = matrix(1, 2)
495 CALL localizer_inertia(matrix, p, m, result)
496 IF (result /= -1) error stop 'Complex symmetric matrix accepted as Hermitian'
497 matrix = 0.0_dp
498 matrix(1, 2) = cmplx(0.0_dp, 1.e-12_dp, dp)
499 matrix(2, 1) = conjg(matrix(1, 2))
500 matrix(3, 3) = 1.0_dp
501 matrix(4, 4) = -1.0_dp
502 CALL localizer_inertia(matrix, p, m, result)
503 IF (result /= 0 .OR. p /= 2 .OR. m /= 2) error stop 'Small imaginary coupling discarded'
504 END SUBROUTINE inertia_tests
505
506! **************************************************************************************************
507!> \brief Cross-check real and complex LDL under congruences, all signatures and permutations.
508! **************************************************************************************************
509 SUBROUTINE real_inertia_tests()
510 COMPLEX(dp), ALLOCATABLE :: cmatrix(:, :), phases(:)
511 INTEGER :: exponent, i, info, j, m, n, negative, p
512 REAL(dp) :: pair(2, 2), scale
513 REAL(dp), ALLOCATABLE :: change(:, :), diagonal(:, :), &
514 matrix(:, :)
515
516 DO n = 1, 12
517 ALLOCATE (diagonal(n, n), change(n, n), matrix(n, n), phases(n), cmatrix(n, n))
518 change = 0.0_dp
519 DO i = 1, n
520 change(i, i) = 0.5_dp + real(i, dp)/n
521 phases(i) = exp(cmplx(0.0_dp, 0.37_dp*i, dp))
522 DO j = i + 1, n
523 change(i, j) = 0.1_dp*sin(real(i + 3*j, dp))
524 END DO
525 END DO
526 DO negative = 0, n
527 diagonal = 0.0_dp
528 DO i = 1, n
529 diagonal(i, i) = 1.0_dp + real(i, dp)/n
530 IF (i <= negative) diagonal(i, i) = -diagonal(i, i)
531 END DO
532 matrix(:, :) = matmul(transpose(change), matmul(diagonal, change))
533 DO i = 1, n
534 DO j = 1, n
535 cmatrix(i, j) = conjg(phases(i))*matrix(i, j)*phases(j)
536 END DO
537 END DO
538 DO exponent = -150, 150, 150
539 scale = 10.0_dp**exponent
540 CALL localizer_inertia(scale*matrix, p, m, info)
541 IF (info /= 0 .OR. p /= n - negative .OR. m /= negative) error stop 'Real LDL signature failed'
542 CALL localizer_inertia(scale*cmatrix, p, m, info)
543 IF (info /= 0 .OR. p /= n - negative .OR. m /= negative) error stop 'Complex LDL signature failed'
544 CALL localizer_inertia(cmplx(scale*matrix, 0.0_dp, dp), p, m, info)
545 IF (info /= 0 .OR. p /= n - negative .OR. m /= negative) error stop 'Real LDL dispatch failed'
546 END DO
547 END DO
548 DEALLOCATE (diagonal, change, matrix, phases, cmatrix)
549 END DO
550 pair = reshape([0.0_dp, 2.0_dp, 2.0_dp, 0.0_dp], [2, 2])
551 CALL localizer_inertia(pair, p, m, info)
552 IF (info /= 0 .OR. p /= 1 .OR. m /= 1) error stop 'Zero-diagonal real 2x2 pivot failed'
553 pair(1, 2) = 0.0_dp
554 CALL localizer_inertia(pair, p, m, info)
555 IF (info /= -1) error stop 'Nonsymmetric real matrix accepted'
556 pair = 0.0_dp
557 CALL localizer_inertia(pair, p, m, info)
558 IF (info /= 1) error stop 'Singular real matrix accepted'
559 END SUBROUTINE real_inertia_tests
560
561! **************************************************************************************************
562!> \brief Generalized gaps have an analytic reference independent of LDL pivot magnitudes.
563! **************************************************************************************************
564 SUBROUTINE gap_tests()
565 COMPLEX(dp) :: change(2, 2), h(2, 2), hc(2, 2), &
566 r(2, 2, 2), s(2, 2), sc(2, 2)
567 INTEGER :: i, index, status
568 REAL(dp) :: crossing(5), expected, gap, ratio, trace
569
570 h = reshape([1.0_dp, 20.0_dp, 20.0_dp, 399.0_dp], [2, 2])
571 s = 0.0_dp
572 s(1, 1) = 2.0_dp
573 s(2, 2) = 2.0_dp
574 r = 0.0_dp
575 trace = real(h(1, 1) + h(2, 2), dp)
576 expected = 1.0_dp/(trace + sqrt(trace**2 + 4.0_dp))
577 CALL localizer_evaluate(h, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 1.0_dp, .false., &
578 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, index, gap, ratio, status)
579 IF (status /= 0 .OR. index /= 0 .OR. abs(gap - expected) > 1.e-12_dp) error stop 'Analytic LDL/gap separation failed'
580 change = reshape([0.5_dp, 0.0_dp, 0.2_dp, 3.0_dp], [2, 2])
581 hc = matmul(transpose(change), matmul(h, change))
582 sc = matmul(transpose(change), matmul(s, change))
583 CALL localizer_evaluate(hc, sc, r, [0.0_dp, 0.0_dp], 0.0_dp, 1.0_dp, .false., &
584 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, index, gap, ratio, status)
585 IF (status /= 0 .OR. abs(gap - expected) > 1.e-12_dp) error stop 'Metric LDL/gap separation failed'
586 crossing = [-1.e-6_dp, -1.e-12_dp, 0.0_dp, 1.e-12_dp, 1.e-6_dp]
587 h = 0.0_dp
588 h(2, 2) = 5.0_dp
589 DO i = 1, SIZE(crossing)
590 h(1, 1) = crossing(i)
591 CALL localizer_evaluate(h, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 1.0_dp, .false., &
592 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, index, gap, ratio, status)
593 IF (i == 1 .OR. i == 5) THEN
594 IF (status /= 0 .OR. abs(gap - 5.e-7_dp) > 1.e-14_dp) error stop 'Resolved gap crossing failed'
595 ELSE
596 IF (status /= 1) error stop 'Unresolved gap assigned an index'
597 END IF
598 END DO
599 END SUBROUTINE gap_tests
600
601! **************************************************************************************************
602!> \brief Independent Pfaffian polynomials distinguish signs despite identical determinants.
603! **************************************************************************************************
604 SUBROUTINE pfaffian_orientation_tests()
605 INTEGER :: exponent, i, j, n, pf_sign, sample, &
606 status
607 INTEGER, ALLOCATABLE :: permutation(:)
608 REAL(dp) :: exact, scale
609 REAL(dp), ALLOCATABLE :: a(:, :), change(:, :), transformed(:, :)
610
611 DO n = 2, 8, 2
612 ALLOCATE (a(n, n), change(n, n), transformed(n, n), permutation(n))
613 DO sample = 1, 8
614 a = 0.0_dp
615 DO j = 2, n
616 DO i = 1, j - 1
617 a(i, j) = sin(real(11*i + 17*j + sample, dp))
618 END DO
619 END DO
620 DO i = 1, n, 2
621 a(i, i + 1) = a(i, i + 1) + 3.0_dp
622 END DO
623 a(:, :) = a - transpose(a)
624 exact = pfaffian_polynomial(a)
625 IF (abs(exact) < 1.e-6_dp) error stop 'Ill-conditioned Pfaffian fixture'
626 permutation(:) = [(i, i=1, n)]
627 permutation(:2) = [2, 1]
628 DO exponent = -150, 150, 150
629 scale = 10.0_dp**exponent
630 CALL localizer_pfaffian_sign(scale*a, pf_sign, status)
631 IF (status /= 0 .OR. pf_sign /= int(sign(1.0_dp, exact))) error stop 'Pfaffian polynomial sign failed'
632 CALL localizer_pfaffian_sign(-scale*a, pf_sign, status)
633 IF (status /= 0 .OR. pf_sign /= (-1)**(n/2)*int(sign(1.0_dp, exact))) error stop 'Negative Pfaffian scale failed'
634 CALL localizer_pfaffian_sign(scale*a(permutation, permutation), pf_sign, status)
635 IF (status /= 0 .OR. pf_sign /= -int(sign(1.0_dp, exact))) error stop 'Equal determinant/opposite Pfaffian failed'
636 END DO
637 change = 0.0_dp
638 DO i = 1, n
639 change(i, i) = 0.8_dp + 0.1_dp*i
640 IF (i < n) change(i, i + 1) = 0.1_dp
641 END DO
642 change(1, :) = -change(1, :)
643 transformed(:, :) = matmul(change, matmul(a, transpose(change)))
644 CALL localizer_pfaffian_sign(transformed, pf_sign, status)
645 IF (status /= 0 .OR. pf_sign /= -int(sign(1.0_dp, exact))) error stop 'Pfaffian congruence orientation failed'
646 END DO
647 DEALLOCATE (a, change, transformed, permutation)
648 END DO
649 ALLOCATE (a(4, 4))
650 a = 0.0_dp
651 a(1, 2) = 1.0_dp
652 a(1, 3) = 1.0_dp
653 a(2, 4) = 1.0_dp
654 a(3, 4) = 1.0_dp
655 a(:, :) = a - transpose(a)
656 CALL localizer_pfaffian_sign(a, pf_sign, status)
657 IF (status /= -2 .OR. pf_sign /= 0) error stop 'Pfaffian polynomial cancellation ignored'
658 END SUBROUTINE pfaffian_orientation_tests
659
660! **************************************************************************************************
661!> \brief Small independent recursive Pfaffian oracle, without pivoting or elimination.
662!> \param a Even-order skew matrix, at most eight rows in these tests
663!> \return Pfaffian from the expansion along the first row
664! **************************************************************************************************
665 RECURSIVE FUNCTION pfaffian_polynomial(a) RESULT(value)
666 REAL(dp), INTENT(IN) :: a(:, :)
667 REAL(dp) :: value
668
669 INTEGER :: i, j, n, keep(SIZE(a, 1) - 2)
670
671 n = SIZE(a, 1)
672 IF (n == 2) THEN
673 value = a(1, 2)
674 RETURN
675 END IF
676 value = 0.0_dp
677 DO j = 2, n
678 keep = pack([(i, i=2, n)], [(i /= j, i=2, n)])
679 value = value + real((-1)**j, dp)*a(1, j)*pfaffian_polynomial(a(keep, keep))
680 END DO
681 END FUNCTION pfaffian_polynomial
682
683! **************************************************************************************************
684!> \brief Apply uniform numerical tolerances in tests.
685!> \param h ...
686!> \param s ...
687!> \param r ...
688!> \param point ...
689!> \param energy ...
690!> \param kappa ...
691!> \param aii ...
692!> \param invariant ...
693!> \param gap ...
694!> \param status ...
695! **************************************************************************************************
696 SUBROUTINE evaluate(h, s, r, point, energy, kappa, aii, invariant, gap, status)
697 COMPLEX(KIND=dp), INTENT(IN) :: h(:, :), s(:, :), r(:, :, :)
698 REAL(KIND=dp), INTENT(IN) :: point(2), energy, kappa
699 LOGICAL, INTENT(IN) :: aii
700 INTEGER, INTENT(OUT) :: invariant
701 REAL(KIND=dp), INTENT(OUT) :: gap
702 INTEGER, INTENT(OUT) :: status
703
704 CALL localizer_evaluate(h, s, r, point, energy, kappa, aii, 1.e-10_dp, 1.e-10_dp, &
705 1.e-9_dp, invariant, gap, ratio, status)
706 END SUBROUTINE evaluate
707
708! **************************************************************************************************
709!> \brief Open square Qi-Wu-Zhang model: sin(kx)sx+sin(ky)sy+(m+cos(kx)+cos(ky))sz.
710!> \param length Number of sites along each direction
711!> \param mass Onsite mass; -1 nontrivial, -3 trivial in the infinite model
712!> \param h Hamiltonian in site-major, two-orbital order
713!> \param s Identity metric
714!> \param r Commuting diagonal positions centered on the flake
715!> \param torus Optional periodic boundary hoppings
716! **************************************************************************************************
717 SUBROUTINE qwz_model(length, mass, h, s, r, torus)
718 INTEGER, INTENT(IN) :: length
719 REAL(KIND=dp), INTENT(IN) :: mass
720 COMPLEX(KIND=dp), ALLOCATABLE, INTENT(OUT) :: h(:, :), s(:, :), r(:, :, :)
721 LOGICAL, INTENT(IN), OPTIONAL :: torus
722
723 COMPLEX(KIND=dp) :: tx(2, 2), ty(2, 2)
724 INTEGER :: i, j, n, orbital, x, y
725
726 n = 2*length**2
727 ALLOCATE (h(n, n), s(n, n), r(n, n, 2))
728 h = 0.0_dp
729 s = 0.0_dp
730 r = 0.0_dp
731 tx = 0.0_dp
732 tx(1, 1) = 0.5_dp
733 tx(2, 2) = -0.5_dp
734 tx(1, 2) = cmplx(0.0_dp, -0.5_dp, dp)
735 tx(2, 1) = tx(1, 2)
736 ty = 0.0_dp
737 ty(1, 1) = 0.5_dp
738 ty(2, 2) = -0.5_dp
739 ty(1, 2) = -0.5_dp
740 ty(2, 1) = 0.5_dp
741 DO y = 1, length
742 DO x = 1, length
743 i = 2*((y - 1)*length + x) - 1
744 h(i, i) = mass
745 h(i + 1, i + 1) = -mass
746 DO orbital = 0, 1
747 s(i + orbital, i + orbital) = 1.0_dp
748 r(i + orbital, i + orbital, 1) = real(x, dp) - real(length + 1, dp)/2.0_dp
749 r(i + orbital, i + orbital, 2) = real(y, dp) - real(length + 1, dp)/2.0_dp
750 END DO
751 IF (x < length) THEN
752 j = i + 2
753 h(i:i + 1, j:j + 1) = tx
754 h(j:j + 1, i:i + 1) = conjg(transpose(tx))
755 END IF
756 IF (y < length) THEN
757 j = i + 2*length
758 h(i:i + 1, j:j + 1) = ty
759 h(j:j + 1, i:i + 1) = conjg(transpose(ty))
760 END IF
761 IF (PRESENT(torus)) THEN
762 IF (torus .AND. x == length) THEN
763 j = 2*(y - 1)*length + 1
764 h(i:i + 1, j:j + 1) = tx
765 h(j:j + 1, i:i + 1) = conjg(transpose(tx))
766 END IF
767 IF (torus .AND. y == length) THEN
768 j = 2*x - 1
769 h(i:i + 1, j:j + 1) = ty
770 h(j:j + 1, i:i + 1) = conjg(transpose(ty))
771 END IF
772 END IF
773 END DO
774 END DO
775 END SUBROUTINE qwz_model
776
777 END SUBROUTINE check_localizer_2d
778
779! **************************************************************************************************
780!> \brief Independent chiral singular values, strong/weak/trivial Wilson-Dirac models and AO covariance.
781!> \author Thomas D. Kuehne, tkuehne@cp2k.org
782! **************************************************************************************************
783 SUBROUTINE check_localizer_3d()
784 USE ieee_arithmetic, ONLY: ieee_value, ieee_quiet_nan
785#if defined(__HAS_IEEE_EXCEPTIONS)
786 USE ieee_exceptions, ONLY: ieee_all, ieee_get_halting_mode, ieee_set_halting_mode
787#endif
788 USE kinds, ONLY: dp
789 USE spectral_localizer, ONLY: localizer_3d_evaluate, localizer_3d_matrix_evaluate
790 IMPLICIT NONE
791 COMPLEX(dp), PARAMETER :: ii = cmplx(0.0_dp, 1.0_dp, dp)
792 INTEGER, PARAMETER :: side = 3, sites = side**3, n = 4*sites
793 REAL(dp), PARAMETER :: masses(3) = [0.0_dp, 4.0_dp, 2.0_dp]
794 COMPLEX(dp) :: h(n, n), s(n, n), r(n, n, 3), u(n, n), hc(n, n), sc(n, n), rc(n, n, 3)
795 COMPLEX(dp) :: pauli(2, 2, 3), a(2*n, 2*n), dummy(1, 1), work(10*n)
796 REAL(dp) :: singular(2*n), rwork(10*n), gap, reference, ratio, mass, shift(3), defect, lower, upper, midpoint
797 INTEGER :: i, j, axis, model, index, expected, status, sign_product, kx, ky, kz, io
798#if defined(__HAS_IEEE_EXCEPTIONS)
799 LOGICAL :: halt(size(ieee_all)), current_halt(size(ieee_all))
800 CALL ieee_get_halting_mode(ieee_all, halt)
801#endif
802
803 io = cp_logger_get_default_io_unit()
804 CALL determinant_controls()
805 pauli = 0.0_dp
806 pauli(1, 2, 1) = 1.0_dp
807 pauli(2, 1, 1) = 1.0_dp
808 pauli(1, 2, 2) = -ii
809 pauli(2, 1, 2) = ii
810 pauli(1, 1, 3) = 1.0_dp
811 pauli(2, 2, 3) = -1.0_dp
812 DO model = 1, 3
813 mass = masses(model)
814 CALL dirac_model(mass)
815 ! Bulk inversion eigenvalues at the eight TRIM: strong, not weak, parity.
816 sign_product = 1
817 DO kx = -1, 1, 2
818 DO ky = -1, 1, 2
819 DO kz = -1, 1, 2
820 sign_product = sign_product*int(sign(1.0_dp, mass + kx + ky + kz))
821 END DO
822 END DO
823 END DO
824 expected = (1 - sign_product)/2
825 CALL evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
826 IF (status /= 0 .OR. index /= expected) error stop '3D strong/trivial parity mismatch'
827 reference = gap
828 a = 0.0_dp
829 a(:n, :n) = h
830 a(n + 1:, n + 1:) = h
831 DO axis = 1, 3
832 DO j = 1, 2
833 DO i = 1, 2
834 a((i - 1)*n + 1:i*n, (j - 1)*n + 1:j*n) = &
835 a((i - 1)*n + 1:i*n, (j - 1)*n + 1:j*n) - 0.5_dp*ii*pauli(i, j, axis)*r(:, :, axis)
836 END DO
837 END DO
838 END DO
839#if defined(__HAS_IEEE_EXCEPTIONS)
840 CALL ieee_set_halting_mode(ieee_all, .false.)
841#endif
842 CALL zgesvd('N', 'N', 2*n, 2*n, a, 2*n, singular, dummy, 1, dummy, 1, work, SIZE(work), rwork, status)
843#if defined(__HAS_IEEE_EXCEPTIONS)
844 CALL ieee_set_halting_mode(ieee_all, halt)
845#endif
846 IF (status /= 0 .OR. abs(reference - minval(singular)) > 2.e-12_dp) error stop '3D singular gap mismatch'
847 IF (io > 0) WRITE (io, '(A,F4.1,I3,2ES18.9)') '3D Dirac mass, index, gap, SVD error: ', &
848 mass, expected, reference, abs(reference - minval(singular))
849 END DO
850
851 ! Reuse the validated mass-2 model for the covariance and origin checks.
852 u = 0.0_dp
853 DO i = 1, n/2
854 u(i, i) = cmplx(0.8_dp + 0.004_dp*i, 0.07_dp*sin(real(i, dp)), dp)
855 u(i + n/2, i + n/2) = conjg(u(i, i))
856 u(i, i + n/2) = cmplx(0.08_dp, 0.03_dp, dp)
857 u(i + n/2, i) = -conjg(u(i, i + n/2))
858 IF (i < n/2) THEN
859 u(i, i + 1) = cmplx(0.02_dp, 0.01_dp, dp)
860 u(i + n/2, i + n/2 + 1) = conjg(u(i, i + 1))
861 END IF
862 END DO
863 hc = matmul(conjg(transpose(u)), matmul(h, u))
864 sc = matmul(conjg(transpose(u)), matmul(s, u))
865 DO axis = 1, 3
866 rc(:, :, axis) = matmul(conjg(transpose(u)), matmul(r(:, :, axis), u))
867 END DO
868 CALL evaluate(hc, sc, rc, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
869 defect = abs(gap - reference)
870 IF (status /= 0 .OR. index /= 1 .OR. defect > 3.e-12_dp) error stop '3D nonorthogonal spin-mixing basis'
871 IF (io > 0) WRITE (io, '(A,ES18.9)') '3D nonorthogonal spin-mixing gap residual: ', defect
872 shift = [1.2_dp, -0.8_dp, 0.4_dp]
873 DO axis = 1, 3
874 rc(:, :, axis) = r(:, :, axis) + shift(axis)*s
875 END DO
876 CALL evaluate(h + 3.0_dp*s, s, rc, shift, 3.0_dp)
877 IF (status /= 0 .OR. index /= 1 .OR. abs(gap - reference) > 3.e-12_dp) error stop '3D origin shifts'
878 rc(:, :, 1) = 0.6_dp*r(:, :, 1) + 0.8_dp*r(:, :, 2)
879 rc(:, :, 2) = -0.8_dp*r(:, :, 1) + 0.6_dp*r(:, :, 2)
880 rc(:, :, 3) = -r(:, :, 3)
881 CALL evaluate(h, s, rc, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
882 IF (status /= 0 .OR. index /= 1 .OR. abs(gap - reference) > 3.e-12_dp) error stop '3D reflection/rotation'
883 ! A too-small finite box/scale is not automatically classified by the bulk parity.
884 CALL localizer_3d_evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp, 0.3_dp, &
885 1.e-10_dp, 1.e-10_dp, 1.e-10_dp, index, gap, ratio, status)
886 IF (status /= 0 .OR. index /= 0) error stop '3D finite-size counterexample lost'
887 lower = 0.3_dp
888 upper = 0.5_dp
889 DO i = 1, 32
890 midpoint = 0.5_dp*(lower + upper)
891 CALL localizer_3d_evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp, midpoint, &
892 1.e-10_dp, 1.e-10_dp, 0.0_dp, index, gap, ratio, status)
893 IF (status == 1) EXIT
894 IF (status /= 0) error stop '3D transition factorization failed'
895 IF (index == 0) THEN
896 lower = midpoint
897 ELSE
898 upper = midpoint
899 END IF
900 END DO
901 CALL localizer_3d_evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp, midpoint, &
902 1.e-10_dp, 1.e-10_dp, 1.e-8_dp, index, gap, ratio, status)
903 IF (status /= 1) error stop '3D closing gap assigned an index'
904 IF (io > 0) WRITE (io, '(A,2ES18.9)') '3D localizer transition scale and unresolved gap: ', midpoint, gap
905 CALL localizer_3d_evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp, 0.5_dp, &
906 1.e-10_dp, -1.0_dp, 1.e-10_dp, index, gap, ratio, status)
907 IF (status /= -1) error stop '3D invalid metric tolerance accepted'
908 CALL localizer_3d_evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp, 0.5_dp, &
909 1.e-10_dp, 1.e-10_dp, ieee_value(0.0_dp, ieee_quiet_nan), index, gap, ratio, status)
910 IF (status /= -1) error stop '3D nonfinite gap tolerance accepted'
911 h(1, 1) = h(1, 1) + 0.1_dp
912 CALL evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
913 IF (status /= -3) error stop '3D broken time reversal accepted'
914 h = 0.0_dp
915 r = 0.0_dp
916 CALL evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
917 IF (status /= 1 .OR. index /= 0) error stop '3D zero gap accepted'
918 s(1, 1) = -1.0_dp
919 s(n/2 + 1, n/2 + 1) = -1.0_dp
920 CALL evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
921 IF (status /= -2) error stop '3D indefinite metric accepted'
922 h(1, 1) = ieee_value(0.0_dp, ieee_quiet_nan)
923 CALL evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
924 IF (status /= -1) error stop '3D nonfinite operator accepted'
925#if defined(__HAS_IEEE_EXCEPTIONS)
926 CALL ieee_get_halting_mode(ieee_all, current_halt)
927 IF (any(halt .NEQV. current_halt)) error stop '3D localizer changed IEEE exception modes'
928#endif
929 IF (io > 0) WRITE (io, *) '3D AII determinant, generalized gap, basis, origin, TR and finite-size tests passed.'
930 CONTAINS
931! **************************************************************************************************
932!> \brief Chiral-block determinant controls independent of the Wilson-Dirac model.
933!> \note Wilkinson's matrix has exponential LU pivot growth despite a large singular gap.
934! **************************************************************************************************
935 SUBROUTINE determinant_controls()
936 INTEGER, PARAMETER :: m = 64, quarter = m/4
937
938 COMPLEX(dp) :: block(m, m), metric(m/2, m/2), &
939 unused(1, 1), workspace(5*m)
940 INTEGER :: col, col_sign, info, orientation, &
941 paired_col, paired_row, RESULT, row, &
942 row_sign
943 REAL(dp) :: matrix(m, m), measured_gap, &
944 metric_condition, real_work(5*m), &
945 svd_gap, values(m)
946
947 metric = 0.0_dp
948 DO row = 1, m/2
949 metric(row, row) = 1.0_dp
950 END DO
951 matrix = 0.0_dp
952 DO row = 1, m
953 matrix(row, row) = 1.0_dp
954 matrix(row, :row - 1) = -1.0_dp
955 matrix(row, m) = 1.0_dp + 0.25_dp*modulo(row, 3)
956 END DO
957 DO orientation = 0, 1
958 ! Positive last-column entries give a positive determinant. One row reversal changes its sign.
959 IF (orientation == 1) matrix(1, :) = -matrix(1, :)
960 ! Exact dyadic entries in Q W Q^dagger retain the pivot-growth fixture after conversion.
961 DO col = 1, m
962 paired_col = (3 - (col - 1)/quarter)*quarter + modulo(col - 1, quarter) + 1
963 col_sign = merge(1, -1, col <= quarter .OR. col > 3*quarter)
964 DO row = 1, m
965 paired_row = (3 - (row - 1)/quarter)*quarter + modulo(row - 1, quarter) + 1
966 row_sign = merge(1, -1, row <= quarter .OR. row > 3*quarter)
967 block(row, col) = 0.5_dp*(matrix(row, col) + row_sign*col_sign*matrix(paired_row, paired_col) + &
968 ii*(col_sign*matrix(row, paired_col) - row_sign*matrix(paired_row, col)))
969 END DO
970 END DO
971 CALL localizer_3d_matrix_evaluate(block, metric, 1.e-10_dp, 1.e-10_dp, 1.e-10_dp, &
972 result, measured_gap, metric_condition, info)
973 IF (info /= 0 .OR. result /= orientation) error stop '3D high-growth determinant sign'
974#if defined(__HAS_IEEE_EXCEPTIONS)
975 CALL ieee_set_halting_mode(ieee_all, .false.)
976#endif
977 CALL zgesvd('N', 'N', m, m, block, m, values, unused, 1, unused, 1, workspace, &
978 SIZE(workspace), real_work, info)
979#if defined(__HAS_IEEE_EXCEPTIONS)
980 CALL ieee_set_halting_mode(ieee_all, halt)
981#endif
982 svd_gap = minval(values)
983 IF (info /= 0 .OR. abs(measured_gap - svd_gap) > 1.e-12_dp) error stop '3D high-growth singular gap'
984 IF (io > 0) WRITE (io, '(A,I2,2ES18.9)') '3D high-growth determinant, gap, SVD error: ', &
985 result, measured_gap, abs(measured_gap - svd_gap)
986 END DO
987 ! Independent atomic limit for both signs of a spin-degenerate energy.
988 DO col = -1, 1, 2
989 block = 0.0_dp
990 DO row = 1, m
991 block(row, row) = real(col, dp)
992 END DO
993 CALL localizer_3d_matrix_evaluate(block, metric, 1.e-10_dp, 1.e-10_dp, 1.e-10_dp, &
994 result, measured_gap, metric_condition, info)
995 IF (info /= 0 .OR. result /= 0 .OR. abs(measured_gap - 1.0_dp) > 1.e-12_dp) THEN
996 error stop '3D atomic reference orientation'
997 END IF
998 END DO
999 END SUBROUTINE determinant_controls
1000
1001! **************************************************************************************************
1002!> \brief Fixed tolerances and localizer scale for the small Wilson-Dirac control.
1003!> \param hh Physical Hamiltonian
1004!> \param ss Physical overlap
1005!> \param rr Three Cartesian moments
1006!> \param point Query position
1007!> \param energy Query energy
1008! **************************************************************************************************
1009 SUBROUTINE evaluate(hh, ss, rr, point, energy)
1010 COMPLEX(dp), INTENT(IN) :: hh(:, :), ss(:, :), rr(:, :, :)
1011 REAL(dp), INTENT(IN) :: point(3), energy
1012
1013 CALL localizer_3d_evaluate(hh, ss, rr, point, energy, 0.5_dp, 1.e-10_dp, 1.e-10_dp, 1.e-10_dp, &
1014 index, gap, ratio, status)
1015 END SUBROUTINE evaluate
1016! **************************************************************************************************
1017!> \brief Open 3x3x3 Wilson-Dirac cube in fixed physical up/down ordering.
1018!> \param m Onsite mass; periodic Bloch mass is m+cos(kx)+cos(ky)+cos(kz)
1019! **************************************************************************************************
1020 SUBROUTINE dirac_model(m)
1021 REAL(dp), INTENT(IN) :: m
1022
1023 INTEGER :: col, direction, neighbor, orbital, &
1024 other, pos(3), row, site, spin, x, y, z
1025
1026 h = 0.0_dp
1027 s = 0.0_dp
1028 r = 0.0_dp
1029 DO z = 0, side - 1
1030 DO y = 0, side - 1
1031 DO x = 0, side - 1
1032 pos = [x, y, z]
1033 site = 1 + x + side*y + side**2*z
1034 DO spin = 1, 2
1035 DO orbital = 1, 2
1036 row = site + sites*(orbital - 1) + 2*sites*(spin - 1)
1037 s(row, row) = 1.0_dp
1038 h(row, row) = m*real(3 - 2*orbital, dp)
1039 r(row, row, :) = real(pos - 1, dp)
1040 DO direction = 1, 3
1041 IF (pos(direction) == side - 1) cycle
1042 neighbor = site + side**(direction - 1)
1043 col = neighbor + sites*(orbital - 1) + 2*sites*(spin - 1)
1044 h(row, col) = 0.5_dp*real(3 - 2*orbital, dp)
1045 h(col, row) = conjg(h(row, col))
1046 DO other = 1, 2
1047 col = neighbor + sites*(2 - orbital) + 2*sites*(other - 1)
1048 h(row, col) = -0.5_dp*ii*pauli(spin, other, direction)
1049 h(col, row) = conjg(h(row, col))
1050 END DO
1051 END DO
1052 END DO
1053 END DO
1054 END DO
1055 END DO
1056 END DO
1057 END SUBROUTINE dirac_model
1058 END SUBROUTINE check_localizer_3d
1059
1060! **************************************************************************************************
1061!> \brief Complex metric-covariant matrix sign against exact spectral constructions.
1062!> \param para_env Parallel environment initialized by the unit-test driver.
1063!> \author Thomas D. Kuehne, tkuehne@cp2k.org
1064! **************************************************************************************************
1065 SUBROUTINE check_localizer_flattening(para_env)
1066 USE cp_dbcsr_api, ONLY: &
1067 dbcsr_create, dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, &
1068 dbcsr_finalize, dbcsr_finalize_lib, dbcsr_init_lib, dbcsr_put_block, dbcsr_type_no_symmetry
1069 USE fft_tools, ONLY: finalize_fft, &
1070 init_fft
1071 USE ieee_arithmetic, ONLY: ieee_quiet_nan, &
1072 ieee_value
1073 USE kinds, ONLY: dp
1074 USE mathconstants, ONLY: twopi
1075 USE message_passing, ONLY: mp_para_env_type
1076 USE spectral_localizer_bloch, ONLY: localizer_flatten_bloch
1077 USE spectral_localizer_dbcsr, ONLY: localizer_dbcsr_dense, &
1078 localizer_dbcsr_release, &
1079 localizer_matrix_type
1080 USE spectral_localizer_flatten, ONLY: localizer_flatten_checked, &
1081 localizer_flatten_dbcsr
1082#include "./base/base_uses.f90"
1083
1084 IMPLICIT NONE
1085 INTEGER, PARAMETER :: n = 12, nb = 6, bs = 2
1086 TYPE(mp_para_env_type), POINTER :: para_env
1087 TYPE(dbcsr_distribution_type) :: distribution
1088 TYPE(localizer_matrix_type) :: h, s, flat
1089 INTEGER, POINTER :: process_grid(:, :), row_dist(:), col_dist(:), sizes(:)
1090 COMPLEX(dp) :: transform(n, n), u(n, n), rotation(n, n), scaled(n, n), hdense(n, n), sdense(n, n), exact(n, n)
1091 COMPLEX(dp), ALLOCATABLE :: result(:, :)
1092 REAL(dp) :: energies(n), residuals(4), difference, maximum_error, energy, scale, gap(2), metric_ratio, exact_gap
1093 INTEGER :: i, j, mode, nr, nc, io, status
1094 CHARACTER(LEN=32) :: argument
1095
1096 io = cp_logger_get_default_io_unit()
1097 CALL dbcsr_init_lib(para_env%get_handle(), io)
1098 CALL init_fft('FFTSG', alltoall=.false., pool_limit=4, wisdom_file='', plan_style=1)
1099 nr = para_env%num_pe
1100 nc = 1
1101 DO i = 1, command_argument_count()
1102 CALL get_command_argument(i, argument)
1103 IF (trim(argument) == '--column-grid') THEN
1104 nr = 1
1105 nc = para_env%num_pe
1106 END IF
1107 END DO
1108 ALLOCATE (process_grid(0:nr - 1, 0:nc - 1), row_dist(nb), col_dist(nb), sizes(nb))
1109 process_grid = reshape([(i, i=0, para_env%num_pe - 1)], [nr, nc])
1110 row_dist = [(modulo(i - 1, nr), i=1, nb)]
1111 col_dist = [(modulo(i - 1, nc), i=1, nb)]
1112 sizes = bs
1113 CALL dbcsr_distribution_new(distribution, group=para_env%get_handle(), pgrid=process_grid, &
1114 row_dist=row_dist, col_dist=col_dist)
1115 energy = 0.002_dp
1116 scale = 1.7_dp
1117 maximum_error = 0.0_dp
1118 DO mode = 0, 5
1119 u = 0.0_dp
1120 transform = 0.0_dp
1121 DO i = 1, n
1122 energies(i) = 0.01_dp*(i - n/2) - 0.005_dp
1123 transform(i, i) = 1.0_dp
1124 DO j = 1, n
1125 u(i, j) = exp(cmplx(0.0_dp, twopi*real((i - 1)*(j - 1), dp)/n, dp))/sqrt(real(n, dp))
1126 END DO
1127 IF (mode == 0) cycle
1128 transform(i, i) = 0.5_dp + 0.1_dp*i
1129 IF (i < n) transform(i, i + 1) = cmplx(0.13_dp, -0.08_dp, dp)
1130 END DO
1131 IF (mode == 2) transform(n, n) = 0.1_dp
1132 IF (mode == 4) energies = energies + 0.5_dp
1133 IF (mode == 5) energies = energies - 0.5_dp
1134 IF (mode == 3) THEN
1135 u = 0.0_dp
1136 rotation = 0.0_dp
1137 transform(:n/2, n/2 + 1:) = 0.0_dp
1138 transform(n/2 + 1:, :n/2) = 0.0_dp
1139 transform(n/2 + 1:, n/2 + 1:) = conjg(transform(:n/2, :n/2))
1140 DO i = 1, n/2
1141 energies(i) = 0.01_dp*(i - n/4) - 0.005_dp
1142 energies(n/2 + i) = energies(i)
1143 rotation(i, i) = cos(0.37_dp)
1144 rotation(n/2 + i, n/2 + i) = cos(0.37_dp)
1145 rotation(i, n/2 + i) = sin(0.37_dp)
1146 rotation(n/2 + i, i) = -sin(0.37_dp)
1147 DO j = 1, n/2
1148 u(i, j) = exp(cmplx(0.0_dp, twopi*real((i - 1)*(j - 1), dp)/(n/2), dp))/sqrt(real(n/2, dp))
1149 END DO
1150 END DO
1151 u(n/2 + 1:, n/2 + 1:) = conjg(u(:n/2, :n/2))
1152 u = matmul(rotation, u)
1153 END IF
1154 scaled = u*spread(energies, 1, n)
1155 hdense = matmul(conjg(transpose(transform)), matmul(matmul(scaled, conjg(transpose(u))), transform))
1156 sdense = matmul(conjg(transpose(transform)), transform)
1157 scaled = u*spread(sign(scale, energies - energy), 1, n)
1158 exact = matmul(conjg(transpose(transform)), matmul(matmul(scaled, conjg(transpose(u))), transform))
1159 CALL make_operator(hdense, h)
1160 CALL make_operator(sdense, s)
1161 CALL localizer_flatten_dbcsr(h, s, energy, scale, 1.e-9_dp, para_env, flat, residuals, status)
1162 IF (status /= 0) THEN
1163 IF (io > 0) WRITE (io, *) 'Flatten failure: ', mode, status, residuals
1164 cpabort('Distributed spectral flattening failed')
1165 END IF
1166 CALL localizer_dbcsr_dense(flat, para_env, result)
1167 difference = maxval(abs(result - exact))
1168 maximum_error = max(maximum_error, difference)
1169 IF (difference > 1.e-9_dp) THEN
1170 cpabort('Covariant sign differs from exact eigensystem')
1171 END IF
1172 IF (mode == 3) THEN
1173 IF (maxval(abs(result(:n/2, :n/2) - conjg(result(n/2 + 1:, n/2 + 1:)))) > 1.e-9_dp) THEN
1174 cpabort('Flattening broke physical TR diagonal blocks')
1175 END IF
1176 IF (maxval(abs(result(:n/2, n/2 + 1:) + conjg(result(n/2 + 1:, :n/2)))) > 1.e-9_dp) THEN
1177 cpabort('Flattening broke physical TR spin mixing')
1178 END IF
1179 END IF
1180 DEALLOCATE (result)
1181 CALL localizer_dbcsr_dense(h, para_env, result)
1182 IF (maxval(abs(result - hdense)) > 1.e-14_dp) THEN
1183 cpabort('Flattening modified input Hamiltonian')
1184 END IF
1185 DEALLOCATE (result)
1186 CALL localizer_dbcsr_dense(s, para_env, result)
1187 IF (maxval(abs(result - sdense)) > 1.e-14_dp) THEN
1188 cpabort('Flattening modified input metric')
1189 END IF
1190 DEALLOCATE (result)
1191 CALL localizer_dbcsr_release(flat)
1192 CALL localizer_flatten_checked(h, s, energy, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1193 para_env, flat, gap, metric_ratio, residuals, status)
1194 IF (status /= 0) cpabort('Checked flattening failed')
1195 exact_gap = minval(abs(energies - energy))
1196 IF (gap(1) > exact_gap + 1.e-10_dp .OR. gap(2) < exact_gap - 1.e-10_dp) THEN
1197 cpabort('Input electronic gap differs from reference')
1198 END IF
1199 IF (metric_ratio <= 1.e-10_dp) cpabort('Invalid positive metric bound')
1200 CALL localizer_dbcsr_dense(flat, para_env, result)
1201 difference = maxval(abs(result - exact))
1202 maximum_error = max(maximum_error, difference)
1203 IF (difference > 1.e-9_dp) THEN
1204 cpabort('Checked flattening changed the exact result')
1205 END IF
1206 DEALLOCATE (result)
1207 CALL localizer_dbcsr_release(flat)
1208 CALL localizer_dbcsr_release(h)
1209 CALL localizer_dbcsr_release(s)
1210 END DO
1211
1212 hdense = 0.0_dp
1213 sdense = 0.0_dp
1214 DO i = 1, n
1215 sdense(i, i) = 1.0_dp
1216 hdense(i, i) = 1.0_dp
1217 END DO
1218 hdense(1, 1) = 0.0_dp
1219 CALL make_operator(hdense, h)
1220 CALL make_operator(sdense, s)
1221 CALL localizer_flatten_dbcsr(h, s, 0.0_dp, scale, 1.e-9_dp, para_env, flat, residuals, status)
1222 IF (status /= 1) cpabort('Singular spectral sign was accepted')
1223 CALL localizer_flatten_checked(h, s, 0.0_dp, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1224 para_env, flat, gap, metric_ratio, residuals, status)
1225 IF (status /= 1) cpabort('Checked flattening accepted a zero spectral mode')
1226 CALL localizer_flatten_dbcsr(h, s, energy, -scale, 1.e-9_dp, para_env, flat, residuals, status)
1227 IF (status /= -1) cpabort('Negative flattening scale was accepted')
1228 CALL localizer_flatten_dbcsr(h, s, ieee_value(energy, ieee_quiet_nan), scale, 1.e-9_dp, &
1229 para_env, flat, residuals, status)
1230 IF (status /= -1) cpabort('Nonfinite energy was accepted')
1231 CALL localizer_dbcsr_release(h)
1232 hdense = 0.0_dp
1233 CALL make_operator(hdense, h)
1234 CALL localizer_flatten_dbcsr(h, s, 0.0_dp, scale, 1.e-9_dp, para_env, flat, residuals, status)
1235 IF (status /= 1) cpabort('Zero Hamiltonian sign was accepted')
1236 CALL localizer_dbcsr_release(h)
1237 hdense(1, 2) = 1.0_dp
1238 CALL make_operator(hdense, h)
1239 CALL localizer_flatten_dbcsr(h, s, energy, scale, 1.e-9_dp, para_env, flat, residuals, status)
1240 IF (status /= -1) cpabort('Non-Hermitian Hamiltonian was accepted')
1241 CALL localizer_dbcsr_release(h)
1242 CALL localizer_dbcsr_release(s)
1243 hdense = sdense
1244 hdense(1, 1) = 1.e-11_dp
1245 CALL make_operator(hdense, h)
1246 CALL make_operator(sdense, s)
1247 CALL localizer_flatten_checked(h, s, 0.0_dp, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1248 para_env, flat, gap, metric_ratio, residuals, status)
1249 IF (status /= 1) cpabort('Unresolved near-zero spectral mode was accepted')
1250 CALL localizer_dbcsr_release(s)
1251 sdense(1, 1) = -1.0_dp
1252 CALL make_operator(sdense, s)
1253 CALL localizer_flatten_checked(h, s, 0.0_dp, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1254 para_env, flat, gap, metric_ratio, residuals, status)
1255 IF (status /= -2) cpabort('Indefinite metric was accepted')
1256 CALL localizer_dbcsr_release(s)
1257 sdense(1, 1) = 1.e-12_dp
1258 CALL make_operator(sdense, s)
1259 CALL localizer_flatten_checked(h, s, 0.0_dp, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1260 para_env, flat, gap, metric_ratio, residuals, status)
1261 IF (status /= -2) cpabort('Ill-conditioned metric was accepted')
1262 CALL localizer_dbcsr_release(h)
1263 CALL localizer_dbcsr_release(s)
1264 CALL check_bloch()
1265 IF (io > 0) WRITE (io, '(A,ES14.5)') 'Spectral flattening exact-matrix maximum error: ', maximum_error
1266 CALL dbcsr_distribution_release(distribution)
1267 DEALLOCATE (process_grid, row_dist, col_dist, sizes)
1268 CALL dbcsr_finalize_lib()
1269 CALL finalize_fft(para_env, '')
1270
1271 CONTAINS
1272
1273! **************************************************************************************************
1274!> \brief Complex nonorthogonal full-band Fourier references in spin/cell/atom order.
1275! **************************************************************************************************
1276 SUBROUTINE check_bloch()
1277 INTEGER :: band, grid(3), k, kr(3), m, ncell, &
1278 nscalar, nsmall, nspin, p, pnext, &
1279 r(3), rnext, row, t, tnext
1280 REAL(dp) :: bloch_error, bloch_gap, norm_error, &
1281 phase, ref_gap
1282
1283 bloch_error = 0.0_dp
1284 DO m = 1, 7
1285 sizes = bs
1286 grid = [3, 1, 1]
1287 nspin = 2
1288 IF (m == 2) grid = [1, 3, 1]
1289 IF (m == 3) grid = [1, 1, 3]
1290 IF (m == 4) THEN
1291 grid = [2, 3, 1]
1292 nspin = 1
1293 END IF
1294 IF (m == 5) THEN
1295 nspin = 1
1296 sizes = [1, 3, 1, 3, 1, 3]
1297 END IF
1298 IF (m == 6) grid = 1
1299 ncell = product(grid)
1300 nscalar = n/(nspin*ncell)
1301 nsmall = nscalar*nspin
1302 u = 0.0_dp
1303 transform = 0.0_dp
1304 DO row = 1, n
1305 t = modulo((row - 1)/nscalar, ncell)
1306 p = ((row - 1)/(ncell*nscalar))*nscalar + modulo(row - 1, nscalar)
1307 r = [modulo(t, grid(1)), modulo(t/grid(1), grid(2)), t/(grid(1)*grid(2))]
1308 transform(row, row) = 1.0_dp + 0.03_dp*p
1309 pnext = modulo(p + 1, nsmall)
1310 rnext = (pnext/nscalar)*ncell*nscalar + t*nscalar + modulo(pnext, nscalar) + 1
1311 transform(row, rnext) = transform(row, rnext) + cmplx(0.07_dp, -0.04_dp, dp)
1312 tnext = modulo(t + 1, ncell)
1313 rnext = (p/nscalar)*ncell*nscalar + tnext*nscalar + modulo(p, nscalar) + 1
1314 ! A 1D cyclic displacement is a lattice translation for these first-axis-fast grids.
1315 IF (m /= 4) transform(row, rnext) = transform(row, rnext) + cmplx(0.05_dp, 0.02_dp, dp)
1316 DO k = 0, ncell - 1
1317 kr = [modulo(k, grid(1)), modulo(k/grid(1), grid(2)), k/(grid(1)*grid(2))]
1318 DO band = 0, nsmall - 1
1319 phase = twopi*(sum(real(kr*r, dp)/real(grid, dp)) + real(p*band, dp)/real(nsmall, dp))
1320 u(row, k*nsmall + band + 1) = exp(cmplx(0.0_dp, phase, dp))/sqrt(real(n, dp))
1321 energies(k*nsmall + band + 1) = 0.2_dp + 0.03_dp*band + 0.01_dp*cos(real(k, dp))
1322 IF (m == 7) energies(k*nsmall + band + 1) = 0.25_dp + 0.01_dp*cos(real(k, dp))
1323 IF (band < nsmall/2) energies(k*nsmall + band + 1) = -energies(k*nsmall + band + 1)
1324 END DO
1325 END DO
1326 END DO
1327 scaled = u*spread(energies, 1, n)
1328 hdense = matmul(conjg(transpose(transform)), matmul(matmul(scaled, conjg(transpose(u))), transform))
1329 sdense = matmul(conjg(transpose(transform)), transform)
1330 scaled = u*spread(sign(scale, energies - energy), 1, n)
1331 exact = matmul(conjg(transpose(transform)), matmul(matmul(scaled, conjg(transpose(u))), transform))
1332 CALL make_operator(hdense, h)
1333 CALL make_operator(sdense, s)
1334 CALL localizer_flatten_bloch(h, s, grid, nspin, energy, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1335 para_env, flat, bloch_gap, metric_ratio, residuals, status)
1336 IF (status /= 0) cpabort('Bloch flattening failed')
1337 CALL localizer_dbcsr_dense(flat, para_env, result)
1338 norm_error = maxval(abs(result - exact))
1339 bloch_error = max(bloch_error, norm_error)
1340 maximum_error = max(maximum_error, norm_error)
1341 ref_gap = minval(abs(energies - energy))
1342 IF (norm_error > 1.e-11_dp .OR. abs(bloch_gap - ref_gap) > 1.e-12_dp) THEN
1343 cpabort('Bloch flattening differs from exact full-torus reference')
1344 END IF
1345 CALL localizer_dbcsr_release(flat)
1346 CALL localizer_flatten_checked(h, s, energy, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1347 para_env, flat, gap, metric_ratio, residuals, status)
1348 IF (status /= 0) cpabort('Full-torus reference flattening failed')
1349 CALL localizer_dbcsr_dense(flat, para_env, result)
1350 IF (maxval(abs(result - exact)) > 1.e-10_dp) THEN
1351 cpabort('Full-torus/Bloch sign mismatch')
1352 END IF
1353 CALL localizer_dbcsr_release(flat)
1354 CALL localizer_dbcsr_release(h)
1355 IF (m /= 6) THEN
1356 hdense(1, 1) = hdense(1, 1) + 0.01_dp
1357 CALL make_operator(hdense, h)
1358 CALL localizer_flatten_bloch(h, s, grid, nspin, energy, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1359 para_env, flat, bloch_gap, metric_ratio, residuals, status)
1360 IF (status /= -7) cpabort('Broken translations accepted')
1361 CALL localizer_dbcsr_release(h)
1362 END IF
1363 CALL make_operator(energy*sdense, h)
1364 CALL localizer_flatten_bloch(h, s, grid, nspin, energy, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1365 para_env, flat, bloch_gap, metric_ratio, residuals, status)
1366 IF (status /= 1) cpabort('Closed Bloch gap accepted')
1367 CALL localizer_dbcsr_release(s)
1368 CALL make_operator(-sdense, s)
1369 CALL localizer_flatten_bloch(h, s, grid, nspin, energy, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1370 para_env, flat, bloch_gap, metric_ratio, residuals, status)
1371 IF (status /= -2) cpabort('Bloch flattening accepted an indefinite metric')
1372 CALL localizer_dbcsr_release(h)
1373 CALL localizer_dbcsr_release(s)
1374 END DO
1375 sizes = bs
1376 ! Each Bloch metric is well conditioned, but their union has an unacceptable ratio.
1377 hdense = 0.0_dp
1378 sdense = 0.0_dp
1379 DO row = 1, n
1380 hdense(row, row) = 1.0_dp
1381 sdense(row, row) = 0.5_dp*(1.0_dp + 1.e-12_dp)
1382 rnext = 1 + modulo(row - 1 + n/2, n)
1383 sdense(row, rnext) = 0.5_dp*(1.0_dp - 1.e-12_dp)
1384 END DO
1385 CALL make_operator(hdense, h)
1386 CALL make_operator(sdense, s)
1387 CALL localizer_flatten_bloch(h, s, [2, 1, 1], 1, 0.0_dp, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1388 para_env, flat, bloch_gap, metric_ratio, residuals, status)
1389 IF (status /= -2) cpabort('Global Bloch metric bound failed')
1390 CALL localizer_flatten_bloch(h, s, [0, 1, 1], 1, 0.0_dp, scale, 1.e-9_dp, 1.e-10_dp, 1.e-8_dp, &
1391 para_env, flat, bloch_gap, metric_ratio, residuals, status)
1392 IF (status /= -1) cpabort('Invalid Bloch grid accepted')
1393 CALL localizer_dbcsr_release(h)
1394 CALL localizer_dbcsr_release(s)
1395 IF (io > 0) WRITE (io, '(A,ES14.5)') 'Seven Bloch flattening pencils; maximum error: ', bloch_error
1396 END SUBROUTINE check_bloch
1397
1398! **************************************************************************************************
1399!> \brief Insert only locally owned dense-reference blocks.
1400!> \param dense Reference matrix
1401!> \param matrix New complex DBCSR operator
1402! **************************************************************************************************
1403 SUBROUTINE make_operator(dense, matrix)
1404 COMPLEX(dp), INTENT(IN) :: dense(:, :)
1405 TYPE(localizer_matrix_type), INTENT(OUT) :: matrix
1406
1407 INTEGER :: col, first, last, left, right, row
1408
1409 CALL dbcsr_create(matrix%real_part, 'FLATTEN REAL', distribution, dbcsr_type_no_symmetry, sizes, sizes)
1410 CALL dbcsr_create(matrix%imaginary_part, template=matrix%real_part)
1411 DO row = 1, nb
1412 first = sum(sizes(:row - 1)) + 1
1413 last = sum(sizes(:row))
1414 DO col = 1, nb
1415 IF (process_grid(row_dist(row), col_dist(col)) /= para_env%mepos) cycle
1416 left = sum(sizes(:col - 1)) + 1
1417 right = sum(sizes(:col))
1418 CALL dbcsr_put_block(matrix%real_part, row, col, real(dense(first:last, left:right), dp))
1419 CALL dbcsr_put_block(matrix%imaginary_part, row, col, aimag(dense(first:last, left:right)))
1420 END DO
1421 END DO
1422 CALL dbcsr_finalize(matrix%real_part)
1423 CALL dbcsr_finalize(matrix%imaginary_part)
1424 END SUBROUTINE make_operator
1425
1426 END SUBROUTINE check_localizer_flattening
1427
1428END MODULE spectral_localizer_test_cases
1429
1430PROGRAM spectral_localizer_unittest
1431 USE cp_log_handling, ONLY: cp_add_default_logger, &
1432 cp_logger_create, &
1433 cp_logger_release, &
1434 cp_logger_type, &
1435 cp_rm_default_logger
1436 USE machine, ONLY: default_output_unit
1437 USE message_passing, ONLY: mp_para_env_type, mp_world_finalize, mp_world_init
1438 USE spectral_localizer_test_cases, ONLY: check_localizer_2d, check_localizer_3d, check_localizer_flattening
1439 IMPLICIT NONE
1440 TYPE(mp_para_env_type), POINTER :: para_env
1441 TYPE(cp_logger_type), POINTER :: logger
1442 INTEGER :: io
1443
1444 ALLOCATE (para_env)
1445 CALL mp_world_init(para_env)
1446 io = -1
1447 IF (para_env%is_source()) io = default_output_unit
1448 CALL cp_logger_create(logger, para_env=para_env, default_global_unit_nr=io, close_global_unit_on_dealloc=.false.)
1449 CALL cp_add_default_logger(logger)
1450 CALL check_localizer_2d()
1451 CALL check_localizer_3d()
1452 CALL check_localizer_flattening(para_env)
1453 CALL cp_rm_default_logger()
1454 CALL cp_logger_release(logger)
1455 DEALLOCATE (para_env)
1456 CALL mp_world_finalize()
1457END PROGRAM spectral_localizer_unittest
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:126
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
spectral_localizer::localizer_inertia
Definition spectral_localizer.F:25
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:643
cp_dbcsr_api::dbcsr_distribution_new
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
Definition cp_dbcsr_api.F:616
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1125
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_rm_default_logger
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
Definition cp_log_handling.F:213
cp_log_handling::cp_logger_release
subroutine, public cp_logger_release(logger)
releases this logger
Definition cp_log_handling.F:429
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_logger_create
subroutine, public cp_logger_create(logger, para_env, print_level, default_global_unit_nr, default_local_unit_nr, global_filename, local_filename, close_global_unit_on_dealloc, iter_info, close_local_unit_on_dealloc, suffix, template_logger)
initializes a logger
Definition cp_log_handling.F:285
cp_log_handling::cp_add_default_logger
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
Definition cp_log_handling.F:188
fft_tools
Definition fft_tools.F:29
fft_tools::finalize_fft
subroutine, public finalize_fft(para_env, wisdom_file, export_wisdom)
does whatever is needed to finalize the current fft setup
Definition fft_tools.F:209
fft_tools::init_fft
subroutine, public init_fft(fftlib, alltoall, pool_limit, wisdom_file, plan_style)
...
Definition fft_tools.F:181
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::default_output_unit
integer, parameter, public default_output_unit
Definition machine.F:48
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
message_passing::mp_world_init
subroutine, public mp_world_init(mp_comm)
initializes the system default communicator
Definition message_passing.F:1389
message_passing::mp_world_finalize
subroutine, public mp_world_finalize()
Finalize the system default communicator and MPI when CP2K owns MPI.
Definition message_passing.F:1481
spectral_localizer_bloch
Exact spectral flattening of a translation-invariant AO torus in primitive Bloch blocks.
Definition spectral_localizer_bloch.F:12
spectral_localizer_bloch::localizer_flatten_bloch
subroutine, public localizer_flatten_bloch(h, s, grid, nspin, energy, scale, tolerance, metric_tolerance, gap_tolerance, para_env, flattened, gap, metric_ratio, residuals, status)
Flatten all Bloch bands, then return to the original distributed real-cell/spin AO basis.
Definition spectral_localizer_bloch.F:56
spectral_localizer_dbcsr
Distributed block-sparse finite localizers, retaining the covariant AO metric.
Definition spectral_localizer_dbcsr.F:12
spectral_localizer_dbcsr::localizer_dbcsr_dense
subroutine, public localizer_dbcsr_dense(matrix, para_env, dense)
Collect the size-limited dense oracle; sparse solvers must not call this routine.
Definition spectral_localizer_dbcsr.F:465
spectral_localizer_dbcsr::localizer_dbcsr_release
subroutine, public localizer_dbcsr_release(matrix)
Release both real components of a complex sparse operator.
Definition spectral_localizer_dbcsr.F:594
spectral_localizer_flatten
Metric-covariant spectral flattening with existing distributed matrix iterations.
Definition spectral_localizer_flatten.F:12
spectral_localizer_flatten::localizer_flatten_checked
subroutine, public localizer_flatten_checked(h, s, energy, scale, tolerance, metric_tolerance, gap_tolerance, para_env, flattened, gap, metric_ratio, residuals, status)
Establish a positive AO metric and electronic gap before iterating the matrix sign.
Definition spectral_localizer_flatten.F:60
spectral_localizer_flatten::localizer_flatten_dbcsr
subroutine, public localizer_flatten_dbcsr(h, s, energy, scale, tolerance, para_env, flattened, residuals, status)
Form scale*S*sign(S^-1*(H-energy*S)) without dense AO matrices or S^-1/2.
Definition spectral_localizer_flatten.F:126
spectral_localizer_test_cases
Definition spectral_localizer_unittest.F:8
spectral_localizer_test_cases::check_localizer_3d
subroutine, public check_localizer_3d()
Independent chiral singular values, strong/weak/trivial Wilson-Dirac models and AO covariance.
Definition spectral_localizer_unittest.F:784
spectral_localizer_test_cases::check_localizer_2d
subroutine, public check_localizer_2d()
Independent finite-model, metric, orientation and Pfaffian tests.
Definition spectral_localizer_unittest.F:21
spectral_localizer_test_cases::check_localizer_flattening
subroutine, public check_localizer_flattening(para_env)
Complex metric-covariant matrix sign against exact spectral constructions.
Definition spectral_localizer_unittest.F:1066
spectral_localizer
Dense reference spectral localizers in a nonorthogonal finite basis.
Definition spectral_localizer.F:12
spectral_localizer::localizer_3d_evaluate
subroutine, public localizer_3d_evaluate(h, s, r, point, energy, kappa, matrix_tol, metric_tol, gap_tol, invariant, gap, metric_ratio, status)
Finite three-dimensional class AII localizer, using all three Cartesian moments.
Definition spectral_localizer.F:49
spectral_localizer::localizer_aii_skew
subroutine, public localizer_aii_skew(localizer, skew, tolerance, residual, status)
Convert the class AII localizer to a real skew matrix using Q=(I-i C)/sqrt(2).
Definition spectral_localizer.F:732
spectral_localizer::localizer_3d_matrix_evaluate
subroutine, public localizer_3d_matrix_evaluate(block, s, matrix_tol, metric_tol, gap_tol, invariant, gap, metric_ratio, status)
Metric gap and real determinant sign of a three-dimensional AII chiral block.
Definition spectral_localizer.F:108
spectral_localizer::localizer_periodic_evaluate
subroutine, public localizer_periodic_evaluate(h, s, cosine, sine, origin, energy, eta, aii, matrix_tol, metric_tol, gap_tol, invariant, gap, metric_ratio, status)
Two-dimensional periodic localizer from covariant trigonometric position operators.
Definition spectral_localizer.F:292
spectral_localizer::localizer_pfaffian_sign
subroutine, public localizer_pfaffian_sign(matrix, pf_sign, status)
Dense pivoted skew elimination, accumulating only the sign of the Pfaffian.
Definition spectral_localizer.F:799
spectral_localizer::localizer_evaluate
subroutine, public localizer_evaluate(h, s, r, point, energy, kappa, aii, matrix_tol, metric_tol, gap_tol, invariant, gap, metric_ratio, status)
Evaluate the 2D class A or AII localizer without projecting onto occupied states.
Definition spectral_localizer.F:364
aii_transform_tests
subroutine aii_transform_tests()
Signed-permutation AII conversion against the explicit dense unitary transform.
Definition spectral_localizer_unittest.F:214
check_bloch
subroutine check_bloch()
Complex nonorthogonal full-band Fourier references in spin/cell/atom order.
Definition spectral_localizer_unittest.F:1277
determinant_controls
subroutine determinant_controls()
Chiral-block determinant controls independent of the Wilson-Dirac model.
Definition spectral_localizer_unittest.F:936
spectral_localizer_unittest
program spectral_localizer_unittest
Definition spectral_localizer_unittest.F:1430
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
orbital
Orbital angular momentum.
Definition grid_common.h:134
spectral_localizer_dbcsr::localizer_matrix_type
Definition spectral_localizer_dbcsr.F:34