(git:f2099e5)
Loading...
Searching...
No Matches
topology_wilson_unittest.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8PROGRAM topology_wilson_unittest
9 USE iso_fortran_env, ONLY: output_unit
10 USE kinds, ONLY: dp
11 USE topology_wilson, ONLY: chern_from_wcc,&
12 surface_resolved,&
13 wcc_distance,&
14 wilson_spectrum,&
15 wilson_step,&
16 z2_from_wcc
17
18 IMPLICIT NONE
19
20 REAL(kind=dp), PARAMETER :: pi = 3.1415926535897932384626433832795_dp, tolerance = 1.e-10_dp
21
22 INTEGER :: io_unit
23
24 io_unit = output_unit
25 CALL test_model(-1.0_dp, 1)
26 CALL test_model(-3.0_dp, 0)
27 CALL test_model(1.0_dp, 1)
28 CALL test_model(3.0_dp, 0)
29 CALL test_checks()
30 CALL test_chern_model(-1.0_dp, 1)
31 CALL test_chern_model(1.0_dp, -1)
32 CALL test_chern_model(-3.0_dp, 0)
33 CALL test_chern_checks()
34 WRITE (io_unit, *) "Wilson-loop, Z2 and Chern unit tests passed."
35
36CONTAINS
37
38! **************************************************************************************************
39!> \brief Known first Chern numbers of the occupied two-band lattice Dirac model.
40!> \param mass Dirac mass in h = (sin kx, sin ky, m+cos kx+cos ky).sigma
41!> \param expected Chern number in the native/Z2Pack Wilson-winding convention
42! **************************************************************************************************
43 SUBROUTINE test_chern_model(mass, expected)
44 REAL(KIND=dp), INTENT(IN) :: mass
45 INTEGER, INTENT(IN) :: expected
46
47 INTEGER, PARAMETER :: nline = 65, npoint = 41
48
49 COMPLEX(KIND=dp) :: h(2, 2), overlap(1, 1), PRODUCT(1, 1), &
50 states(2, npoint), work(8)
51 INTEGER :: i, invariant, j, line, status
52 REAL(KIND=dp) :: berry, d, doubled(2, nline), &
53 eigenvalues(2), kx, ky, minimum_sv, &
54 rwork(4), wcc(1, nline), winding
55
56 DO line = 1, nline
57 ky = 2.0_dp*pi*real(line - 1, dp)/real(nline - 1, dp)
58 DO i = 1, npoint
59 kx = 2.0_dp*pi*real(i - 1, dp)/real(npoint, dp)
60 d = mass + cos(kx) + cos(ky)
61 h(1, 1) = d
62 h(2, 2) = -d
63 h(1, 2) = cmplx(sin(kx), -sin(ky), dp)
64 h(2, 1) = conjg(h(1, 2))
65 CALL zheev('V', 'U', 2, h, 2, eigenvalues, work, SIZE(work), rwork, status)
66 IF (status /= 0) error stop "Chern model diagonalization failed"
67 states(:, i) = h(:, 1)
68 END DO
69 product(:, :) = 1.0_dp
70 DO i = 1, npoint
71 j = mod(i, npoint) + 1
72 overlap(1, 1) = dot_product(states(:, i), states(:, j))
73 CALL wilson_step(product, overlap, minimum_sv, status, tolerance)
74 IF (status /= 0) error stop "Chern model overlap is singular"
75 END DO
76 CALL wilson_spectrum(product, wcc(:, line), berry, status)
77 IF (status /= 0) error stop "Chern model spectrum failed"
78 END DO
79 CALL chern_from_wcc(wcc, invariant, winding, status, tolerance)
80 IF (status /= 0 .OR. invariant /= expected) error stop "Incorrect model Chern number"
81 IF (abs(winding - real(expected, dp)) > tolerance) error stop "Incorrect Chern winding"
82 CALL chern_from_wcc(wcc(:, nline:1:-1), invariant, winding, status, tolerance)
83 IF (status /= 0 .OR. invariant /= -expected) error stop "Chern orientation reversal failed"
84 doubled(1, :) = wcc(1, :)
85 doubled(2, :) = wcc(1, :)
86 CALL chern_from_wcc(doubled, invariant, winding, status, tolerance)
87 IF (status /= 0 .OR. invariant /= 2*expected) error stop "Direct-sum Chern additivity failed"
88 END SUBROUTINE test_chern_model
89
90! **************************************************************************************************
91!> \brief Reject unclosed surfaces and unresolved determinant-phase steps.
92! **************************************************************************************************
93 SUBROUTINE test_chern_checks()
94 INTEGER :: invariant, status
95 REAL(KIND=dp) :: surface(1, 5), winding
96
97 surface(1, :) = [0.0_dp, 0.1_dp, 0.2_dp, 0.3_dp, 0.4_dp]
98 CALL chern_from_wcc(surface, invariant, winding, status, tolerance)
99 IF (status /= -2) error stop "Unclosed Chern surface was accepted"
100 surface(1, :) = [0.0_dp, 0.3_dp, 0.6_dp, 0.9_dp, 1.0_dp]
101 CALL chern_from_wcc(surface, invariant, winding, status, tolerance)
102 IF (status /= -3) error stop "Unresolved Chern phase step was accepted"
103 surface(:, :) = 0.1_dp
104 CALL chern_from_wcc(surface, invariant, winding, status, tolerance)
105 IF (status /= 0 .OR. invariant /= 0) error stop "Constant Chern surface was rejected"
106 END SUBROUTINE test_chern_checks
107
108! **************************************************************************************************
109!> \brief Known trivial/nontrivial lattice BHZ models, non-Abelian gauge changes and loop reversal.
110!> \param mass Dirac mass; |mass| < 2 (excluding zero) is topological, |mass| > 2 is trivial
111!> \param expected expected Z2 invariant
112! **************************************************************************************************
113 SUBROUTINE test_model(mass, expected)
114 REAL(KIND=dp), INTENT(IN) :: mass
115 INTEGER, INTENT(IN) :: expected
116
117 INTEGER, PARAMETER :: nline = 31, npoint = 41
118
119 COMPLEX(KIND=dp) :: gauged(2, 2), gauges(2, 2, npoint), links(2, 2, npoint), overlap(2, 2), &
120 PRODUCT(2, 2), reversed(2, 2), states(4, 2, npoint)
121 INTEGER :: i, invariant, j, line, status
122 REAL(KIND=dp) :: berry, gauge_wcc(2, nline), kx, ky, &
123 minimum_sv, reverse_wcc(2, nline), &
124 wcc(2, nline)
125
126 DO line = 1, nline
127 ky = pi*real(line - 1, dp)/real(nline - 1, dp)
128 DO i = 1, npoint
129 kx = 2.0_dp*pi*real(i - 1, dp)/real(npoint, dp)
130 CALL occupied_states(kx, ky, mass, states(:, :, i))
131 CALL gauge_matrix(kx, ky, gauges(:, :, i))
132 END DO
133 product(:, :) = 0.0_dp
134 product(1, 1) = 1.0_dp
135 product(2, 2) = 1.0_dp
136 gauged(:, :) = product
137 reversed(:, :) = product
138 DO i = 1, npoint
139 j = mod(i, npoint) + 1
140 links(:, :, i) = matmul(conjg(transpose(states(:, :, i))), states(:, :, j))
141 CALL wilson_step(product, links(:, :, i), minimum_sv, status, tolerance)
142 IF (status /= 0) error stop "BHZ overlap is singular"
143 overlap(:, :) = matmul(conjg(transpose(gauges(:, :, i))), &
144 matmul(links(:, :, i), gauges(:, :, j)))
145 CALL wilson_step(gauged, overlap, minimum_sv, status, tolerance)
146 IF (status /= 0) error stop "Gauge-transformed BHZ overlap is singular"
147 END DO
148 DO i = npoint, 1, -1
149 overlap(:, :) = conjg(transpose(links(:, :, i)))
150 CALL wilson_step(reversed, overlap, minimum_sv, status, tolerance)
151 IF (status /= 0) error stop "Reversed BHZ overlap is singular"
152 END DO
153 CALL wilson_spectrum(product, wcc(:, line), berry, status)
154 IF (status /= 0) error stop "BHZ spectrum failed"
155 CALL wilson_spectrum(gauged, gauge_wcc(:, line), berry, status)
156 IF (status /= 0) error stop "Gauge-transformed spectrum failed"
157 CALL wilson_spectrum(reversed, reverse_wcc(:, line), berry, status)
158 IF (status /= 0) error stop "Reversed spectrum failed"
159 IF (wcc_distance(wcc(:, line), gauge_wcc(:, line)) > tolerance) THEN
160 error stop "Wilson spectrum depends on occupied-space gauge"
161 END IF
162 IF (wcc_distance(wcc(:, line), -reverse_wcc(:, line)) > tolerance) THEN
163 error stop "Loop reversal did not conjugate Wilson eigenvalues"
164 END IF
165 END DO
166 CALL z2_from_wcc(wcc, invariant, status, tolerance)
167 IF (status /= 0 .OR. invariant /= expected) error stop "Incorrect BHZ Z2 invariant"
168 CALL z2_from_wcc(gauge_wcc, invariant, status, tolerance)
169 IF (status /= 0 .OR. invariant /= expected) error stop "Gauge-dependent Z2 invariant"
170 CALL z2_from_wcc(reverse_wcc, invariant, status, tolerance)
171 IF (status /= 0 .OR. invariant /= expected) error stop "Reversal-dependent Z2 invariant"
172 END SUBROUTINE test_model
173
174! **************************************************************************************************
175!> \brief Occupied eigenvectors of diag(h(k), h*(-k)), with h = (sin kx, sin ky, m+cos kx+cos ky).sigma.
176!> \param kx first reduced angle
177!> \param ky second reduced angle
178!> \param mass Dirac mass
179!> \param states two occupied eigenvectors
180! **************************************************************************************************
181 SUBROUTINE occupied_states(kx, ky, mass, states)
182 REAL(KIND=dp), INTENT(IN) :: kx, ky, mass
183 COMPLEX(KIND=dp), INTENT(OUT) :: states(4, 2)
184
185 COMPLEX(KIND=dp) :: h(4, 4), work(16)
186 INTEGER :: status
187 REAL(KIND=dp) :: d, eigenvalues(4), rwork(10)
188
189 d = mass + cos(kx) + cos(ky)
190 h(:, :) = 0.0_dp
191 h(1, 1) = d
192 h(2, 2) = -d
193 h(1, 2) = cmplx(sin(kx), -sin(ky), dp)
194 h(2, 1) = conjg(h(1, 2))
195 h(3, 3) = d
196 h(4, 4) = -d
197 h(3, 4) = cmplx(-sin(kx), -sin(ky), dp)
198 h(4, 3) = conjg(h(3, 4))
199 CALL zheev('V', 'U', 4, h, 4, eigenvalues, work, SIZE(work), rwork, status)
200 IF (status /= 0) error stop "BHZ diagonalization failed"
201 states(:, :) = h(:, 1:2)
202 END SUBROUTINE occupied_states
203
204! **************************************************************************************************
205!> \brief Deterministic k-dependent U(2) rotations, mixing the two occupied states.
206!> \param kx first reduced angle
207!> \param ky second reduced angle
208!> \param matrix gauge matrix
209! **************************************************************************************************
210 SUBROUTINE gauge_matrix(kx, ky, matrix)
211 REAL(KIND=dp), INTENT(IN) :: kx, ky
212 COMPLEX(KIND=dp), INTENT(OUT) :: matrix(2, 2)
213
214 COMPLEX(KIND=dp) :: determinant_phase, phase
215 REAL(KIND=dp) :: angle
216
217 angle = 0.37_dp + 0.61_dp*sin(3.0_dp*kx + ky)
218 phase = exp(cmplx(0.0_dp, 0.73_dp*cos(kx - 2.0_dp*ky), dp))
219 determinant_phase = exp(cmplx(0.0_dp, 0.29_dp*sin(kx + ky), dp))
220 matrix(1, 1) = cos(angle)
221 matrix(1, 2) = sin(angle)*phase
222 matrix(2, 1) = -sin(angle)*conjg(phase)
223 matrix(2, 2) = cos(angle)
224 matrix(:, 1) = matrix(:, 1)*determinant_phase
225 END SUBROUTINE gauge_matrix
226
227! **************************************************************************************************
228!> \brief Analytic polar/Berry phases and singular-link, Kramers-pair and surface-resolution checks.
229! **************************************************************************************************
230 SUBROUTINE test_checks()
231 COMPLEX(KIND=dp) :: overlap(2, 2), PRODUCT(2, 2)
232 INTEGER :: invariant, status
233 REAL(KIND=dp) :: berry, minimum_sv, surface(2, 3), wcc(2)
234
235 product(:, :) = 0.0_dp
236 product(1, 1) = 1.0_dp
237 product(2, 2) = 1.0_dp
238 overlap(:, :) = 0.0_dp
239 CALL wilson_step(product, overlap, minimum_sv, status, tolerance)
240 IF (status /= -2) error stop "Singular overlap was accepted"
241
242 overlap(1, 1) = 0.8_dp*exp(cmplx(0.0_dp, 0.2_dp*pi, dp))
243 overlap(2, 2) = 0.6_dp*exp(cmplx(0.0_dp, 0.6_dp*pi, dp))
244 CALL wilson_step(product, overlap, minimum_sv, status, tolerance)
245 IF (status /= 0 .OR. abs(minimum_sv - 0.6_dp) > tolerance) THEN
246 error stop "Polar factor singular value is incorrect"
247 END IF
248 CALL wilson_spectrum(product, wcc, berry, status)
249 IF (status /= 0) error stop "Analytic spectrum failed"
250 IF (wcc_distance(wcc, [0.1_dp, 0.3_dp]) > tolerance) error stop "Incorrect analytic WCC"
251 IF (abs(berry - 0.8_dp*pi) > tolerance) error stop "Incorrect analytic Berry phase"
252
253 surface(:, 1) = wcc
254 surface(:, 2) = wcc
255 surface(:, 3) = wcc
256 CALL z2_from_wcc(surface, invariant, status, tolerance)
257 IF (status /= -2 .OR. invariant /= -1) error stop "Missing Kramers pairs were accepted"
258 surface(:, :) = 0.1_dp
259 surface(:, 2) = 0.6_dp
260 IF (surface_resolved(surface)) error stop "Unresolved surface was accepted"
261 surface(:, :) = 0.1_dp
262 IF (.NOT. surface_resolved(surface)) error stop "Constant surface was rejected"
263 CALL z2_from_wcc(surface, invariant, status, tolerance)
264 IF (status /= 0 .OR. invariant /= 0) error stop "Incorrect constant-surface parity"
265 IF (wcc_distance([0.99_dp, 0.2_dp], [0.2_dp, -0.01_dp]) > tolerance) THEN
266 error stop "Circular WCC matching failed"
267 END IF
268 END SUBROUTINE test_checks
269
270END PROGRAM topology_wilson_unittest
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
topology_wilson
Local Wilson-loop linear algebra, independent of the electronic-structure representation....
Definition topology_wilson.F:12
topology_wilson::wcc_distance
real(kind=dp) function, public wcc_distance(a, b)
Minimum maximal cyclic matching distance between two WCC sets.
Definition topology_wilson.F:229
topology_wilson::wilson_step
subroutine, public wilson_step(product, overlap, minimum_sv, status, sv_tol)
Multiply by the unitary polar factor of an overlap. Reject rank-deficient links.
Definition topology_wilson.F:115
topology_wilson::z2_from_wcc
subroutine, public z2_from_wcc(wcc, invariant, status, pair_tol)
Largest-gap crossing parity for an ordered time-reversal half-surface. Endpoint degeneracy is necessa...
Definition topology_wilson.F:183
topology_wilson::surface_resolved
logical function, public surface_resolved(wcc)
Conservative neighbouring-line movement and largest-gap separation checks.
Definition topology_wilson.F:76
topology_wilson::wilson_spectrum
subroutine, public wilson_spectrum(product, wcc, berry, status)
Sorted Wilson eigenphases in reduced units and total Berry phase in radians.
Definition topology_wilson.F:151
topology_wilson::chern_from_wcc
subroutine, public chern_from_wcc(wcc, invariant, winding, status, closure_tol)
First Chern number from determinant Wilson-phase winding on a closed surface. Both transverse endpoin...
Definition topology_wilson.F:37
test_chern_model
subroutine test_chern_model(mass, expected)
Known first Chern numbers of the occupied two-band lattice Dirac model.
Definition topology_wilson_unittest.F:44
topology_wilson_unittest
program topology_wilson_unittest
Definition topology_wilson_unittest.F:8