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 : MODULE 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
14 : CONTAINS
15 :
16 : ! **************************************************************************************************
17 : !> \brief Independent finite-model, metric, orientation and Pfaffian tests.
18 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
19 : ! **************************************************************************************************
20 2 : 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 2 : COMPLEX(KIND=dp), ALLOCATABLE :: h(:, :), s(:, :), r(:, :, :), transform(:, :), &
36 2 : 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 2 : 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 2 : CALL inertia_tests()
45 2 : CALL real_inertia_tests()
46 2 : CALL gap_tests()
47 2 : CALL pfaffian_orientation_tests()
48 2 : CALL aii_transform_tests()
49 2 : CALL aii_roundoff_tests()
50 2 : CALL periodic_tests()
51 : ! Independent Pfaffian polynomial for size four, including odd permutations.
52 2 : skew = 0.0_dp
53 2 : skew(1, 2) = 0.3_dp
54 2 : skew(1, 3) = 2.0_dp
55 2 : skew(1, 4) = -1.0_dp
56 2 : skew(2, 3) = 0.7_dp
57 2 : skew(2, 4) = 3.0_dp
58 2 : skew(3, 4) = 0.8_dp
59 82 : skew = skew - TRANSPOSE(skew)
60 2 : exact = skew(1, 2)*skew(3, 4) - skew(1, 3)*skew(2, 4) + skew(1, 4)*skew(2, 3)
61 2 : CALL localizer_pfaffian_sign(skew, pf, status)
62 2 : IF (status /= 0 .OR. pf /= INT(SIGN(1.0_dp, exact))) ERROR STOP 'Pfaffian polynomial failed'
63 2 : permutation = [2, 1, 3, 4]
64 42 : CALL localizer_pfaffian_sign(skew(permutation, permutation), pf, status)
65 2 : IF (status /= 0 .OR. pf /= -INT(SIGN(1.0_dp, exact))) ERROR STOP 'Pfaffian orientation failed'
66 2 : skew = 0.0_dp
67 2 : CALL localizer_pfaffian_sign(skew, pf, status)
68 2 : IF (status == 0 .OR. pf /= 0) ERROR STOP 'Singular Pfaffian accepted'
69 :
70 2 : ALLOCATE (h(1, 1), s(1, 1), r(1, 1, 2))
71 6 : h = 4.0_dp
72 6 : s = 2.0_dp
73 6 : r(1, 1, :) = [6.0_dp, 8.0_dp]
74 2 : CALL evaluate(h, s, r, [0.0_dp, 0.0_dp], 1.0_dp, 0.2_dp, .FALSE., index, gap, status)
75 2 : IF (status /= 0 .OR. index /= 0 .OR. ABS(gap - SQRT(2.0_dp)) > 1.e-12_dp) THEN
76 0 : ERROR STOP 'Analytic generalized localizer failed'
77 : END IF
78 2 : CALL evaluate(h, s, r, [3.0_dp, 4.0_dp], 2.0_dp, 0.2_dp, .FALSE., index, gap, status)
79 2 : IF (status /= 1) ERROR STOP 'Zero gap was assigned an index'
80 2 : DEALLOCATE (h, s, r)
81 :
82 2 : ALLOCATE (h(2, 2), s(2, 2), r(2, 2, 2))
83 2 : h = 0.0_dp
84 2 : s = 0.0_dp
85 2 : r = 0.0_dp
86 2 : s(1, 1) = 1.e-9_dp
87 2 : s(2, 2) = 1.0_dp
88 2 : h(1, 1) = 1.e-16_dp
89 2 : h(2, 2) = 1.0_dp
90 2 : CALL evaluate(h, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.2_dp, .FALSE., index, gap, status)
91 2 : IF (status /= 1) ERROR STOP 'Metric-amplified roundoff ignored'
92 2 : DEALLOCATE (h, s, r)
93 :
94 2 : CALL qwz_model(7, -1.0_dp, h, s, r)
95 2 : n = SIZE(h, 1)
96 2 : CALL evaluate(h, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .FALSE., reference, gap_reference, status)
97 2 : IF (status /= 0 .OR. reference /= -1 .OR. ABS(gap_reference - 0.740327333453389_dp) > 1.e-11_dp) THEN
98 0 : ERROR STOP 'Nontrivial open QWZ model failed'
99 : END IF
100 38814 : CALL evaluate(h, s, r(:, :, [2, 1]), [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .FALSE., index, gap, status)
101 2 : IF (status /= 0 .OR. index /= -reference .OR. ABS(gap - gap_reference) > 1.e-11_dp) THEN
102 0 : ERROR STOP 'Oriented-plane reversal failed'
103 : END IF
104 19406 : CALL evaluate(CONJG(h), s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .FALSE., index, gap, status)
105 2 : IF (status /= 0 .OR. index /= -reference) ERROR STOP 'Time reversed Chern sign failed'
106 2 : CALL evaluate(h, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.15_dp, .FALSE., index, gap, status)
107 2 : IF (status /= 0 .OR. index /= reference) ERROR STOP 'Kappa plateau failed'
108 2 : CALL evaluate(h, s, r, [20.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .FALSE., index, gap, status)
109 2 : IF (status /= 0 .OR. index /= 0) ERROR STOP 'Outside-vacuum limit failed'
110 :
111 24 : ALLOCATE (transform(n, n), hc(n, n), sc(n, n), rc(n, n, 2))
112 2 : transform = 0.0_dp
113 198 : DO i = 1, n
114 196 : 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 198 : IF (i < n) transform(i, i + 1) = CMPLX(0.05_dp, 0.03_dp, dp)
116 : END DO
117 19406 : hc(:, :) = MATMUL(CONJG(TRANSPOSE(transform)), MATMUL(h, transform))
118 19406 : sc(:, :) = MATMUL(CONJG(TRANSPOSE(transform)), MATMUL(s, transform))
119 6 : DO axis = 1, 2
120 38814 : rc(:, :, axis) = MATMUL(CONJG(TRANSPOSE(transform)), MATMUL(r(:, :, axis), transform))
121 : END DO
122 2 : CALL evaluate(hc, sc, rc, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .FALSE., index, gap, status)
123 2 : IF (status /= 0 .OR. index /= reference .OR. ABS(gap - gap_reference) > 1.e-11_dp) THEN
124 0 : ERROR STOP 'Complex nonorthogonal basis invariance failed'
125 : END IF
126 19406 : rc(:, :, 1) = rc(:, :, 1) + 13.0_dp*sc
127 19406 : rc(:, :, 2) = rc(:, :, 2) - 7.0_dp*sc
128 19406 : hc(:, :) = hc + 11.0_dp*sc
129 2 : CALL evaluate(hc, sc, rc, [13.0_dp, -7.0_dp], 11.0_dp, 0.25_dp, .FALSE., index, gap, status)
130 2 : IF (status /= 0 .OR. index /= reference .OR. ABS(gap - gap_reference) > 1.e-11_dp) THEN
131 0 : ERROR STOP 'Origin/energy shift covariance failed'
132 : END IF
133 2 : sc = 0.0_dp
134 198 : DO i = 1, n
135 198 : sc(i, i) = 1.0_dp
136 : END DO
137 2 : sc(1, 1) = 1.e-14_dp
138 2 : CALL evaluate(hc, sc, rc, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .FALSE., index, gap, status)
139 2 : IF (status /= -2) ERROR STOP 'Ill-conditioned metric accepted'
140 2 : sc(1, 1) = -1.0_dp
141 2 : CALL evaluate(hc, sc, rc, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .FALSE., index, gap, status)
142 2 : IF (status /= -2) ERROR STOP 'Indefinite metric accepted'
143 2 : hc(1, 2) = hc(1, 2) + 1.0_dp
144 2 : CALL evaluate(hc, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .FALSE., index, gap, status)
145 2 : IF (status /= -1) ERROR STOP 'Non-Hermitian Hamiltonian accepted'
146 :
147 : ! Two time-reversed Chern copies, with TR-preserving, spin-nonconserving mixing.
148 20 : ALLOCATE (hs(2*n, 2*n), ss(2*n, 2*n), rs(2*n, 2*n, 2))
149 2 : hs = 0.0_dp
150 2 : ss = 0.0_dp
151 2 : rs = 0.0_dp
152 19406 : hs(:n, :n) = h
153 19406 : hs(n + 1:, n + 1:) = CONJG(h)
154 19406 : ss(:n, :n) = s
155 19406 : ss(n + 1:, n + 1:) = s
156 38814 : rs(:n, :n, :) = r
157 38814 : rs(n + 1:, n + 1:, :) = r
158 2 : DO i = 1, n, 2
159 98 : hs(i, n + i + 1) = 0.1_dp
160 98 : hs(i + 1, n + i) = -0.1_dp
161 : END DO
162 38810 : hs(n + 1:, :n) = CONJG(TRANSPOSE(hs(:n, n + 1:)))
163 2 : CALL evaluate(hs, ss, rs, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .TRUE., index, gap, status)
164 2 : IF (status /= 0 .OR. index /= 1) ERROR STOP 'Spin-mixed nontrivial AII model failed'
165 2 : gap_reference = gap
166 : ! Real, nonorthogonal AO change applied to both spin components preserves fixed Theta.
167 19406 : transform = REAL(transform, dp)
168 198 : DO i = 1, n
169 198 : transform(i, i) = 1.0_dp + REAL(MOD(i, 5), dp)/5.0_dp
170 : END DO
171 6 : DO axis = 1, 2
172 14 : DO i = 1, 2
173 : hs((i - 1)*n + 1:i*n, (axis - 1)*n + 1:axis*n) = &
174 77628 : 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 4 : sc(:, :) = MATMUL(TRANSPOSE(transform), MATMUL(s, transform))
178 19406 : ss(:n, :n) = sc
179 19406 : ss(n + 1:, n + 1:) = sc
180 6 : DO axis = 1, 2
181 8 : rc(:, :, axis) = MATMUL(TRANSPOSE(transform), MATMUL(r(:, :, axis), transform))
182 38812 : rs(:n, :n, axis) = rc(:, :, axis)
183 38814 : rs(n + 1:, n + 1:, axis) = rc(:, :, axis)
184 : END DO
185 2 : CALL evaluate(hs, ss, rs, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .TRUE., index, gap, status)
186 2 : IF (status /= 0 .OR. index /= 1 .OR. ABS(gap - gap_reference) > 1.e-11_dp) THEN
187 0 : ERROR STOP 'AII metric covariance failed'
188 : END IF
189 2 : hs(1, 1) = hs(1, 1) + 1.e-12_dp
190 2 : CALL evaluate(hs, ss, rs, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .TRUE., index, gap, status)
191 2 : IF (status /= 0 .OR. index /= 1) ERROR STOP 'Roundoff-level AII residual failed'
192 2 : hs(1, 1) = hs(1, 1) + 1.0_dp
193 2 : CALL evaluate(hs, ss, rs, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .TRUE., index, gap, status)
194 2 : IF (status /= -3) ERROR STOP 'Broken physical time reversal accepted'
195 77226 : hs(:, :) = ss
196 2 : CALL evaluate(hs, ss, rs, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .TRUE., index, gap, status)
197 2 : IF (status /= 0 .OR. index /= 0) ERROR STOP 'Atomic AII calibration failed'
198 :
199 2 : DEALLOCATE (h, s, r)
200 2 : CALL qwz_model(7, -3.0_dp, h, s, r)
201 2 : CALL evaluate(h, s, r, [0.0_dp, 0.0_dp], 0.0_dp, 0.25_dp, .FALSE., index, gap, status)
202 2 : 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 2 : SUBROUTINE aii_transform_tests()
214 2 : COMPLEX(dp), ALLOCATABLE :: matrix(:, :), q(:, :), reference(:, :)
215 : INTEGER :: i, j, m, n, status
216 : REAL(dp) :: error, largest, residual
217 2 : REAL(dp), ALLOCATABLE :: actual(:, :), expected(:, :)
218 :
219 2 : largest = 0.0_dp
220 28 : DO n = 1, 13
221 26 : m = 4*n
222 338 : ALLOCATE (q(m, m), matrix(m, m), reference(m, m), expected(m, m), actual(m, m))
223 26 : q = 0.0_dp
224 754 : DO i = 1, m
225 728 : q(i, i) = 1.0_dp
226 26962 : DO j = 1, m
227 26936 : expected(i, j) = SIN(REAL(i + 2*j, dp)) + COS(REAL(3*i - j, dp))
228 : END DO
229 : END DO
230 53898 : expected(:, :) = expected - TRANSPOSE(expected)
231 208 : DO i = 1, n
232 182 : q(i, 3*n + i) = CMPLX(0.0_dp, -1.0_dp, dp)
233 182 : q(n + i, 2*n + i) = CMPLX(0.0_dp, 1.0_dp, dp)
234 182 : q(2*n + i, n + i) = CMPLX(0.0_dp, 1.0_dp, dp)
235 208 : q(3*n + i, i) = CMPLX(0.0_dp, -1.0_dp, dp)
236 : END DO
237 26962 : q = q/SQRT(2.0_dp)
238 687730 : matrix(:, :) = CMPLX(0.0_dp, -1.0_dp, dp)*MATMUL(q, MATMUL(expected, CONJG(TRANSPOSE(q))))
239 433116 : reference(:, :) = MATMUL(CONJG(TRANSPOSE(q)), MATMUL(matrix, q))
240 26 : CALL localizer_aii_skew(matrix, actual, 1.e-10_dp, residual, status)
241 53924 : error = MAX(MAXVAL(ABS(actual - expected)), MAXVAL(ABS(actual + AIMAG(reference))))
242 26 : IF (status /= 0 .OR. error > 1.e-13_dp) ERROR STOP 'AII sparse-unitary transform failed'
243 26 : largest = MAX(largest, error)
244 754 : DO i = 1, m
245 754 : matrix(i, i) = matrix(i, i) - CMPLX(0.0_dp, 5.e-12_dp, dp)
246 : END DO
247 26 : CALL localizer_aii_skew(matrix, actual, 1.e-10_dp, residual, status)
248 26962 : IF (status /= 0 .OR. MAXVAL(ABS(actual - expected)) > 1.e-13_dp .OR. residual < 4.9e-12_dp) THEN
249 0 : ERROR STOP 'AII skew projection or discarded-error bound failed'
250 : END IF
251 26962 : IF (MAXVAL(ABS(actual + TRANSPOSE(actual))) /= 0.0_dp) ERROR STOP 'AII projection is not exactly skew'
252 754 : DO i = 1, m
253 754 : matrix(i, i) = matrix(i, i) + CMPLX(0.0_dp, 5.e-12_dp, dp)
254 : END DO
255 26 : matrix(1, 1) = matrix(1, 1) + 0.1_dp
256 26 : CALL localizer_aii_skew(matrix, actual, 1.e-10_dp, residual, status)
257 26 : IF (status == 0) ERROR STOP 'Broken AII transform accepted'
258 26 : CALL localizer_aii_skew(matrix(:m - 1, :m - 1), actual, 1.e-10_dp, residual, status)
259 26 : IF (status == 0) ERROR STOP 'Invalid AII dimensions accepted'
260 28 : DEALLOCATE (q, matrix, reference, expected, actual)
261 : END DO
262 2 : IF (io > 0) WRITE (io, '(A,ES14.5)') 'AII signed-permutation maximum error: ', largest
263 4 : END SUBROUTINE aii_transform_tests
264 :
265 : ! **************************************************************************************************
266 : !> \brief Accepted Hermiticity roundoff must not fail the stricter Pfaffian skew check.
267 : ! **************************************************************************************************
268 2 : 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 2 : h0 = 0.0_dp
274 2 : s0 = 0.0_dp
275 2 : r0 = 0.0_dp
276 10 : DO i = 1, 4
277 8 : s0(i, i) = 1.0_dp
278 10 : h0(i, i) = REAL(3*MOD(i + 1, 2) - 1, dp)
279 : END DO
280 2 : h0(1, 2) = 5.e-12_dp
281 2 : h0(3, 4) = 5.e-12_dp
282 2 : CALL evaluate(h0, s0, r0, [0.0_dp, 0.0_dp], 0.0_dp, 0.2_dp, .TRUE., index, gap, status)
283 2 : IF (status /= 0 .OR. index /= 0 .OR. ABS(gap - 1.0_dp) > 1.e-12_dp) ERROR STOP 'AII roundoff rejected'
284 2 : h0(1, 1) = -1.e-10_dp
285 2 : h0(3, 3) = -1.e-10_dp
286 2 : CALL evaluate(h0, s0, r0, [0.0_dp, 0.0_dp], 0.0_dp, 0.2_dp, .TRUE., index, gap, status)
287 2 : IF (status /= 1) ERROR STOP 'AII unresolved gap assigned an index'
288 2 : h0(1, 2) = 1.e-5_dp
289 2 : CALL evaluate(h0, s0, r0, [0.0_dp, 0.0_dp], 0.0_dp, 0.2_dp, .TRUE., index, gap, status)
290 2 : IF (status /= -1) ERROR STOP 'AII non-Hermitian input projected away'
291 2 : END SUBROUTINE aii_roundoff_tests
292 :
293 : ! **************************************************************************************************
294 : !> \brief Torus models with known Chern/Z2 phases, origin periodicity and metric covariance.
295 : ! **************************************************************************************************
296 2 : SUBROUTINE periodic_tests()
297 : INTEGER, PARAMETER :: length = 8
298 :
299 2 : COMPLEX(dp), ALLOCATABLE :: co(:, :, :), ct(:, :, :), h0(:, :), ht(:, :), pos(:, :, :), &
300 2 : s0(:, :), si(:, :, :), sn(:, :, :), spin_c(:, :, :), spin_h(:, :), spin_s(:, :), &
301 2 : 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 2 : two_pi = 8.0_dp*ATAN(1.0_dp)
307 8 : DO phase = 1, 3
308 6 : mass = -1.0_dp
309 6 : expected = -1
310 6 : IF (phase == 2) THEN
311 2 : mass = 1.0_dp
312 2 : expected = 1
313 4 : ELSE IF (phase == 3) THEN
314 2 : mass = -3.0_dp
315 2 : expected = 0
316 : END IF
317 6 : CALL qwz_model(length, mass, h0, s0, pos, torus=.TRUE.)
318 6 : n = SIZE(h0, 1)
319 42 : ALLOCATE (co(n, n, 2), si(n, n, 2))
320 6 : co = 0.0_dp
321 6 : si = 0.0_dp
322 18 : DO j = 1, 2
323 1554 : DO i = 1, n
324 1536 : co(i, i, j) = COS(two_pi*REAL(pos(i, i, j), dp)/length)
325 1548 : 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 6 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g0, cond, info)
330 6 : 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 18 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
333 6 : IF (info /= 0 .OR. idx /= expected .OR. ABS(g - g0) > 1.e-11_dp) THEN
334 0 : 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 18 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
339 6 : IF (info /= 0 .OR. idx /= expected .OR. ABS(g - g0) > 1.e-11_dp) THEN
340 0 : 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 396318 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
345 6 : IF (info /= 0 .OR. idx /= -expected .OR. ABS(g - g0) > 1.e-11_dp) THEN
346 0 : ERROR STOP 'Periodic orientation reversal failed'
347 : END IF
348 :
349 72 : 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 6 : spin_h = 0.0_dp
351 6 : spin_s = 0.0_dp
352 6 : spin_c = 0.0_dp
353 6 : spin_sin = 0.0_dp
354 99078 : spin_h(:n, :n) = h0
355 99078 : spin_h(n + 1:, n + 1:) = CONJG(h0)
356 6 : DO i = 1, n, 2
357 384 : spin_h(i, n + i + 1) = 0.1_dp
358 384 : spin_h(i + 1, n + i) = -0.1_dp
359 : END DO
360 198150 : spin_h(n + 1:, :n) = CONJG(TRANSPOSE(spin_h(:n, n + 1:)))
361 99078 : spin_s(:n, :n) = s0
362 99078 : spin_s(n + 1:, n + 1:) = s0
363 198162 : spin_c(:n, :n, :) = co
364 198162 : spin_c(n + 1:, n + 1:, :) = co
365 198162 : spin_sin(:n, :n, :) = si
366 198162 : 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 6 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
370 6 : IF (info /= 0 .OR. idx /= ABS(expected)) ERROR STOP 'Periodic spin-mixed Z2 phase failed'
371 394758 : 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 6 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
375 6 : IF (info /= 0 .OR. idx /= 0) ERROR STOP 'Periodic atomic Pfaffian calibration failed'
376 6 : DEALLOCATE (spin_h, spin_s, spin_c, spin_sin)
377 :
378 6 : IF (phase == 1) THEN
379 24 : 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 2 : t = 0.0_dp
382 18 : DO ky = 0, length - 1
383 146 : DO kx = 0, length - 1
384 128 : col = 2*(ky*length + kx) + 1
385 1168 : DO y = 0, length - 1
386 9344 : DO x = 0, length - 1
387 8192 : row = 2*(y*length + x) + 1
388 25600 : DO orbital = 0, 1
389 : t(row + orbital, col + orbital) = &
390 24576 : 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 33026 : ht(:, :) = MATMUL(CONJG(TRANSPOSE(t)), MATMUL(h0, t))
397 33026 : st(:, :) = s0
398 18 : DO ky = 0, length - 1
399 146 : DO kx = 0, length - 1
400 128 : col = 2*(ky*length + kx) + 1
401 384 : wave = two_pi*[REAL(kx, dp), REAL(ky, dp)]/length
402 384 : IF (ABS(ht(col, col) - (mass + SUM(COS(wave)))) > 1.e-12_dp) THEN
403 0 : ERROR STOP 'Periodic real-space/k-space Hamiltonian mismatch'
404 : END IF
405 144 : IF (ABS(ht(col, col + 1) - CMPLX(SIN(wave(1)), -SIN(wave(2)), dp)) > 1.e-12_dp) THEN
406 0 : ERROR STOP 'Periodic Bloch Hamiltonian offdiagonal mismatch'
407 : END IF
408 : END DO
409 : END DO
410 6 : DO j = 1, 2
411 66052 : ct(:, :, j) = MATMUL(CONJG(TRANSPOSE(t)), MATMUL(co(:, :, j), t))
412 66054 : 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 2 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
417 2 : IF (info /= 0 .OR. idx /= expected .OR. ABS(g - g0) > 1.e-11_dp) THEN
418 0 : ERROR STOP 'Full k-mesh/real-space periodic localizer mismatch'
419 : END IF
420 2 : t = 0.0_dp
421 258 : DO i = 1, n
422 256 : t(i, i) = (1.0_dp + MOD(i, 5)/5.0_dp)*EXP(CMPLX(0.0_dp, 0.17_dp*i, dp))
423 258 : IF (i < n) t(i, i + 1) = CMPLX(0.07_dp, -0.04_dp, dp)
424 : END DO
425 33026 : ht(:, :) = MATMUL(CONJG(TRANSPOSE(t)), MATMUL(h0, t))
426 33026 : st(:, :) = MATMUL(CONJG(TRANSPOSE(t)), t)
427 6 : DO j = 1, 2
428 66052 : ct(:, :, j) = MATMUL(CONJG(TRANSPOSE(t)), MATMUL(co(:, :, j), t))
429 66054 : 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 33026 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
434 2 : IF (info /= 0 .OR. idx /= expected .OR. ABS(g - g0) > 1.e-11_dp) THEN
435 0 : 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 33026 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
440 2 : IF (info /= 0 .OR. idx /= expected .OR. ABS(g - 2*g0) > 1.e-11_dp) THEN
441 0 : ERROR STOP 'Periodic energy-scale covariance failed'
442 : END IF
443 2 : 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 6 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
447 6 : IF (info /= -1) ERROR STOP 'Zero torus energy scale accepted'
448 6 : 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 6 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
451 6 : IF (info /= -1) ERROR STOP 'Non-Hermitian torus cosine accepted'
452 6 : h0 = 0.0_dp
453 6 : si = 0.0_dp
454 99078 : co(:, :, 1) = s0
455 99078 : 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 6 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, idx, g, cond, info)
458 6 : IF (info /= 1) ERROR STOP 'Closed periodic localizer assigned an index'
459 8 : DEALLOCATE (h0, s0, pos, co, si)
460 : END DO
461 2 : END SUBROUTINE periodic_tests
462 :
463 : ! **************************************************************************************************
464 : !> \brief Exercise 2x2 pivots, complex congruences, scale changes and singularity.
465 : ! **************************************************************************************************
466 2 : 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 2 : matrix = 0.0_dp
472 2 : matrix(1, 2) = CMPLX(0.0_dp, 2.0_dp, dp)
473 2 : matrix(2, 1) = CONJG(matrix(1, 2))
474 2 : matrix(3, 3) = 1.e-12_dp
475 2 : matrix(4, 4) = 3.0_dp
476 2 : change = 0.0_dp
477 10 : DO i = 1, 4
478 8 : change(i, i) = CMPLX(1.0_dp + 0.2_dp*i, 0.1_dp*i, dp)
479 22 : DO j = i + 1, 4
480 20 : change(i, j) = CMPLX(0.03_dp*j, -0.07_dp*i, dp)
481 : END DO
482 : END DO
483 338 : transformed = MATMUL(CONJG(TRANSPOSE(change)), MATMUL(matrix, change))
484 6 : DO exponent = -150, 150, 150
485 126 : CALL localizer_inertia((10.0_dp**exponent)*transformed, p, m, RESULT)
486 8 : IF (RESULT /= 0 .OR. p /= 3 .OR. m /= 1) ERROR STOP 'LDL block inertia or congruence failed'
487 : END DO
488 2 : matrix(3, 3) = -1.e-12_dp
489 2 : CALL localizer_inertia(matrix, p, m, RESULT)
490 2 : IF (RESULT /= 0 .OR. p /= 2 .OR. m /= 2) ERROR STOP 'LDL gap crossing failed'
491 2 : matrix(3, 3) = 0.0_dp
492 2 : CALL localizer_inertia(matrix, p, m, RESULT)
493 2 : IF (RESULT /= 1) ERROR STOP 'Singular LDL accepted'
494 2 : matrix(2, 1) = matrix(1, 2)
495 2 : CALL localizer_inertia(matrix, p, m, RESULT)
496 2 : IF (RESULT /= -1) ERROR STOP 'Complex symmetric matrix accepted as Hermitian'
497 2 : matrix = 0.0_dp
498 2 : matrix(1, 2) = CMPLX(0.0_dp, 1.e-12_dp, dp)
499 2 : matrix(2, 1) = CONJG(matrix(1, 2))
500 2 : matrix(3, 3) = 1.0_dp
501 2 : matrix(4, 4) = -1.0_dp
502 2 : CALL localizer_inertia(matrix, p, m, RESULT)
503 2 : IF (RESULT /= 0 .OR. p /= 2 .OR. m /= 2) ERROR STOP 'Small imaginary coupling discarded'
504 2 : END SUBROUTINE inertia_tests
505 :
506 : ! **************************************************************************************************
507 : !> \brief Cross-check real and complex LDL under congruences, all signatures and permutations.
508 : ! **************************************************************************************************
509 2 : SUBROUTINE real_inertia_tests()
510 2 : COMPLEX(dp), ALLOCATABLE :: cmatrix(:, :), phases(:)
511 : INTEGER :: exponent, i, info, j, m, n, negative, p
512 : REAL(dp) :: pair(2, 2), scale
513 2 : REAL(dp), ALLOCATABLE :: change(:, :), diagonal(:, :), &
514 2 : matrix(:, :)
515 :
516 26 : DO n = 1, 12
517 312 : ALLOCATE (diagonal(n, n), change(n, n), matrix(n, n), phases(n), cmatrix(n, n))
518 24 : change = 0.0_dp
519 180 : DO i = 1, n
520 156 : change(i, i) = 0.5_dp + REAL(i, dp)/n
521 156 : phases(i) = EXP(CMPLX(0.0_dp, 0.37_dp*i, dp))
522 752 : DO j = i + 1, n
523 728 : change(i, j) = 0.1_dp*SIN(REAL(i + 3*j, dp))
524 : END DO
525 : END DO
526 204 : DO negative = 0, n
527 180 : diagonal = 0.0_dp
528 1636 : DO i = 1, n
529 1456 : diagonal(i, i) = 1.0_dp + REAL(i, dp)/n
530 1636 : IF (i <= negative) diagonal(i, i) = -diagonal(i, i)
531 : END DO
532 564920 : matrix(:, :) = MATMUL(TRANSPOSE(change), MATMUL(diagonal, change))
533 1636 : DO i = 1, n
534 15104 : DO j = 1, n
535 14924 : cmatrix(i, j) = CONJG(phases(i))*matrix(i, j)*phases(j)
536 : END DO
537 : END DO
538 564 : DO exponent = -150, 150, 150
539 540 : scale = 10.0_dp**exponent
540 45312 : CALL localizer_inertia(scale*matrix, p, m, info)
541 540 : IF (info /= 0 .OR. p /= n - negative .OR. m /= negative) ERROR STOP 'Real LDL signature failed'
542 45312 : CALL localizer_inertia(scale*cmatrix, p, m, info)
543 540 : IF (info /= 0 .OR. p /= n - negative .OR. m /= negative) ERROR STOP 'Complex LDL signature failed'
544 45312 : CALL localizer_inertia(CMPLX(scale*matrix, 0.0_dp, dp), p, m, info)
545 720 : IF (info /= 0 .OR. p /= n - negative .OR. m /= negative) ERROR STOP 'Real LDL dispatch failed'
546 : END DO
547 : END DO
548 26 : DEALLOCATE (diagonal, change, matrix, phases, cmatrix)
549 : END DO
550 2 : pair = RESHAPE([0.0_dp, 2.0_dp, 2.0_dp, 0.0_dp], [2, 2])
551 2 : CALL localizer_inertia(pair, p, m, info)
552 2 : IF (info /= 0 .OR. p /= 1 .OR. m /= 1) ERROR STOP 'Zero-diagonal real 2x2 pivot failed'
553 2 : pair(1, 2) = 0.0_dp
554 2 : CALL localizer_inertia(pair, p, m, info)
555 2 : IF (info /= -1) ERROR STOP 'Nonsymmetric real matrix accepted'
556 2 : pair = 0.0_dp
557 2 : CALL localizer_inertia(pair, p, m, info)
558 2 : IF (info /= 1) ERROR STOP 'Singular real matrix accepted'
559 2 : END SUBROUTINE real_inertia_tests
560 :
561 : ! **************************************************************************************************
562 : !> \brief Generalized gaps have an analytic reference independent of LDL pivot magnitudes.
563 : ! **************************************************************************************************
564 2 : 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 2 : h = RESHAPE([1.0_dp, 20.0_dp, 20.0_dp, 399.0_dp], [2, 2])
571 2 : s = 0.0_dp
572 2 : s(1, 1) = 2.0_dp
573 2 : s(2, 2) = 2.0_dp
574 2 : r = 0.0_dp
575 2 : trace = REAL(h(1, 1) + h(2, 2), dp)
576 2 : 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 2 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, index, gap, ratio, status)
579 2 : IF (status /= 0 .OR. index /= 0 .OR. ABS(gap - expected) > 1.e-12_dp) ERROR STOP 'Analytic LDL/gap separation failed'
580 2 : change = RESHAPE([0.5_dp, 0.0_dp, 0.2_dp, 3.0_dp], [2, 2])
581 58 : hc = MATMUL(TRANSPOSE(change), MATMUL(h, change))
582 58 : 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 2 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, index, gap, ratio, status)
585 2 : IF (status /= 0 .OR. ABS(gap - expected) > 1.e-12_dp) ERROR STOP 'Metric LDL/gap separation failed'
586 2 : crossing = [-1.e-6_dp, -1.e-12_dp, 0.0_dp, 1.e-12_dp, 1.e-6_dp]
587 2 : h = 0.0_dp
588 2 : h(2, 2) = 5.0_dp
589 12 : DO i = 1, SIZE(crossing)
590 10 : 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 10 : 1.e-10_dp, 1.e-10_dp, 1.e-9_dp, index, gap, ratio, status)
593 12 : IF (i == 1 .OR. i == 5) THEN
594 4 : IF (status /= 0 .OR. ABS(gap - 5.e-7_dp) > 1.e-14_dp) ERROR STOP 'Resolved gap crossing failed'
595 : ELSE
596 6 : IF (status /= 1) ERROR STOP 'Unresolved gap assigned an index'
597 : END IF
598 : END DO
599 2 : END SUBROUTINE gap_tests
600 :
601 : ! **************************************************************************************************
602 : !> \brief Independent Pfaffian polynomials distinguish signs despite identical determinants.
603 : ! **************************************************************************************************
604 2 : SUBROUTINE pfaffian_orientation_tests()
605 : INTEGER :: exponent, i, j, n, pf_sign, sample, &
606 : status
607 2 : INTEGER, ALLOCATABLE :: permutation(:)
608 : REAL(dp) :: exact, scale
609 2 : REAL(dp), ALLOCATABLE :: a(:, :), change(:, :), transformed(:, :)
610 :
611 2 : DO n = 2, 8, 2
612 80 : ALLOCATE (a(n, n), change(n, n), transformed(n, n), permutation(n))
613 72 : DO sample = 1, 8
614 64 : a = 0.0_dp
615 320 : DO j = 2, n
616 1120 : DO i = 1, j - 1
617 1056 : a(i, j) = SIN(REAL(11*i + 17*j + sample, dp))
618 : END DO
619 : END DO
620 64 : DO i = 1, n, 2
621 160 : a(i, i + 1) = a(i, i + 1) + 3.0_dp
622 : END DO
623 4544 : a(:, :) = a - TRANSPOSE(a)
624 64 : exact = pfaffian_polynomial(a)
625 64 : IF (ABS(exact) < 1.e-6_dp) ERROR STOP 'Ill-conditioned Pfaffian fixture'
626 704 : permutation(:) = [(i, i=1, n)]
627 192 : permutation(:2) = [2, 1]
628 192 : DO exponent = -150, 150, 150
629 192 : scale = 10.0_dp**exponent
630 6912 : CALL localizer_pfaffian_sign(scale*a, pf_sign, status)
631 192 : IF (status /= 0 .OR. pf_sign /= INT(SIGN(1.0_dp, exact))) ERROR STOP 'Pfaffian polynomial sign failed'
632 6912 : CALL localizer_pfaffian_sign(-scale*a, pf_sign, status)
633 192 : IF (status /= 0 .OR. pf_sign /= (-1)**(n/2)*INT(SIGN(1.0_dp, exact))) ERROR STOP 'Negative Pfaffian scale failed'
634 6912 : CALL localizer_pfaffian_sign(scale*a(permutation, permutation), pf_sign, status)
635 256 : IF (status /= 0 .OR. pf_sign /= -INT(SIGN(1.0_dp, exact))) ERROR STOP 'Equal determinant/opposite Pfaffian failed'
636 : END DO
637 64 : change = 0.0_dp
638 384 : DO i = 1, n
639 320 : change(i, i) = 0.8_dp + 0.1_dp*i
640 384 : IF (i < n) change(i, i + 1) = 0.1_dp
641 : END DO
642 384 : change(1, :) = -change(1, :)
643 55936 : transformed(:, :) = MATMUL(change, MATMUL(a, TRANSPOSE(change)))
644 64 : CALL localizer_pfaffian_sign(transformed, pf_sign, status)
645 72 : IF (status /= 0 .OR. pf_sign /= -INT(SIGN(1.0_dp, exact))) ERROR STOP 'Pfaffian congruence orientation failed'
646 : END DO
647 8 : DEALLOCATE (a, change, transformed, permutation)
648 : END DO
649 2 : ALLOCATE (a(4, 4))
650 2 : a = 0.0_dp
651 2 : a(1, 2) = 1.0_dp
652 2 : a(1, 3) = 1.0_dp
653 2 : a(2, 4) = 1.0_dp
654 2 : a(3, 4) = 1.0_dp
655 82 : a(:, :) = a - TRANSPOSE(a)
656 2 : CALL localizer_pfaffian_sign(a, pf_sign, status)
657 2 : IF (status /= -2 .OR. pf_sign /= 0) ERROR STOP 'Pfaffian polynomial cancellation ignored'
658 4 : 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 2784 : RECURSIVE FUNCTION pfaffian_polynomial(a) RESULT(value)
666 : REAL(dp), INTENT(IN) :: a(:, :)
667 : REAL(dp) :: value
668 :
669 5568 : INTEGER :: i, j, n, keep(SIZE(a, 1) - 2)
670 :
671 2784 : n = SIZE(a, 1)
672 2784 : IF (n == 2) THEN
673 1984 : value = a(1, 2)
674 1984 : RETURN
675 : END IF
676 : value = 0.0_dp
677 3520 : DO j = 2, n
678 42272 : keep = PACK([(i, i=2, n)], [(i /= j, i=2, n)])
679 32832 : 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 44 : 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 44 : 1.e-9_dp, invariant, gap, ratio, status)
706 44 : 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 10 : 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 10 : n = 2*length**2
727 100 : ALLOCATE (h(n, n), s(n, n), r(n, n, 2))
728 10 : h = 0.0_dp
729 10 : s = 0.0_dp
730 10 : r = 0.0_dp
731 : tx = 0.0_dp
732 10 : tx(1, 1) = 0.5_dp
733 10 : tx(2, 2) = -0.5_dp
734 10 : tx(1, 2) = CMPLX(0.0_dp, -0.5_dp, dp)
735 10 : tx(2, 1) = tx(1, 2)
736 : ty = 0.0_dp
737 10 : ty(1, 1) = 0.5_dp
738 10 : ty(2, 2) = -0.5_dp
739 10 : ty(1, 2) = -0.5_dp
740 10 : ty(2, 1) = 0.5_dp
741 86 : DO y = 1, length
742 666 : DO x = 1, length
743 580 : i = 2*((y - 1)*length + x) - 1
744 580 : h(i, i) = mass
745 580 : h(i + 1, i + 1) = -mass
746 1740 : DO orbital = 0, 1
747 1160 : s(i + orbital, i + orbital) = 1.0_dp
748 1160 : r(i + orbital, i + orbital, 1) = REAL(x, dp) - REAL(length + 1, dp)/2.0_dp
749 1740 : r(i + orbital, i + orbital, 2) = REAL(y, dp) - REAL(length + 1, dp)/2.0_dp
750 : END DO
751 580 : IF (x < length) THEN
752 504 : j = i + 2
753 3528 : h(i:i + 1, j:j + 1) = tx
754 3528 : h(j:j + 1, i:i + 1) = CONJG(TRANSPOSE(tx))
755 : END IF
756 580 : IF (y < length) THEN
757 504 : j = i + 2*length
758 3528 : h(i:i + 1, j:j + 1) = ty
759 3528 : h(j:j + 1, i:i + 1) = CONJG(TRANSPOSE(ty))
760 : END IF
761 656 : IF (PRESENT(torus)) THEN
762 384 : IF (torus .AND. x == length) THEN
763 48 : j = 2*(y - 1)*length + 1
764 336 : h(i:i + 1, j:j + 1) = tx
765 336 : h(j:j + 1, i:i + 1) = CONJG(TRANSPOSE(tx))
766 : END IF
767 384 : IF (torus .AND. y == length) THEN
768 48 : j = 2*x - 1
769 336 : h(i:i + 1, j:j + 1) = ty
770 336 : h(j:j + 1, i:i + 1) = CONJG(TRANSPOSE(ty))
771 : END IF
772 : END IF
773 : END DO
774 : END DO
775 10 : 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 2 : 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 2 : io = cp_logger_get_default_io_unit()
804 2 : CALL determinant_controls()
805 2 : pauli = 0.0_dp
806 2 : pauli(1, 2, 1) = 1.0_dp
807 2 : pauli(2, 1, 1) = 1.0_dp
808 2 : pauli(1, 2, 2) = -ii
809 2 : pauli(2, 1, 2) = ii
810 2 : pauli(1, 1, 3) = 1.0_dp
811 2 : pauli(2, 2, 3) = -1.0_dp
812 8 : DO model = 1, 3
813 6 : mass = masses(model)
814 6 : CALL dirac_model(mass)
815 : ! Bulk inversion eigenvalues at the eight TRIM: strong, not weak, parity.
816 6 : sign_product = 1
817 18 : DO kx = -1, 1, 2
818 42 : DO ky = -1, 1, 2
819 36 : DO kz = -1, 1, 2
820 48 : sign_product = sign_product*INT(SIGN(1.0_dp, mass + kx + ky + kz))
821 : END DO
822 : END DO
823 : END DO
824 6 : expected = (1 - sign_product)/2
825 6 : CALL evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
826 6 : IF (status /= 0 .OR. index /= expected) ERROR STOP '3D strong/trivial parity mismatch'
827 6 : reference = gap
828 6 : a = 0.0_dp
829 70638 : a(:n, :n) = h
830 70638 : a(n + 1:, n + 1:) = h
831 24 : DO axis = 1, 3
832 60 : DO j = 1, 2
833 126 : DO i = 1, 2
834 : a((i - 1)*n + 1:i*n, (j - 1)*n + 1:j*n) = &
835 847692 : 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 6 : 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 1302 : IF (status /= 0 .OR. ABS(reference - MINVAL(singular)) > 2.e-12_dp) ERROR STOP '3D singular gap mismatch'
847 9 : IF (io > 0) WRITE (io, '(A,F4.1,I3,2ES18.9)') '3D Dirac mass, index, gap, SVD error: ', &
848 656 : 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 2 : u = 0.0_dp
853 110 : DO i = 1, n/2
854 108 : u(i, i) = CMPLX(0.8_dp + 0.004_dp*i, 0.07_dp*SIN(REAL(i, dp)), dp)
855 108 : u(i + n/2, i + n/2) = CONJG(u(i, i))
856 108 : u(i, i + n/2) = CMPLX(0.08_dp, 0.03_dp, dp)
857 108 : u(i + n/2, i) = -CONJG(u(i, i + n/2))
858 110 : IF (i < n/2) THEN
859 106 : u(i, i + 1) = CMPLX(0.02_dp, 0.01_dp, dp)
860 106 : u(i + n/2, i + n/2 + 1) = CONJG(u(i, i + 1))
861 : END IF
862 : END DO
863 23546 : hc = MATMUL(CONJG(TRANSPOSE(u)), MATMUL(h, u))
864 23546 : sc = MATMUL(CONJG(TRANSPOSE(u)), MATMUL(s, u))
865 8 : DO axis = 1, 3
866 70640 : rc(:, :, axis) = MATMUL(CONJG(TRANSPOSE(u)), MATMUL(r(:, :, axis), u))
867 : END DO
868 2 : CALL evaluate(hc, sc, rc, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
869 2 : defect = ABS(gap - reference)
870 2 : IF (status /= 0 .OR. index /= 1 .OR. defect > 3.e-12_dp) ERROR STOP '3D nonorthogonal spin-mixing basis'
871 2 : IF (io > 0) WRITE (io, '(A,ES18.9)') '3D nonorthogonal spin-mixing gap residual: ', defect
872 2 : shift = [1.2_dp, -0.8_dp, 0.4_dp]
873 8 : DO axis = 1, 3
874 70640 : rc(:, :, axis) = r(:, :, axis) + shift(axis)*s
875 : END DO
876 23546 : CALL evaluate(h + 3.0_dp*s, s, rc, shift, 3.0_dp)
877 2 : IF (status /= 0 .OR. index /= 1 .OR. ABS(gap - reference) > 3.e-12_dp) ERROR STOP '3D origin shifts'
878 23546 : rc(:, :, 1) = 0.6_dp*r(:, :, 1) + 0.8_dp*r(:, :, 2)
879 23546 : rc(:, :, 2) = -0.8_dp*r(:, :, 1) + 0.6_dp*r(:, :, 2)
880 23546 : rc(:, :, 3) = -r(:, :, 3)
881 2 : CALL evaluate(h, s, rc, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
882 2 : 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 2 : 1.e-10_dp, 1.e-10_dp, 1.e-10_dp, index, gap, ratio, status)
886 2 : IF (status /= 0 .OR. index /= 0) ERROR STOP '3D finite-size counterexample lost'
887 : lower = 0.3_dp
888 : upper = 0.5_dp
889 62 : DO i = 1, 32
890 62 : 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 62 : 1.e-10_dp, 1.e-10_dp, 0.0_dp, index, gap, ratio, status)
893 62 : IF (status == 1) EXIT
894 60 : IF (status /= 0) ERROR STOP '3D transition factorization failed'
895 62 : IF (index == 0) THEN
896 : lower = midpoint
897 : ELSE
898 40 : 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 2 : 1.e-10_dp, 1.e-10_dp, 1.e-8_dp, index, gap, ratio, status)
903 2 : IF (status /= 1) ERROR STOP '3D closing gap assigned an index'
904 2 : 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 2 : 1.e-10_dp, -1.0_dp, 1.e-10_dp, index, gap, ratio, status)
907 2 : 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 2 : 1.e-10_dp, 1.e-10_dp, ieee_value(0.0_dp, ieee_quiet_nan), index, gap, ratio, status)
910 2 : IF (status /= -1) ERROR STOP '3D nonfinite gap tolerance accepted'
911 2 : h(1, 1) = h(1, 1) + 0.1_dp
912 2 : CALL evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
913 2 : IF (status /= -3) ERROR STOP '3D broken time reversal accepted'
914 2 : h = 0.0_dp
915 2 : r = 0.0_dp
916 2 : CALL evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
917 2 : IF (status /= 1 .OR. index /= 0) ERROR STOP '3D zero gap accepted'
918 2 : s(1, 1) = -1.0_dp
919 2 : s(n/2 + 1, n/2 + 1) = -1.0_dp
920 2 : CALL evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
921 2 : IF (status /= -2) ERROR STOP '3D indefinite metric accepted'
922 2 : h(1, 1) = ieee_value(0.0_dp, ieee_quiet_nan)
923 2 : CALL evaluate(h, s, r, [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp)
924 2 : 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 3 : 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 2 : 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 2 : metric = 0.0_dp
948 66 : DO row = 1, m/2
949 66 : metric(row, row) = 1.0_dp
950 : END DO
951 2 : matrix = 0.0_dp
952 130 : DO row = 1, m
953 128 : matrix(row, row) = 1.0_dp
954 4160 : matrix(row, :row - 1) = -1.0_dp
955 130 : matrix(row, m) = 1.0_dp + 0.25_dp*MODULO(row, 3)
956 : END DO
957 6 : DO orientation = 0, 1
958 : ! Positive last-column entries give a positive determinant. One row reversal changes its sign.
959 134 : IF (orientation == 1) matrix(1, :) = -matrix(1, :)
960 : ! Exact dyadic entries in Q W Q^dagger retain the pivot-growth fixture after conversion.
961 260 : DO col = 1, m
962 256 : paired_col = (3 - (col - 1)/quarter)*quarter + MODULO(col - 1, quarter) + 1
963 256 : col_sign = MERGE(1, -1, col <= quarter .OR. col > 3*quarter)
964 16644 : DO row = 1, m
965 16384 : paired_row = (3 - (row - 1)/quarter)*quarter + MODULO(row - 1, quarter) + 1
966 16384 : 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 16640 : 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 4 : RESULT, measured_gap, metric_condition, info)
973 4 : 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 4 : SIZE(workspace), real_work, info)
979 : #if defined(__HAS_IEEE_EXCEPTIONS)
980 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
981 : #endif
982 260 : svd_gap = MINVAL(values)
983 4 : IF (info /= 0 .OR. ABS(measured_gap - svd_gap) > 1.e-12_dp) ERROR STOP '3D high-growth singular gap'
984 6 : IF (io > 0) WRITE (io, '(A,I2,2ES18.9)') '3D high-growth determinant, gap, SVD error: ', &
985 6 : RESULT, measured_gap, ABS(measured_gap - svd_gap)
986 : END DO
987 : ! Independent atomic limit for both signs of a spin-degenerate energy.
988 4 : DO col = -1, 1, 2
989 4 : block = 0.0_dp
990 260 : DO row = 1, m
991 260 : 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 4 : RESULT, measured_gap, metric_condition, info)
995 6 : IF (info /= 0 .OR. RESULT /= 0 .OR. ABS(measured_gap - 1.0_dp) > 1.e-12_dp) THEN
996 0 : ERROR STOP '3D atomic reference orientation'
997 : END IF
998 : END DO
999 2 : 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 20 : 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 20 : index, gap, ratio, status)
1015 20 : 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 6 : 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 6 : h = 0.0_dp
1027 6 : s = 0.0_dp
1028 6 : r = 0.0_dp
1029 24 : DO z = 0, side - 1
1030 78 : DO y = 0, side - 1
1031 234 : DO x = 0, side - 1
1032 648 : pos = [x, y, z]
1033 162 : site = 1 + x + side*y + side**2*z
1034 540 : DO spin = 1, 2
1035 1134 : DO orbital = 1, 2
1036 648 : row = site + sites*(orbital - 1) + 2*sites*(spin - 1)
1037 648 : s(row, row) = 1.0_dp
1038 648 : h(row, row) = m*REAL(3 - 2*orbital, dp)
1039 2592 : r(row, row, :) = REAL(pos - 1, dp)
1040 2916 : DO direction = 1, 3
1041 1944 : IF (pos(direction) == side - 1) CYCLE
1042 1296 : neighbor = site + side**(direction - 1)
1043 1296 : col = neighbor + sites*(orbital - 1) + 2*sites*(spin - 1)
1044 1296 : h(row, col) = 0.5_dp*REAL(3 - 2*orbital, dp)
1045 1296 : h(col, row) = CONJG(h(row, col))
1046 4536 : DO other = 1, 2
1047 2592 : col = neighbor + sites*(2 - orbital) + 2*sites*(other - 1)
1048 2592 : h(row, col) = -0.5_dp*ii*pauli(spin, other, direction)
1049 4536 : 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 6 : 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 2 : 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 2 : 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 2 : 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 4 : io = cp_logger_get_default_io_unit()
1097 2 : CALL dbcsr_init_lib(para_env%get_handle(), io)
1098 2 : CALL init_fft('FFTSG', alltoall=.FALSE., pool_limit=4, wisdom_file='', plan_style=1)
1099 2 : nr = para_env%num_pe
1100 2 : nc = 1
1101 4 : DO i = 1, command_argument_count()
1102 2 : CALL get_command_argument(i, argument)
1103 4 : IF (TRIM(argument) == '--column-grid') THEN
1104 0 : nr = 1
1105 0 : nc = para_env%num_pe
1106 : END IF
1107 : END DO
1108 8 : ALLOCATE (process_grid(0:nr - 1, 0:nc - 1), row_dist(nb), col_dist(nb), sizes(nb))
1109 14 : process_grid = RESHAPE([(i, i=0, para_env%num_pe - 1)], [nr, nc])
1110 26 : row_dist = [(MODULO(i - 1, nr), i=1, nb)]
1111 26 : col_dist = [(MODULO(i - 1, nc), i=1, nb)]
1112 14 : sizes = bs
1113 : CALL dbcsr_distribution_new(distribution, group=para_env%get_handle(), pgrid=process_grid, &
1114 2 : row_dist=row_dist, col_dist=col_dist)
1115 2 : energy = 0.002_dp
1116 2 : scale = 1.7_dp
1117 2 : maximum_error = 0.0_dp
1118 14 : DO mode = 0, 5
1119 12 : u = 0.0_dp
1120 12 : transform = 0.0_dp
1121 156 : DO i = 1, n
1122 144 : energies(i) = 0.01_dp*(i - n/2) - 0.005_dp
1123 144 : transform(i, i) = 1.0_dp
1124 1872 : DO j = 1, n
1125 1872 : 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 144 : IF (mode == 0) CYCLE
1128 120 : transform(i, i) = 0.5_dp + 0.1_dp*i
1129 132 : IF (i < n) transform(i, i + 1) = CMPLX(0.13_dp, -0.08_dp, dp)
1130 : END DO
1131 12 : IF (mode == 2) transform(n, n) = 0.1_dp
1132 36 : IF (mode == 4) energies = energies + 0.5_dp
1133 36 : IF (mode == 5) energies = energies - 0.5_dp
1134 12 : IF (mode == 3) THEN
1135 2 : u = 0.0_dp
1136 2 : rotation = 0.0_dp
1137 86 : transform(:n/2, n/2 + 1:) = 0.0_dp
1138 86 : transform(n/2 + 1:, :n/2) = 0.0_dp
1139 86 : transform(n/2 + 1:, n/2 + 1:) = CONJG(transform(:n/2, :n/2))
1140 14 : DO i = 1, n/2
1141 12 : energies(i) = 0.01_dp*(i - n/4) - 0.005_dp
1142 12 : energies(n/2 + i) = energies(i)
1143 12 : rotation(i, i) = COS(0.37_dp)
1144 12 : rotation(n/2 + i, n/2 + i) = COS(0.37_dp)
1145 12 : rotation(i, n/2 + i) = SIN(0.37_dp)
1146 12 : rotation(n/2 + i, i) = -SIN(0.37_dp)
1147 86 : DO j = 1, n/2
1148 84 : 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 86 : u(n/2 + 1:, n/2 + 1:) = CONJG(u(:n/2, :n/2))
1152 3770 : u = MATMUL(rotation, u)
1153 : END IF
1154 1884 : scaled = u*SPREAD(energies, 1, n)
1155 67836 : hdense = MATMUL(CONJG(TRANSPOSE(transform)), MATMUL(MATMUL(scaled, CONJG(TRANSPOSE(u))), transform))
1156 22620 : sdense = MATMUL(CONJG(TRANSPOSE(transform)), transform)
1157 2028 : scaled = u*SPREAD(SIGN(scale, energies - energy), 1, n)
1158 67836 : exact = MATMUL(CONJG(TRANSPOSE(transform)), MATMUL(MATMUL(scaled, CONJG(TRANSPOSE(u))), transform))
1159 12 : CALL make_operator(hdense, h)
1160 12 : CALL make_operator(sdense, s)
1161 12 : CALL localizer_flatten_dbcsr(h, s, energy, scale, 1.e-9_dp, para_env, flat, residuals, status)
1162 12 : IF (status /= 0) THEN
1163 0 : IF (io > 0) WRITE (io, *) 'Flatten failure: ', mode, status, residuals
1164 0 : CPABORT('Distributed spectral flattening failed')
1165 : END IF
1166 12 : CALL localizer_dbcsr_dense(flat, para_env, RESULT)
1167 1884 : difference = MAXVAL(ABS(RESULT - exact))
1168 12 : maximum_error = MAX(maximum_error, difference)
1169 12 : IF (difference > 1.e-9_dp) THEN
1170 0 : CPABORT('Covariant sign differs from exact eigensystem')
1171 : END IF
1172 12 : IF (mode == 3) THEN
1173 86 : IF (MAXVAL(ABS(RESULT(:n/2, :n/2) - CONJG(RESULT(n/2 + 1:, n/2 + 1:)))) > 1.e-9_dp) THEN
1174 0 : CPABORT('Flattening broke physical TR diagonal blocks')
1175 : END IF
1176 86 : IF (MAXVAL(ABS(RESULT(:n/2, n/2 + 1:) + CONJG(RESULT(n/2 + 1:, :n/2)))) > 1.e-9_dp) THEN
1177 0 : CPABORT('Flattening broke physical TR spin mixing')
1178 : END IF
1179 : END IF
1180 12 : DEALLOCATE (RESULT)
1181 12 : CALL localizer_dbcsr_dense(h, para_env, RESULT)
1182 1884 : IF (MAXVAL(ABS(RESULT - hdense)) > 1.e-14_dp) THEN
1183 0 : CPABORT('Flattening modified input Hamiltonian')
1184 : END IF
1185 12 : DEALLOCATE (RESULT)
1186 12 : CALL localizer_dbcsr_dense(s, para_env, RESULT)
1187 1884 : IF (MAXVAL(ABS(RESULT - sdense)) > 1.e-14_dp) THEN
1188 0 : CPABORT('Flattening modified input metric')
1189 : END IF
1190 12 : DEALLOCATE (RESULT)
1191 12 : 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 12 : para_env, flat, gap, metric_ratio, residuals, status)
1194 12 : IF (status /= 0) CPABORT('Checked flattening failed')
1195 156 : exact_gap = MINVAL(ABS(energies - energy))
1196 12 : IF (gap(1) > exact_gap + 1.e-10_dp .OR. gap(2) < exact_gap - 1.e-10_dp) THEN
1197 0 : CPABORT('Input electronic gap differs from reference')
1198 : END IF
1199 12 : IF (metric_ratio <= 1.e-10_dp) CPABORT('Invalid positive metric bound')
1200 12 : CALL localizer_dbcsr_dense(flat, para_env, RESULT)
1201 1884 : difference = MAXVAL(ABS(RESULT - exact))
1202 12 : maximum_error = MAX(maximum_error, difference)
1203 12 : IF (difference > 1.e-9_dp) THEN
1204 0 : CPABORT('Checked flattening changed the exact result')
1205 : END IF
1206 12 : DEALLOCATE (RESULT)
1207 12 : CALL localizer_dbcsr_release(flat)
1208 12 : CALL localizer_dbcsr_release(h)
1209 14 : CALL localizer_dbcsr_release(s)
1210 : END DO
1211 :
1212 2 : hdense = 0.0_dp
1213 2 : sdense = 0.0_dp
1214 26 : DO i = 1, n
1215 24 : sdense(i, i) = 1.0_dp
1216 26 : hdense(i, i) = 1.0_dp
1217 : END DO
1218 2 : hdense(1, 1) = 0.0_dp
1219 2 : CALL make_operator(hdense, h)
1220 2 : CALL make_operator(sdense, s)
1221 2 : CALL localizer_flatten_dbcsr(h, s, 0.0_dp, scale, 1.e-9_dp, para_env, flat, residuals, status)
1222 2 : 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 2 : para_env, flat, gap, metric_ratio, residuals, status)
1225 2 : IF (status /= 1) CPABORT('Checked flattening accepted a zero spectral mode')
1226 2 : CALL localizer_flatten_dbcsr(h, s, energy, -scale, 1.e-9_dp, para_env, flat, residuals, status)
1227 2 : 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 2 : para_env, flat, residuals, status)
1230 2 : IF (status /= -1) CPABORT('Nonfinite energy was accepted')
1231 2 : CALL localizer_dbcsr_release(h)
1232 2 : hdense = 0.0_dp
1233 2 : CALL make_operator(hdense, h)
1234 2 : CALL localizer_flatten_dbcsr(h, s, 0.0_dp, scale, 1.e-9_dp, para_env, flat, residuals, status)
1235 2 : IF (status /= 1) CPABORT('Zero Hamiltonian sign was accepted')
1236 2 : CALL localizer_dbcsr_release(h)
1237 2 : hdense(1, 2) = 1.0_dp
1238 2 : CALL make_operator(hdense, h)
1239 2 : CALL localizer_flatten_dbcsr(h, s, energy, scale, 1.e-9_dp, para_env, flat, residuals, status)
1240 2 : IF (status /= -1) CPABORT('Non-Hermitian Hamiltonian was accepted')
1241 2 : CALL localizer_dbcsr_release(h)
1242 2 : CALL localizer_dbcsr_release(s)
1243 2 : hdense = sdense
1244 2 : hdense(1, 1) = 1.e-11_dp
1245 2 : CALL make_operator(hdense, h)
1246 2 : 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 2 : para_env, flat, gap, metric_ratio, residuals, status)
1249 2 : IF (status /= 1) CPABORT('Unresolved near-zero spectral mode was accepted')
1250 2 : CALL localizer_dbcsr_release(s)
1251 2 : sdense(1, 1) = -1.0_dp
1252 2 : 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 2 : para_env, flat, gap, metric_ratio, residuals, status)
1255 2 : IF (status /= -2) CPABORT('Indefinite metric was accepted')
1256 2 : CALL localizer_dbcsr_release(s)
1257 2 : sdense(1, 1) = 1.e-12_dp
1258 2 : 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 2 : para_env, flat, gap, metric_ratio, residuals, status)
1261 2 : IF (status /= -2) CPABORT('Ill-conditioned metric was accepted')
1262 2 : CALL localizer_dbcsr_release(h)
1263 2 : CALL localizer_dbcsr_release(s)
1264 2 : CALL check_bloch()
1265 2 : IF (io > 0) WRITE (io, '(A,ES14.5)') 'Spectral flattening exact-matrix maximum error: ', maximum_error
1266 2 : CALL dbcsr_distribution_release(distribution)
1267 2 : DEALLOCATE (process_grid, row_dist, col_dist, sizes)
1268 2 : CALL dbcsr_finalize_lib()
1269 4 : 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 2 : 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 2 : bloch_error = 0.0_dp
1284 16 : DO m = 1, 7
1285 98 : sizes = bs
1286 14 : grid = [3, 1, 1]
1287 14 : nspin = 2
1288 14 : IF (m == 2) grid = [1, 3, 1]
1289 14 : IF (m == 3) grid = [1, 1, 3]
1290 14 : IF (m == 4) THEN
1291 2 : grid = [2, 3, 1]
1292 2 : nspin = 1
1293 : END IF
1294 14 : IF (m == 5) THEN
1295 2 : nspin = 1
1296 14 : sizes = [1, 3, 1, 3, 1, 3]
1297 : END IF
1298 20 : IF (m == 6) grid = 1
1299 56 : ncell = PRODUCT(grid)
1300 14 : nscalar = n/(nspin*ncell)
1301 14 : nsmall = nscalar*nspin
1302 14 : u = 0.0_dp
1303 14 : transform = 0.0_dp
1304 182 : DO row = 1, n
1305 168 : t = MODULO((row - 1)/nscalar, ncell)
1306 168 : p = ((row - 1)/(ncell*nscalar))*nscalar + MODULO(row - 1, nscalar)
1307 672 : r = [MODULO(t, grid(1)), MODULO(t/grid(1), grid(2)), t/(grid(1)*grid(2))]
1308 168 : transform(row, row) = 1.0_dp + 0.03_dp*p
1309 168 : pnext = MODULO(p + 1, nsmall)
1310 168 : rnext = (pnext/nscalar)*ncell*nscalar + t*nscalar + MODULO(pnext, nscalar) + 1
1311 168 : transform(row, rnext) = transform(row, rnext) + CMPLX(0.07_dp, -0.04_dp, dp)
1312 168 : tnext = MODULO(t + 1, ncell)
1313 168 : 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 168 : IF (m /= 4) transform(row, rnext) = transform(row, rnext) + CMPLX(0.05_dp, 0.02_dp, dp)
1316 710 : DO k = 0, ncell - 1
1317 2112 : kr = [MODULO(k, grid(1)), MODULO(k/grid(1), grid(2)), k/(grid(1)*grid(2))]
1318 2712 : DO band = 0, nsmall - 1
1319 8064 : phase = twopi*(SUM(REAL(kr*r, dp)/REAL(grid, dp)) + REAL(p*band, dp)/REAL(nsmall, dp))
1320 2016 : u(row, k*nsmall + band + 1) = EXP(CMPLX(0.0_dp, phase, dp))/SQRT(REAL(n, dp))
1321 2016 : energies(k*nsmall + band + 1) = 0.2_dp + 0.03_dp*band + 0.01_dp*COS(REAL(k, dp))
1322 2016 : IF (m == 7) energies(k*nsmall + band + 1) = 0.25_dp + 0.01_dp*COS(REAL(k, dp))
1323 2544 : IF (band < nsmall/2) energies(k*nsmall + band + 1) = -energies(k*nsmall + band + 1)
1324 : END DO
1325 : END DO
1326 : END DO
1327 2198 : scaled = u*SPREAD(energies, 1, n)
1328 79142 : hdense = MATMUL(CONJG(TRANSPOSE(transform)), MATMUL(MATMUL(scaled, CONJG(TRANSPOSE(u))), transform))
1329 26390 : sdense = MATMUL(CONJG(TRANSPOSE(transform)), transform)
1330 2366 : scaled = u*SPREAD(SIGN(scale, energies - energy), 1, n)
1331 79142 : exact = MATMUL(CONJG(TRANSPOSE(transform)), MATMUL(MATMUL(scaled, CONJG(TRANSPOSE(u))), transform))
1332 14 : CALL make_operator(hdense, h)
1333 14 : 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 14 : para_env, flat, bloch_gap, metric_ratio, residuals, status)
1336 14 : IF (status /= 0) CPABORT('Bloch flattening failed')
1337 14 : CALL localizer_dbcsr_dense(flat, para_env, RESULT)
1338 2198 : norm_error = MAXVAL(ABS(RESULT - exact))
1339 14 : bloch_error = MAX(bloch_error, norm_error)
1340 14 : maximum_error = MAX(maximum_error, norm_error)
1341 182 : ref_gap = MINVAL(ABS(energies - energy))
1342 14 : IF (norm_error > 1.e-11_dp .OR. ABS(bloch_gap - ref_gap) > 1.e-12_dp) THEN
1343 0 : CPABORT('Bloch flattening differs from exact full-torus reference')
1344 : END IF
1345 14 : 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 14 : para_env, flat, gap, metric_ratio, residuals, status)
1348 14 : IF (status /= 0) CPABORT('Full-torus reference flattening failed')
1349 14 : CALL localizer_dbcsr_dense(flat, para_env, RESULT)
1350 2198 : IF (MAXVAL(ABS(RESULT - exact)) > 1.e-10_dp) THEN
1351 0 : CPABORT('Full-torus/Bloch sign mismatch')
1352 : END IF
1353 14 : CALL localizer_dbcsr_release(flat)
1354 14 : CALL localizer_dbcsr_release(h)
1355 14 : IF (m /= 6) THEN
1356 12 : hdense(1, 1) = hdense(1, 1) + 0.01_dp
1357 12 : 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 12 : para_env, flat, bloch_gap, metric_ratio, residuals, status)
1360 12 : IF (status /= -7) CPABORT('Broken translations accepted')
1361 12 : CALL localizer_dbcsr_release(h)
1362 : END IF
1363 2198 : 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 14 : para_env, flat, bloch_gap, metric_ratio, residuals, status)
1366 14 : IF (status /= 1) CPABORT('Closed Bloch gap accepted')
1367 14 : CALL localizer_dbcsr_release(s)
1368 2198 : 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 14 : para_env, flat, bloch_gap, metric_ratio, residuals, status)
1371 14 : IF (status /= -2) CPABORT('Bloch flattening accepted an indefinite metric')
1372 14 : CALL localizer_dbcsr_release(h)
1373 16 : CALL localizer_dbcsr_release(s)
1374 : END DO
1375 14 : sizes = bs
1376 : ! Each Bloch metric is well conditioned, but their union has an unacceptable ratio.
1377 2 : hdense = 0.0_dp
1378 2 : sdense = 0.0_dp
1379 26 : DO row = 1, n
1380 24 : hdense(row, row) = 1.0_dp
1381 24 : sdense(row, row) = 0.5_dp*(1.0_dp + 1.e-12_dp)
1382 24 : rnext = 1 + MODULO(row - 1 + n/2, n)
1383 26 : sdense(row, rnext) = 0.5_dp*(1.0_dp - 1.e-12_dp)
1384 : END DO
1385 2 : CALL make_operator(hdense, h)
1386 2 : 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 2 : para_env, flat, bloch_gap, metric_ratio, residuals, status)
1389 2 : 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 2 : para_env, flat, bloch_gap, metric_ratio, residuals, status)
1392 2 : IF (status /= -1) CPABORT('Invalid Bloch grid accepted')
1393 2 : CALL localizer_dbcsr_release(h)
1394 2 : CALL localizer_dbcsr_release(s)
1395 2 : IF (io > 0) WRITE (io, '(A,ES14.5)') 'Seven Bloch flattening pencils; maximum error: ', bloch_error
1396 2 : 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 112 : 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 112 : CALL dbcsr_create(matrix%real_part, 'FLATTEN REAL', distribution, dbcsr_type_no_symmetry, sizes, sizes)
1410 112 : CALL dbcsr_create(matrix%imaginary_part, template=matrix%real_part)
1411 784 : DO row = 1, nb
1412 2352 : first = SUM(sizes(:row - 1)) + 1
1413 3024 : last = SUM(sizes(:row))
1414 4816 : DO col = 1, nb
1415 4032 : IF (process_grid(row_dist(row), col_dist(col)) /= para_env%mepos) CYCLE
1416 7056 : left = SUM(sizes(:col - 1)) + 1
1417 9072 : right = SUM(sizes(:col))
1418 14112 : CALL dbcsr_put_block(matrix%real_part, row, col, REAL(dense(first:last, left:right), dp))
1419 16800 : CALL dbcsr_put_block(matrix%imaginary_part, row, col, AIMAG(dense(first:last, left:right)))
1420 : END DO
1421 : END DO
1422 112 : CALL dbcsr_finalize(matrix%real_part)
1423 112 : CALL dbcsr_finalize(matrix%imaginary_part)
1424 112 : END SUBROUTINE make_operator
1425 :
1426 : END SUBROUTINE check_localizer_flattening
1427 :
1428 338 : END MODULE spectral_localizer_test_cases
1429 :
1430 2 : PROGRAM spectral_localizer_unittest
1431 : USE cp_log_handling, ONLY: cp_add_default_logger, &
1432 : cp_logger_create, &
1433 : cp_logger_release, &
1434 2 : 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 2 : ALLOCATE (para_env)
1445 2 : CALL mp_world_init(para_env)
1446 2 : io = -1
1447 2 : IF (para_env%is_source()) io = default_output_unit
1448 2 : CALL cp_logger_create(logger, para_env=para_env, default_global_unit_nr=io, close_global_unit_on_dealloc=.FALSE.)
1449 2 : CALL cp_add_default_logger(logger)
1450 2 : CALL check_localizer_2d()
1451 2 : CALL check_localizer_3d()
1452 2 : CALL check_localizer_flattening(para_env)
1453 2 : CALL cp_rm_default_logger()
1454 2 : CALL cp_logger_release(logger)
1455 2 : DEALLOCATE (para_env)
1456 2 : CALL mp_world_finalize()
1457 2 : END PROGRAM spectral_localizer_unittest
|
