Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Local Wilson-loop linear algebra, independent of the electronic-structure representation.
10 : !> WCC are arg(eigenvalue)/(2*pi), as in Z2Pack. The occupied subspace must be isolated.
11 : ! **************************************************************************************************
12 : MODULE topology_wilson
13 : USE ieee_arithmetic, ONLY: ieee_is_finite
14 : USE kinds, ONLY: dp
15 :
16 : IMPLICIT NONE
17 : PRIVATE
18 : PUBLIC :: wilson_step, wilson_spectrum, z2_from_wcc, wcc_distance, surface_resolved, chern_from_wcc
19 : REAL(KIND=dp), PARAMETER :: two_pi = 6.283185307179586476925286766559_dp
20 : ! Reduced-coordinate tolerance for an ambiguous largest-gap crossing.
21 : REAL(KIND=dp), PARAMETER :: crossing_tol = 1.e-12_dp
22 : ! Fraction of the adjacent largest gap allowed for motion and gap separation.
23 : REAL(KIND=dp), PARAMETER :: gap_fraction = 0.3_dp
24 : CONTAINS
25 :
26 : ! **************************************************************************************************
27 : !> \brief First Chern number from determinant Wilson-phase winding on a closed surface.
28 : !> Both transverse endpoints must be present. This is not a Z2 half-surface.
29 : !> The sign follows Z2Pack: increasing transverse coordinate, fixed loop orientation.
30 : !> \param wcc Wilson centres, band index first and transverse line second
31 : !> \param invariant candidate integer, meaningful only for status zero
32 : !> \param winding unrounded winding
33 : !> \param status 0 success, -1 invalid data, -2 unclosed surface, -3 unresolved phase step
34 : !> \param closure_tol maximal endpoint WCC mismatch
35 : ! **************************************************************************************************
36 28 : SUBROUTINE chern_from_wcc(wcc, invariant, winding, status, closure_tol)
37 : REAL(KIND=dp), INTENT(IN) :: wcc(:, :)
38 : INTEGER, INTENT(OUT) :: invariant
39 : REAL(KIND=dp), INTENT(OUT) :: winding
40 : INTEGER, INTENT(OUT) :: status
41 : REAL(KIND=dp), INTENT(IN) :: closure_tol
42 :
43 : REAL(KIND=dp), PARAMETER :: max_phase_step = 0.25_dp
44 :
45 : INTEGER :: i, nline
46 : REAL(KIND=dp) :: delta
47 :
48 28 : invariant = 0
49 28 : winding = 0.0_dp
50 28 : status = -1
51 28 : nline = SIZE(wcc, 2)
52 28 : IF (SIZE(wcc, 1) < 1 .OR. nline < 3 .OR. closure_tol <= 0.0_dp) RETURN
53 2850 : IF (.NOT. ALL(ieee_is_finite(wcc))) RETURN
54 28 : status = -2
55 28 : IF (wcc_distance(wcc(:, 1), wcc(:, nline)) > closure_tol) RETURN
56 26 : status = -3
57 1198 : DO i = 2, nline
58 : ! Avoid accepting a small aliased determinant step when many centres move.
59 1174 : IF (SIZE(wcc, 1)*wcc_distance(wcc(:, i), wcc(:, i - 1)) >= max_phase_step) RETURN
60 4284 : delta = MODULO(SUM(wcc(:, i)) - SUM(wcc(:, i - 1)) + 0.5_dp, 1.0_dp) - 0.5_dp
61 1172 : IF (ABS(delta) >= max_phase_step) RETURN
62 1196 : winding = winding + delta
63 : END DO
64 24 : status = -2
65 24 : IF (ABS(winding - NINT(winding)) > closure_tol) RETURN
66 24 : invariant = NINT(winding)
67 24 : status = 0
68 : END SUBROUTINE chern_from_wcc
69 :
70 : ! **************************************************************************************************
71 : !> \brief Conservative neighbouring-line movement and largest-gap separation checks.
72 : !> \param wcc centres, band index first and surface-line index second
73 : !> \return whether the sampled surface is locally resolved (not a proof between samples)
74 : ! **************************************************************************************************
75 6 : FUNCTION surface_resolved(wcc) RESULT(resolved)
76 : REAL(KIND=dp), INTENT(IN) :: wcc(:, :)
77 : LOGICAL :: resolved
78 :
79 : INTEGER :: i, j, n
80 : REAL(KIND=dp) :: delta, left_gap, left_size, right_gap, &
81 : right_size
82 6 : REAL(KIND=dp), ALLOCATABLE :: sorted(:, :)
83 :
84 6 : resolved = .FALSE.
85 6 : n = SIZE(wcc, 1)
86 6 : IF (n < 1 .OR. SIZE(wcc, 2) < 2) RETURN
87 24 : ALLOCATE (sorted(n, SIZE(wcc, 2)))
88 132 : sorted(:, :) = MODULO(wcc, 1.0_dp)
89 28 : DO i = 1, SIZE(wcc, 2)
90 28 : CALL sort_wcc(sorted(:, i))
91 : END DO
92 18 : DO i = 2, SIZE(wcc, 2)
93 14 : CALL largest_gap(sorted(:, i - 1), left_gap, left_size)
94 14 : CALL largest_gap(sorted(:, i), right_gap, right_size)
95 14 : IF (wcc_distance(sorted(:, i - 1), sorted(:, i)) >= gap_fraction*MIN(left_size, right_size)) RETURN
96 88 : DO j = 1, n
97 72 : delta = ABS(sorted(j, i) - left_gap)
98 72 : IF (MIN(delta, 1.0_dp - delta) <= gap_fraction*left_size) RETURN
99 72 : delta = ABS(sorted(j, i - 1) - right_gap)
100 84 : IF (MIN(delta, 1.0_dp - delta) <= gap_fraction*right_size) RETURN
101 : END DO
102 : END DO
103 6 : resolved = .TRUE.
104 6 : END FUNCTION surface_resolved
105 :
106 : ! **************************************************************************************************
107 : !> \brief Multiply by the unitary polar factor of an overlap. Reject rank-deficient links.
108 : !> \param product running Wilson matrix
109 : !> \param overlap overlap between adjacent occupied subspaces
110 : !> \param minimum_sv smallest singular value of this link
111 : !> \param status zero on success; negative for invalid/singular input; positive LAPACK failure
112 : !> \param sv_tol rank tolerance
113 : ! **************************************************************************************************
114 89956 : SUBROUTINE wilson_step(product, overlap, minimum_sv, status, sv_tol)
115 : COMPLEX(KIND=dp), INTENT(INOUT) :: product(:, :)
116 : COMPLEX(KIND=dp), INTENT(IN) :: overlap(:, :)
117 : REAL(KIND=dp), INTENT(OUT) :: minimum_sv
118 : INTEGER, INTENT(OUT) :: status
119 : REAL(KIND=dp), INTENT(IN) :: sv_tol
120 :
121 89956 : COMPLEX(KIND=dp), ALLOCATABLE :: a(:, :), u(:, :), vh(:, :), work(:)
122 : INTEGER :: n
123 89956 : REAL(KIND=dp), ALLOCATABLE :: rwork(:), sv(:)
124 :
125 89956 : n = SIZE(overlap, 1)
126 89956 : status = -1
127 89956 : minimum_sv = 0.0_dp
128 269868 : IF (n < 1 .OR. SIZE(overlap, 2) /= n .OR. ANY(SHAPE(product) /= [n, n])) RETURN
129 571732 : IF (.NOT. ALL(ieee_is_finite(REAL(overlap, dp)))) RETURN
130 571732 : IF (.NOT. ALL(ieee_is_finite(AIMAG(overlap)))) RETURN
131 1259384 : ALLOCATE (a(n, n), u(n, n), vh(n, n), sv(n), rwork(5*n), work(MAX(1, 4*n)))
132 571732 : a(:, :) = overlap
133 89956 : CALL zgesvd('A', 'A', n, n, a, n, sv, u, n, vh, n, work, SIZE(work), rwork, status)
134 89956 : IF (status /= 0) RETURN
135 254286 : minimum_sv = MINVAL(sv)
136 89956 : IF (minimum_sv <= sv_tol) THEN
137 4 : status = -2
138 4 : RETURN
139 : END IF
140 4262586 : product = MATMUL(product, MATMUL(u, vh))
141 89956 : END SUBROUTINE wilson_step
142 :
143 : ! **************************************************************************************************
144 : !> \brief Sorted Wilson eigenphases in reduced units and total Berry phase in radians.
145 : !> \param product Wilson matrix
146 : !> \param wcc sorted centres in [0,1)
147 : !> \param berry Berry phase in [-pi,pi)
148 : !> \param status zero on success
149 : ! **************************************************************************************************
150 1186 : SUBROUTINE wilson_spectrum(product, wcc, berry, status)
151 : COMPLEX(KIND=dp), INTENT(IN) :: product(:, :)
152 : REAL(KIND=dp), INTENT(OUT) :: wcc(:), berry
153 : INTEGER, INTENT(OUT) :: status
154 :
155 : COMPLEX(KIND=dp) :: dummy(1, 1)
156 1186 : COMPLEX(KIND=dp), ALLOCATABLE :: a(:, :), eig(:), work(:)
157 : INTEGER :: n
158 1186 : REAL(KIND=dp), ALLOCATABLE :: rwork(:)
159 :
160 1186 : n = SIZE(product, 1)
161 1186 : status = -1
162 1186 : berry = 0.0_dp
163 1186 : IF (n < 1 .OR. SIZE(product, 2) /= n .OR. SIZE(wcc) /= n) RETURN
164 11860 : ALLOCATE (a(n, n), eig(n), work(MAX(1, 4*n)), rwork(2*n))
165 7662 : a(:, :) = product
166 1186 : CALL zgeev('N', 'N', n, a, n, eig, dummy, 1, dummy, 1, work, SIZE(work), rwork, status)
167 1186 : IF (status /= 0) RETURN
168 3230 : wcc = MODULO(ATAN2(AIMAG(eig), REAL(eig, dp))/two_pi, 1.0_dp)
169 1186 : CALL sort_wcc(wcc)
170 3230 : berry = two_pi*(MODULO(SUM(wcc) + 0.5_dp, 1.0_dp) - 0.5_dp)
171 1186 : END SUBROUTINE wilson_spectrum
172 :
173 : ! **************************************************************************************************
174 : !> \brief Largest-gap crossing parity for an ordered time-reversal half-surface.
175 : !> Endpoint degeneracy is necessary but does not establish time-reversal symmetry.
176 : !> Sampling convergence and an isolated fixed-rank subspace must be checked by the caller.
177 : !> \param wcc Wilson centres, band index first and surface-line index second
178 : !> \param invariant Z2 parity, or -1 if checks fail
179 : !> \param status 0 success, -1 invalid input, -2 missing Kramers pairs, -3 ambiguous crossing
180 : !> \param pair_tol tolerance on boundary Kramers degeneracy, in reduced units
181 : ! **************************************************************************************************
182 32 : SUBROUTINE z2_from_wcc(wcc, invariant, status, pair_tol)
183 : REAL(KIND=dp), INTENT(IN) :: wcc(:, :)
184 : INTEGER, INTENT(OUT) :: invariant, status
185 : REAL(KIND=dp), INTENT(IN) :: pair_tol
186 :
187 : INTEGER :: crossings, i, j, lines, n
188 : REAL(KIND=dp) :: gap_size, hi, lo
189 32 : REAL(KIND=dp), ALLOCATABLE :: gaps(:), sorted(:, :)
190 :
191 32 : invariant = -1
192 32 : status = -1
193 32 : n = SIZE(wcc, 1)
194 32 : lines = SIZE(wcc, 2)
195 36 : IF (n < 2 .OR. MOD(n, 2) /= 0 .OR. lines < 2) RETURN
196 2444 : IF (.NOT. ALL(ieee_is_finite(wcc))) RETURN
197 192 : ALLOCATE (sorted(n, lines), gaps(lines))
198 2444 : sorted(:, :) = MODULO(wcc, 1.0_dp)
199 804 : DO i = 1, lines
200 772 : CALL sort_wcc(sorted(:, i))
201 804 : CALL largest_gap(sorted(:, i), gaps(i), gap_size)
202 : END DO
203 32 : status = -2
204 32 : IF (.NOT. kramers_pairs(sorted(:, 1), pair_tol)) RETURN
205 30 : IF (.NOT. kramers_pairs(sorted(:, lines), pair_tol)) RETURN
206 30 : crossings = 0
207 766 : DO i = 2, lines
208 736 : lo = MIN(gaps(i - 1), gaps(i))
209 736 : hi = MAX(gaps(i - 1), gaps(i))
210 2310 : DO j = 1, n
211 1544 : IF (ABS(sorted(j, i) - gaps(i - 1)) < crossing_tol) THEN
212 0 : status = -3
213 0 : RETURN
214 : END IF
215 2280 : IF (sorted(j, i) > lo .AND. sorted(j, i) < hi) crossings = crossings + 1
216 : END DO
217 : END DO
218 30 : invariant = MOD(crossings, 2)
219 30 : status = 0
220 36 : END SUBROUTINE z2_from_wcc
221 :
222 : ! **************************************************************************************************
223 : !> \brief Minimum maximal cyclic matching distance between two WCC sets.
224 : !> \param a first WCC set
225 : !> \param b second WCC set
226 : !> \return distance, or huge if sizes differ
227 : ! **************************************************************************************************
228 1758 : FUNCTION wcc_distance(a, b) RESULT(distance)
229 : REAL(KIND=dp), INTENT(IN) :: a(:), b(:)
230 : REAL(KIND=dp) :: distance
231 :
232 : INTEGER :: i, j, n, shift
233 : REAL(KIND=dp) :: delta, error
234 1758 : REAL(KIND=dp), ALLOCATABLE :: aa(:), bb(:)
235 :
236 1758 : distance = HUGE(1.0_dp)
237 1758 : n = SIZE(a)
238 1758 : IF (n /= SIZE(b) .OR. n < 1) RETURN
239 7032 : ALLOCATE (aa(n), bb(n))
240 4566 : aa(:) = MODULO(a, 1.0_dp)
241 4566 : bb(:) = MODULO(b, 1.0_dp)
242 1758 : CALL sort_wcc(aa)
243 1758 : CALL sort_wcc(bb)
244 4566 : DO shift = 0, n - 1
245 : error = 0.0_dp
246 8640 : DO i = 1, n
247 5832 : j = MOD(i - 1 + shift, n) + 1
248 5832 : delta = ABS(aa(i) - bb(j))
249 8640 : error = MAX(error, MIN(delta, 1.0_dp - delta))
250 : END DO
251 4566 : distance = MIN(distance, error)
252 : END DO
253 1758 : END FUNCTION wcc_distance
254 :
255 : ! **************************************************************************************************
256 : !> \brief Test circularly adjacent Kramers pairs (including pairs across the branch cut).
257 : !> \param wcc sorted centres
258 : !> \param tolerance degeneracy tolerance
259 : !> \return whether one of the two cyclic pairings succeeds
260 : ! **************************************************************************************************
261 62 : FUNCTION kramers_pairs(wcc, tolerance) RESULT(paired)
262 : REAL(KIND=dp), INTENT(IN) :: wcc(:), tolerance
263 : LOGICAL :: paired
264 :
265 : INTEGER :: i, j, k, n, offset
266 : LOGICAL :: candidate
267 : REAL(KIND=dp) :: delta
268 :
269 62 : paired = .FALSE.
270 62 : n = SIZE(wcc)
271 186 : DO offset = 0, 1
272 124 : candidate = .TRUE.
273 124 : DO i = 1, n, 2
274 172 : j = MOD(i - 1 + offset, n) + 1
275 172 : k = MOD(i + offset, n) + 1
276 172 : delta = ABS(wcc(j) - wcc(k))
277 208 : candidate = candidate .AND. MIN(delta, 1.0_dp - delta) <= tolerance
278 : END DO
279 186 : paired = paired .OR. candidate
280 : END DO
281 62 : END FUNCTION kramers_pairs
282 :
283 : ! **************************************************************************************************
284 : !> \brief Find the midpoint of the largest gap between sorted centres.
285 : !> \param wcc sorted centres
286 : !> \param centre gap midpoint
287 : !> \param width gap width
288 : ! **************************************************************************************************
289 800 : SUBROUTINE largest_gap(wcc, centre, width)
290 : REAL(KIND=dp), INTENT(IN) :: wcc(:)
291 : REAL(KIND=dp), INTENT(OUT) :: centre, width
292 :
293 : INTEGER :: i, n
294 : REAL(KIND=dp) :: delta
295 :
296 800 : n = SIZE(wcc)
297 800 : width = -1.0_dp
298 2592 : DO i = 1, n
299 1792 : IF (i < n) THEN
300 992 : delta = wcc(i + 1) - wcc(i)
301 : ELSE
302 800 : delta = wcc(1) + 1.0_dp - wcc(n)
303 : END IF
304 2592 : IF (delta > width) THEN
305 1020 : width = delta
306 1020 : centre = MODULO(wcc(i) + delta/2.0_dp, 1.0_dp)
307 : END IF
308 : END DO
309 800 : END SUBROUTINE largest_gap
310 :
311 : ! **************************************************************************************************
312 : !> \brief In-place insertion sort of Wilson centres.
313 : !> \param values values to sort
314 : ! **************************************************************************************************
315 5496 : SUBROUTINE sort_wcc(values)
316 : REAL(KIND=dp), INTENT(INOUT) :: values(:)
317 :
318 : INTEGER :: i, j
319 : REAL(KIND=dp) :: value
320 :
321 9404 : DO i = 2, SIZE(values)
322 3908 : value = values(i)
323 3908 : j = i - 1
324 4772 : DO WHILE (j >= 1)
325 4116 : IF (values(j) <= value) EXIT
326 864 : values(j + 1) = values(j)
327 4116 : j = j - 1
328 : END DO
329 9404 : values(j + 1) = value
330 : END DO
331 5496 : END SUBROUTINE sort_wcc
332 89952 : END MODULE topology_wilson
|
