Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 2 : PROGRAM topology_wilson_unittest
9 : USE iso_fortran_env, ONLY: output_unit
10 2 : 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 2 : io_unit = output_unit
25 2 : CALL test_model(-1.0_dp, 1)
26 2 : CALL test_model(-3.0_dp, 0)
27 2 : CALL test_model(1.0_dp, 1)
28 2 : CALL test_model(3.0_dp, 0)
29 2 : CALL test_checks()
30 2 : CALL test_chern_model(-1.0_dp, 1)
31 2 : CALL test_chern_model(1.0_dp, -1)
32 2 : CALL test_chern_model(-3.0_dp, 0)
33 2 : CALL test_chern_checks()
34 2 : WRITE (io_unit, *) "Wilson-loop, Z2 and Chern unit tests passed."
35 :
36 : CONTAINS
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 6 : 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 396 : DO line = 1, nline
57 390 : ky = 2.0_dp*pi*REAL(line - 1, dp)/REAL(nline - 1, dp)
58 16380 : DO i = 1, npoint
59 15990 : kx = 2.0_dp*pi*REAL(i - 1, dp)/REAL(npoint, dp)
60 15990 : d = mass + COS(kx) + COS(ky)
61 15990 : h(1, 1) = d
62 15990 : h(2, 2) = -d
63 15990 : h(1, 2) = CMPLX(SIN(kx), -SIN(ky), dp)
64 15990 : h(2, 1) = CONJG(h(1, 2))
65 15990 : CALL zheev('V', 'U', 2, h, 2, eigenvalues, work, SIZE(work), rwork, status)
66 15990 : IF (status /= 0) ERROR STOP "Chern model diagonalization failed"
67 48360 : states(:, i) = h(:, 1)
68 : END DO
69 1170 : PRODUCT(:, :) = 1.0_dp
70 16380 : DO i = 1, npoint
71 15990 : j = MOD(i, npoint) + 1
72 47970 : overlap(1, 1) = DOT_PRODUCT(states(:, i), states(:, j))
73 15990 : CALL wilson_step(product, overlap, minimum_sv, status, tolerance)
74 16380 : IF (status /= 0) ERROR STOP "Chern model overlap is singular"
75 : END DO
76 390 : CALL wilson_spectrum(product, wcc(:, line), berry, status)
77 396 : IF (status /= 0) ERROR STOP "Chern model spectrum failed"
78 : END DO
79 6 : CALL chern_from_wcc(wcc, invariant, winding, status, tolerance)
80 6 : IF (status /= 0 .OR. invariant /= expected) ERROR STOP "Incorrect model Chern number"
81 6 : IF (ABS(winding - REAL(expected, dp)) > tolerance) ERROR STOP "Incorrect Chern winding"
82 6 : CALL chern_from_wcc(wcc(:, nline:1:-1), invariant, winding, status, tolerance)
83 6 : IF (status /= 0 .OR. invariant /= -expected) ERROR STOP "Chern orientation reversal failed"
84 396 : doubled(1, :) = wcc(1, :)
85 396 : doubled(2, :) = wcc(1, :)
86 6 : CALL chern_from_wcc(doubled, invariant, winding, status, tolerance)
87 6 : IF (status /= 0 .OR. invariant /= 2*expected) ERROR STOP "Direct-sum Chern additivity failed"
88 2 : END SUBROUTINE test_chern_model
89 :
90 : ! **************************************************************************************************
91 : !> \brief Reject unclosed surfaces and unresolved determinant-phase steps.
92 : ! **************************************************************************************************
93 2 : SUBROUTINE test_chern_checks()
94 : INTEGER :: invariant, status
95 : REAL(KIND=dp) :: surface(1, 5), winding
96 :
97 2 : surface(1, :) = [0.0_dp, 0.1_dp, 0.2_dp, 0.3_dp, 0.4_dp]
98 2 : CALL chern_from_wcc(surface, invariant, winding, status, tolerance)
99 2 : IF (status /= -2) ERROR STOP "Unclosed Chern surface was accepted"
100 2 : surface(1, :) = [0.0_dp, 0.3_dp, 0.6_dp, 0.9_dp, 1.0_dp]
101 2 : CALL chern_from_wcc(surface, invariant, winding, status, tolerance)
102 2 : IF (status /= -3) ERROR STOP "Unresolved Chern phase step was accepted"
103 22 : surface(:, :) = 0.1_dp
104 2 : CALL chern_from_wcc(surface, invariant, winding, status, tolerance)
105 2 : IF (status /= 0 .OR. invariant /= 0) ERROR STOP "Constant Chern surface was rejected"
106 2 : 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 8 : 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 256 : DO line = 1, nline
127 248 : ky = pi*REAL(line - 1, dp)/REAL(nline - 1, dp)
128 10416 : DO i = 1, npoint
129 10168 : kx = 2.0_dp*pi*REAL(i - 1, dp)/REAL(npoint, dp)
130 10168 : CALL occupied_states(kx, ky, mass, states(:, :, i))
131 10416 : CALL gauge_matrix(kx, ky, gauges(:, :, i))
132 : END DO
133 248 : PRODUCT(:, :) = 0.0_dp
134 248 : PRODUCT(1, 1) = 1.0_dp
135 248 : PRODUCT(2, 2) = 1.0_dp
136 248 : gauged(:, :) = product
137 248 : reversed(:, :) = product
138 10416 : DO i = 1, npoint
139 10168 : j = MOD(i, npoint) + 1
140 305040 : links(:, :, i) = MATMUL(CONJG(TRANSPOSE(states(:, :, i))), states(:, :, j))
141 10168 : CALL wilson_step(product, links(:, :, i), minimum_sv, status, tolerance)
142 10168 : IF (status /= 0) ERROR STOP "BHZ overlap is singular"
143 : overlap(:, :) = MATMUL(CONJG(TRANSPOSE(gauges(:, :, i))), &
144 294872 : MATMUL(links(:, :, i), gauges(:, :, j)))
145 10168 : CALL wilson_step(gauged, overlap, minimum_sv, status, tolerance)
146 10416 : IF (status /= 0) ERROR STOP "Gauge-transformed BHZ overlap is singular"
147 : END DO
148 10416 : DO i = npoint, 1, -1
149 71176 : overlap(:, :) = CONJG(TRANSPOSE(links(:, :, i)))
150 10168 : CALL wilson_step(reversed, overlap, minimum_sv, status, tolerance)
151 10416 : IF (status /= 0) ERROR STOP "Reversed BHZ overlap is singular"
152 : END DO
153 248 : CALL wilson_spectrum(product, wcc(:, line), berry, status)
154 248 : IF (status /= 0) ERROR STOP "BHZ spectrum failed"
155 248 : CALL wilson_spectrum(gauged, gauge_wcc(:, line), berry, status)
156 248 : IF (status /= 0) ERROR STOP "Gauge-transformed spectrum failed"
157 248 : CALL wilson_spectrum(reversed, reverse_wcc(:, line), berry, status)
158 248 : IF (status /= 0) ERROR STOP "Reversed spectrum failed"
159 248 : IF (wcc_distance(wcc(:, line), gauge_wcc(:, line)) > tolerance) THEN
160 0 : ERROR STOP "Wilson spectrum depends on occupied-space gauge"
161 : END IF
162 752 : IF (wcc_distance(wcc(:, line), -reverse_wcc(:, line)) > tolerance) THEN
163 0 : ERROR STOP "Loop reversal did not conjugate Wilson eigenvalues"
164 : END IF
165 : END DO
166 8 : CALL z2_from_wcc(wcc, invariant, status, tolerance)
167 8 : IF (status /= 0 .OR. invariant /= expected) ERROR STOP "Incorrect BHZ Z2 invariant"
168 8 : CALL z2_from_wcc(gauge_wcc, invariant, status, tolerance)
169 8 : IF (status /= 0 .OR. invariant /= expected) ERROR STOP "Gauge-dependent Z2 invariant"
170 8 : CALL z2_from_wcc(reverse_wcc, invariant, status, tolerance)
171 8 : IF (status /= 0 .OR. invariant /= expected) ERROR STOP "Reversal-dependent Z2 invariant"
172 8 : 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 10168 : 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 10168 : d = mass + COS(kx) + COS(ky)
190 10168 : h(:, :) = 0.0_dp
191 10168 : h(1, 1) = d
192 10168 : h(2, 2) = -d
193 10168 : h(1, 2) = CMPLX(SIN(kx), -SIN(ky), dp)
194 10168 : h(2, 1) = CONJG(h(1, 2))
195 10168 : h(3, 3) = d
196 10168 : h(4, 4) = -d
197 10168 : h(3, 4) = CMPLX(-SIN(kx), -SIN(ky), dp)
198 10168 : h(4, 3) = CONJG(h(3, 4))
199 10168 : CALL zheev('V', 'U', 4, h, 4, eigenvalues, work, SIZE(work), rwork, status)
200 10168 : IF (status /= 0) ERROR STOP "BHZ diagonalization failed"
201 111848 : states(:, :) = h(:, 1:2)
202 10168 : 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 10168 : 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 10168 : angle = 0.37_dp + 0.61_dp*SIN(3.0_dp*kx + ky)
218 10168 : phase = EXP(CMPLX(0.0_dp, 0.73_dp*COS(kx - 2.0_dp*ky), dp))
219 10168 : determinant_phase = EXP(CMPLX(0.0_dp, 0.29_dp*SIN(kx + ky), dp))
220 10168 : matrix(1, 1) = COS(angle)
221 10168 : matrix(1, 2) = SIN(angle)*phase
222 10168 : matrix(2, 1) = -SIN(angle)*CONJG(phase)
223 10168 : matrix(2, 2) = COS(angle)
224 30504 : matrix(:, 1) = matrix(:, 1)*determinant_phase
225 10168 : END SUBROUTINE gauge_matrix
226 :
227 : ! **************************************************************************************************
228 : !> \brief Analytic polar/Berry phases and singular-link, Kramers-pair and surface-resolution checks.
229 : ! **************************************************************************************************
230 2 : 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 2 : PRODUCT(:, :) = 0.0_dp
236 2 : PRODUCT(1, 1) = 1.0_dp
237 2 : PRODUCT(2, 2) = 1.0_dp
238 2 : overlap(:, :) = 0.0_dp
239 2 : CALL wilson_step(product, overlap, minimum_sv, status, tolerance)
240 2 : IF (status /= -2) ERROR STOP "Singular overlap was accepted"
241 :
242 2 : overlap(1, 1) = 0.8_dp*EXP(CMPLX(0.0_dp, 0.2_dp*pi, dp))
243 2 : overlap(2, 2) = 0.6_dp*EXP(CMPLX(0.0_dp, 0.6_dp*pi, dp))
244 2 : CALL wilson_step(product, overlap, minimum_sv, status, tolerance)
245 2 : IF (status /= 0 .OR. ABS(minimum_sv - 0.6_dp) > tolerance) THEN
246 0 : ERROR STOP "Polar factor singular value is incorrect"
247 : END IF
248 2 : CALL wilson_spectrum(product, wcc, berry, status)
249 2 : IF (status /= 0) ERROR STOP "Analytic spectrum failed"
250 2 : IF (wcc_distance(wcc, [0.1_dp, 0.3_dp]) > tolerance) ERROR STOP "Incorrect analytic WCC"
251 2 : IF (ABS(berry - 0.8_dp*pi) > tolerance) ERROR STOP "Incorrect analytic Berry phase"
252 :
253 6 : surface(:, 1) = wcc
254 6 : surface(:, 2) = wcc
255 6 : surface(:, 3) = wcc
256 2 : CALL z2_from_wcc(surface, invariant, status, tolerance)
257 2 : IF (status /= -2 .OR. invariant /= -1) ERROR STOP "Missing Kramers pairs were accepted"
258 20 : surface(:, :) = 0.1_dp
259 6 : surface(:, 2) = 0.6_dp
260 2 : IF (surface_resolved(surface)) ERROR STOP "Unresolved surface was accepted"
261 20 : surface(:, :) = 0.1_dp
262 2 : IF (.NOT. surface_resolved(surface)) ERROR STOP "Constant surface was rejected"
263 2 : CALL z2_from_wcc(surface, invariant, status, tolerance)
264 2 : IF (status /= 0 .OR. invariant /= 0) ERROR STOP "Incorrect constant-surface parity"
265 2 : IF (wcc_distance([0.99_dp, 0.2_dp], [0.2_dp, -0.01_dp]) > tolerance) THEN
266 0 : ERROR STOP "Circular WCC matching failed"
267 : END IF
268 2 : END SUBROUTINE test_checks
269 :
270 : END PROGRAM topology_wilson_unittest
|
