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 : MODULE kubo_test_cases
8 : IMPLICIT NONE
9 : PRIVATE
10 : PUBLIC :: check_projected_current, check_bloch_current, check_hall_response
11 : CONTAINS
12 :
13 : ! **************************************************************************************************
14 : !> \brief Validate projected currents and their AO-metric covariance.
15 : ! **************************************************************************************************
16 2 : SUBROUTINE check_projected_current()
17 : USE kinds, ONLY: dp
18 : USE kubo_projected, ONLY: projected_current
19 : COMPLEX(dp) :: b(4, 4), h(4, 4), h2(4, 4), j(4, 4, 3), j2(4, 4, 3), norm(3, 3), norm2(3, 3), &
20 : s(4, 4), s2(4, 4), x(4, 4, 3), x2(4, 4, 3)
21 : INTEGER :: a, c, i, status
22 : REAL(dp) :: e(4), e2(4)
23 :
24 2 : h = 0.0_dp
25 2 : s = 0.0_dp
26 2 : x = 0.0_dp
27 10 : DO i = 1, 4
28 10 : s(i, i) = 1.0_dp
29 : END DO
30 : ! Kramers-degenerate spinor model with imaginary SOC matrix elements.
31 2 : h(2, 2) = 0.8_dp
32 2 : h(4, 4) = 0.8_dp
33 2 : h(1, 2) = CMPLX(0.0_dp, 0.3_dp, dp)
34 2 : h(2, 1) = CONJG(h(1, 2))
35 2 : h(3, 4) = CONJG(h(1, 2))
36 2 : h(4, 3) = h(1, 2)
37 2 : x(1, 2, 1) = 1.0_dp
38 2 : x(2, 1, 1) = 1.0_dp
39 2 : x(3, 4, 1) = 1.0_dp
40 2 : x(4, 3, 1) = 1.0_dp
41 2 : x(1, 1, 2) = -0.5_dp
42 2 : x(2, 2, 2) = 0.5_dp
43 2 : x(3, 3, 2) = -0.5_dp
44 2 : x(4, 4, 2) = 0.5_dp
45 42 : x(:, :, 3) = 2*x(:, :, 1) - x(:, :, 2)
46 2 : CALL projected_current(h, s, x, e, j, status)
47 2 : IF (status /= 0) ERROR STOP 'Spinor current failed'
48 10 : IF (MAXVAL(ABS(e - [-0.1_dp, -0.1_dp, 0.9_dp, 0.9_dp])) > 1.e-12_dp) ERROR STOP 'Spinor spectrum failed'
49 8 : DO a = 1, 3
50 126 : IF (MAXVAL(ABS(j(:, :, a) + CONJG(TRANSPOSE(j(:, :, a))))) > 1.e-12_dp) ERROR STOP 'Current is not anti-Hermitian'
51 26 : DO c = 1, 3
52 384 : norm(a, c) = SUM(CONJG(j(:, :, a))*j(:, :, c))
53 : END DO
54 : END DO
55 2 : IF (ABS(norm(1, 1) - 4.0_dp) > 1.e-12_dp) ERROR STOP 'Analytic current strength failed'
56 : ! Complex nonorthogonal basis mixing spin and orbitals, including a degenerate subspace.
57 2 : b = s
58 2 : b(1, 1) = 1.7_dp
59 2 : b(2, 2) = 0.8_dp
60 2 : b(3, 1) = CMPLX(0.2_dp, 0.4_dp, dp)
61 2 : b(1, 4) = CMPLX(-0.3_dp, 0.1_dp, dp)
62 2 : b(2, 3) = 0.4_dp
63 338 : h2 = MATMUL(CONJG(TRANSPOSE(b)), MATMUL(h, b))
64 170 : s2 = MATMUL(CONJG(TRANSPOSE(b)), b)
65 8 : DO a = 1, 3
66 1256 : x2(:, :, a) = MATMUL(CONJG(TRANSPOSE(b)), MATMUL(x(:, :, a), b)) + REAL(a, dp)*s2
67 : END DO
68 2 : CALL projected_current(h2, s2, x2, e2, j2, status)
69 10 : IF (status /= 0 .OR. MAXVAL(ABS(e - e2)) > 1.e-12_dp) ERROR STOP 'Nonorthogonal current spectrum failed'
70 8 : DO a = 1, 3
71 26 : DO c = 1, 3
72 384 : norm2(a, c) = SUM(CONJG(j2(:, :, a))*j2(:, :, c))
73 : END DO
74 : END DO
75 26 : IF (MAXVAL(ABS(norm - norm2)) > 1.e-11_dp) ERROR STOP 'Basis/origin/gauge invariance failed'
76 2 : s(1, 1) = -1.0_dp
77 2 : CALL projected_current(h, s, x, e, j, status)
78 2 : IF (status /= -2) ERROR STOP 'Indefinite overlap accepted'
79 2 : s(1, 1) = 1.0_dp
80 2 : x(1, 2, 1) = 2.0_dp
81 2 : CALL projected_current(h, s, x, e, j, status)
82 2 : IF (status /= -1) ERROR STOP 'Non-Hermitian position accepted'
83 2 : PRINT *, 'Projected spinor current, metric and gauge tests passed.'
84 2 : END SUBROUTINE check_projected_current
85 :
86 : ! **************************************************************************************************
87 : !> \brief Validate Bloch currents, degenerate blocks and k-dependent AO gauges.
88 : ! **************************************************************************************************
89 2 : SUBROUTINE check_bloch_current()
90 2 : USE kinds, ONLY: dp
91 : USE kubo_bloch, ONLY: bloch_current, &
92 : bloch_current_from_frame, &
93 : fermi_slope
94 : USE kubo_projected, ONLY: projected_current
95 : COMPLEX(dp), PARAMETER :: iu = (0.0_dp, 1.0_dp)
96 :
97 : COMPLEX(dp) :: b(2, 2), d(2, 2, 3), d2(2, 2, 3), db(2, 2), dh(2, 2, 3), dh2(2, 2, 3), &
98 : ds(2, 2, 3), ds2(2, 2, 3), expected(2, 2), frame(2, 2), h(2, 2), h2(2, 2), j(2, 2, 3), &
99 : j2(2, 2, 3), s(2, 2), s2(2, 2)
100 : INTEGER :: i, ik, status
101 : REAL(dp) :: e(2), e2(2), k, residuals(2), root, &
102 : shift, speed, strength, strength2
103 :
104 2 : s = 0.0_dp
105 2 : s(1, 1) = 1.0_dp
106 2 : s(2, 2) = 1.0_dp
107 32 : DO ik = -7, 7
108 30 : k = REAL(ik, dp)/5
109 210 : h = k*k*s
110 30 : h(1, 1) = h(1, 1) + 0.3_dp
111 30 : h(2, 2) = h(2, 2) - 0.3_dp
112 30 : h(1, 2) = 0.4_dp*SIN(k) - iu*0.2_dp*COS(k)
113 30 : h(2, 1) = CONJG(h(1, 2))
114 30 : dh = 0.0_dp
115 210 : dh(:, :, 1) = 2*k*s
116 30 : dh(1, 2, 1) = 0.4_dp*COS(k) + iu*0.2_dp*SIN(k)
117 30 : dh(2, 1, 1) = CONJG(dh(1, 2, 1))
118 30 : ds = 0.0_dp
119 30 : d = 0.0_dp
120 30 : CALL bloch_current(h, s, dh, ds, d, e, j, status, frame)
121 30 : IF (status /= 0) ERROR STOP 'Bloch current failed'
122 30 : CALL bloch_current_from_frame(h, s, dh, ds, d, e, frame, j2, 1.e-9_dp, residuals, status)
123 660 : IF (status /= 0 .OR. MAXVAL(ABS(j2 - j)) > 1.e-12_dp) ERROR STOP 'Reused eigenframe failed'
124 90 : CALL bloch_current_from_frame(h, s, dh, ds, d, e + 0.1_dp, frame, j2, 1.e-9_dp, residuals, status)
125 30 : IF (status /= -4) ERROR STOP 'Wrong frame eigenvalues accepted'
126 210 : CALL bloch_current_from_frame(h, s, dh, ds, d, e, 1.1_dp*frame, j2, 1.e-9_dp, residuals, status)
127 30 : IF (status /= -4) ERROR STOP 'Nonorthonormal frame accepted'
128 30 : root = SQRT(0.09_dp + 0.16_dp*SIN(k)**2 + 0.04_dp*COS(k)**2)
129 90 : IF (MAXVAL(ABS(e - [k*k - root, k*k + root])) > 1.e-12_dp) ERROR STOP 'Analytic bands failed'
130 30 : speed = 0.12_dp*SIN(k)*COS(k)/root
131 30 : IF (ABS(iu*j(1, 1, 1) - (2*k - speed)) > 1.e-12_dp) ERROR STOP 'Lower band velocity failed'
132 30 : IF (ABS(iu*j(2, 2, 1) - (2*k + speed)) > 1.e-12_dp) ERROR STOP 'Upper band velocity failed'
133 210 : IF (MAXVAL(ABS(j(:, :, 1) + CONJG(TRANSPOSE(j(:, :, 1))))) > 1.e-12_dp) THEN
134 0 : ERROR STOP 'Bloch current not anti-Hermitian'
135 : END IF
136 :
137 : ! A k-dependent nonunitary AO change requires both overlap and connection derivatives.
138 30 : db = 0.0_dp
139 30 : db(1, 1) = 0.13_dp
140 30 : db(2, 1) = CMPLX(0.17_dp, 0.23_dp, dp)
141 30 : db(1, 2) = CMPLX(-0.09_dp, 0.11_dp, dp)
142 210 : b = s + k*db
143 870 : h2 = MATMUL(CONJG(TRANSPOSE(b)), MATMUL(h, b))
144 450 : s2 = MATMUL(CONJG(TRANSPOSE(b)), b)
145 30 : dh2 = 0.0_dp
146 30 : ds2 = 0.0_dp
147 30 : d2 = 0.0_dp
148 : dh2(:, :, 1) = MATMUL(CONJG(TRANSPOSE(b)), MATMUL(dh(:, :, 1), b)) + &
149 2910 : MATMUL(CONJG(TRANSPOSE(db)), MATMUL(h, b)) + MATMUL(CONJG(TRANSPOSE(b)), MATMUL(h, db))
150 1230 : ds2(:, :, 1) = MATMUL(CONJG(TRANSPOSE(db)), b) + MATMUL(CONJG(TRANSPOSE(b)), db)
151 810 : d2(:, :, 1) = -iu*MATMUL(CONJG(TRANSPOSE(db)), b) + 7*s2
152 30 : CALL bloch_current(h2, s2, dh2, ds2, d2, e2, j2, status)
153 90 : IF (status /= 0 .OR. MAXVAL(ABS(e - e2)) > 1.e-12_dp) ERROR STOP 'Bloch basis covariance failed'
154 90 : DO i = 1, 2
155 90 : IF (ABS(j(i, i, 1) - j2(i, i, 1)) > 1.e-11_dp) ERROR STOP 'Band velocity basis dependence'
156 : END DO
157 660 : IF (MAXVAL(ABS(ABS(j) - ABS(j2))) > 1.e-11_dp) ERROR STOP 'Interband velocity gauge dependence'
158 : ! An energy-zero shift must include the k derivative of the metric.
159 30 : shift = 4.2_dp
160 840 : CALL bloch_current(h2 + shift*s2, s2, dh2 + shift*ds2, ds2, d2, e2, j2, status)
161 90 : IF (status /= 0 .OR. MAXVAL(ABS(e + shift - e2)) > 1.e-11_dp) ERROR STOP 'Energy shift failed'
162 660 : IF (MAXVAL(ABS(ABS(j) - ABS(j2))) > 1.e-10_dp) ERROR STOP 'Velocity energy-zero dependence'
163 210 : d2(:, :, 1) = 0.0_dp
164 30 : CALL bloch_current(h2, s2, dh2, ds2, d2, e2, j2, status)
165 32 : IF (status /= -3) ERROR STOP 'Inconsistent connection accepted'
166 : END DO
167 :
168 2 : dh = 0.0_dp
169 2 : ds = 0.0_dp
170 2 : d = 0.0_dp
171 2 : d(1, 2, 1) = iu
172 2 : d(2, 1, 1) = -iu
173 2 : CALL bloch_current(h, s, dh, ds, d, e, j, status)
174 2 : IF (status /= 0) ERROR STOP 'Finite Bloch limit failed'
175 2 : CALL projected_current(h, s, d, e2, j2, status)
176 44 : IF (status /= 0 .OR. MAXVAL(ABS(j - j2)) > 1.e-12_dp) ERROR STOP 'Finite commutator limit failed'
177 :
178 : ! Every matrix element within an exactly degenerate block has the same Fermi weight.
179 2 : h = 0.0_dp
180 2 : dh = 0.0_dp
181 2 : d = 0.0_dp
182 14 : dh(:, :, 1) = s
183 2 : dh(1, 2, 1) = iu
184 2 : dh(2, 1, 1) = -iu
185 2 : CALL bloch_current(h, s, dh, ds, d, e, j, status)
186 2 : IF (status /= 0) ERROR STOP 'Degenerate Bloch block failed'
187 14 : strength = SUM(ABS(j(:, :, 1))**2)*(-fermi_slope(0.0_dp, 0.0_dp, 0.0_dp, 0.1_dp))
188 2 : IF (ABS(strength - 10.0_dp) > 1.e-12_dp) ERROR STOP 'Intraband Drude weight failed'
189 2 : b = RESHAPE([iu, (1.0_dp, 0.0_dp), (1.0_dp, 0.0_dp), iu], [2, 2])/SQRT(2.0_dp)
190 58 : expected = MATMUL(CONJG(TRANSPOSE(b)), MATMUL(dh(:, :, 1), b))
191 14 : dh(:, :, 1) = expected
192 2 : CALL bloch_current(h, s, dh, ds, d, e, j, status)
193 14 : strength2 = SUM(ABS(j(:, :, 1))**2)*(-fermi_slope(0.0_dp, 0.0_dp, 0.0_dp, 0.1_dp))
194 2 : IF (status /= 0 .OR. ABS(strength - strength2) > 1.e-12_dp) ERROR STOP 'Degenerate response gauge failed'
195 2 : IF (ABS(fermi_slope(1.e-9_dp, -1.e-9_dp, 0.0_dp, 0.1_dp) + 2.5_dp) > 1.e-12_dp) THEN
196 0 : ERROR STOP 'Fermi near-degeneracy limit failed'
197 : END IF
198 2 : IF (ABS(fermi_slope(0.1_dp, -0.1_dp, 0.0_dp, 0.1_dp) + TANH(0.5_dp)/0.2_dp) > 1.e-12_dp) THEN
199 0 : ERROR STOP 'Fermi divided difference failed'
200 : END IF
201 2 : PRINT *, 'Bloch metric/connection, band velocity, gauge, degeneracy and Drude tests passed.'
202 2 : END SUBROUTINE check_bloch_current
203 :
204 : ! **************************************************************************************************
205 : !> \brief Validate Hall signs against Liouvillian and Berry-curvature references.
206 : ! **************************************************************************************************
207 2 : SUBROUTINE check_hall_response()
208 2 : USE ieee_arithmetic, ONLY: ieee_quiet_nan, &
209 : ieee_value
210 : USE kinds, ONLY: dp
211 : USE kubo_bloch, ONLY: bloch_current, &
212 : fermi_slope
213 : USE kubo_hall, ONLY: hall_response
214 :
215 : IMPLICIT NONE
216 : COMPLEX(dp), PARAMETER :: iu = (0.0_dp, 1.0_dp)
217 : REAL(dp), PARAMETER :: pi = 3.1415926535897932384626433832795_dp
218 : REAL(dp) :: energies(4), hall(3, 3), check(3, 3), reference(3, 3), rotation(3, 3)
219 : REAL(dp) :: analytic, error, gauge_error, mass, kx, ky, radius, omega, integrated, curvature, gamma
220 : REAL(dp) :: e2(2), masses(3), expected(3), unit(3)
221 : COMPLEX(dp) :: current(4, 4, 3), changed(4, 4, 3), u(4, 4)
222 : COMPLEX(dp) :: h2(2, 2), s2(2, 2), dh2(2, 2, 3), ds2(2, 2, 3), d2(2, 2, 3), j2(2, 2, 3)
223 : INTEGER :: a, b, i, j, imass, ix, iy, status
224 : INTEGER, PARAMETER :: mesh = 48
225 :
226 : ! A two-level response fixes the sign, occupation factor and relaxation denominator.
227 2 : j2 = 0.0_dp
228 2 : j2(1, 2, 1) = -iu
229 2 : j2(2, 1, 1) = -iu
230 2 : j2(1, 2, 2) = -1.0_dp
231 2 : j2(2, 1, 2) = 1.0_dp
232 2 : CALL hall_response([-0.7_dp, 0.7_dp], j2, 0.0_dp, 0.2_dp, 0.3_dp, hall, status)
233 2 : analytic = 2.0_dp*TANH(0.7_dp/0.4_dp)/(1.4_dp**2 + 0.3_dp**2)
234 2 : IF (status /= 0 .OR. ABS(hall(1, 2) - analytic) > 1.e-13_dp) ERROR STOP 'Two-level Hall sign or scale'
235 26 : IF (MAXVAL(ABS(hall + TRANSPOSE(hall))) /= 0.0_dp) ERROR STOP 'Hall not antisymmetric'
236 2 : CALL hall_response([0.0_dp, 0.0_dp], j2, 0.0_dp, 0.2_dp, 0.3_dp, hall, status)
237 26 : IF (status /= 0 .OR. MAXVAL(ABS(hall)) /= 0.0_dp) ERROR STOP 'Degenerate Hall block'
238 14 : DO i = 2, 14, 2
239 14 : radius = 10.0_dp**(-i)
240 42 : CALL hall_response([-radius, radius], j2, 0.0_dp, 0.2_dp, 0.3_dp, hall, status)
241 14 : analytic = 2.0_dp*TANH(radius/0.4_dp)/(4*radius**2 + 0.3_dp**2)
242 16 : IF (status /= 0 .OR. ABS(hall(1, 2)/analytic - 1.0_dp) > 1.e-12_dp) THEN
243 0 : ERROR STOP 'Near-degenerate Hall limit'
244 : END IF
245 : END DO
246 :
247 2 : energies = [-0.9_dp, -0.9_dp, 0.8_dp, 0.8_dp]
248 8 : DO a = 1, 3
249 30 : DO j = 1, 4
250 126 : DO i = 1, 4
251 120 : current(i, j, a) = CMPLX(SIN(REAL(i + 3*j + a, dp)), COS(REAL(2*i - j + a, dp)), dp)
252 : END DO
253 : END DO
254 248 : current(:, :, a) = current(:, :, a) - CONJG(TRANSPOSE(current(:, :, a)))
255 : END DO
256 2 : CALL hall_response(energies, current, 0.1_dp, 0.15_dp, 0.07_dp, hall, status)
257 26 : IF (status /= 0 .OR. MAXVAL(ABS(hall)) < 0.1_dp) ERROR STOP 'Nontrivial four-band response'
258 2 : CALL liouville_reference(energies, current, 0.1_dp, 0.15_dp, 0.07_dp, reference)
259 26 : error = MAXVAL(ABS(hall - 0.5_dp*(reference - TRANSPOSE(reference))))
260 2 : IF (error > 2.e-12_dp) ERROR STOP 'Full Liouvillian reference mismatch'
261 2 : PRINT '(A,ES14.6)', 'Hall Liouvillian residual: ', error
262 8 : DO a = 1, 3
263 6 : unit = 0.0_dp
264 6 : unit(a) = 1.0_dp
265 98 : IF (DOT_PRODUCT(unit, MATMUL(reference, unit)) < -1.e-12_dp) ERROR STOP 'Dissipation sign'
266 : END DO
267 :
268 2 : u = 0.0_dp
269 14 : u(1:2, 1:2) = RESHAPE([iu, (1.0_dp, 0.0_dp), (1.0_dp, 0.0_dp), iu], [2, 2])/SQRT(2.0_dp)
270 : u(3:4, 3:4) = RESHAPE([(0.6_dp, 0.0_dp), (0.0_dp, 0.8_dp), (0.0_dp, 0.8_dp), &
271 14 : (0.6_dp, 0.0_dp)], [2, 2])
272 8 : DO a = 1, 3
273 1136 : changed(:, :, a) = MATMUL(CONJG(TRANSPOSE(u)), MATMUL(current(:, :, a), u))
274 : END DO
275 2 : CALL hall_response(energies, changed, 0.1_dp, 0.15_dp, 0.07_dp, check, status)
276 26 : IF (status /= 0 .OR. MAXVAL(ABS(hall - check)) > 2.e-12_dp) ERROR STOP 'Degenerate gauge covariance'
277 10 : CALL hall_response(energies + 2.3_dp, current, 2.4_dp, 0.15_dp, 0.07_dp, check, status)
278 26 : IF (status /= 0 .OR. MAXVAL(ABS(hall - check)) > 2.e-12_dp) ERROR STOP 'Energy zero covariance'
279 128 : changed = current([3, 1, 4, 2], [3, 1, 4, 2], :)
280 10 : CALL hall_response(energies([3, 1, 4, 2]), changed, 0.1_dp, 0.15_dp, 0.07_dp, check, status)
281 26 : IF (status /= 0 .OR. MAXVAL(ABS(hall - check)) > 2.e-12_dp) ERROR STOP 'Band permutation covariance'
282 128 : CALL hall_response(energies, CONJG(current), 0.1_dp, 0.15_dp, 0.07_dp, check, status)
283 26 : IF (status /= 0 .OR. MAXVAL(ABS(hall + check)) > 2.e-12_dp) ERROR STOP 'Time-reversed Hall sign'
284 : rotation = RESHAPE([0.6_dp, 0.8_dp, 0.0_dp, -0.8_dp, 0.6_dp, 0.0_dp, &
285 2 : 0.0_dp, 0.0_dp, -1.0_dp], [3, 3])
286 2 : changed = 0.0_dp
287 8 : DO a = 1, 3
288 26 : DO b = 1, 3
289 384 : changed(:, :, a) = changed(:, :, a) + rotation(a, b)*current(:, :, b)
290 : END DO
291 : END DO
292 2 : CALL hall_response(energies, changed, 0.1_dp, 0.15_dp, 0.07_dp, check, status)
293 158 : reference = MATMUL(rotation, MATMUL(hall, TRANSPOSE(rotation)))
294 26 : IF (status /= 0 .OR. MAXVAL(ABS(reference - check)) > 2.e-12_dp) ERROR STOP 'Cartesian tensor covariance'
295 :
296 2 : CALL hall_response(energies, current, 0.0_dp, 0.0_dp, 0.1_dp, check, status)
297 2 : IF (status == 0) ERROR STOP 'Zero temperature accepted'
298 2 : CALL hall_response(energies, current, 0.0_dp, 0.1_dp, 0.0_dp, check, status)
299 2 : IF (status == 0) ERROR STOP 'Zero relaxation accepted'
300 2 : CALL hall_response(energies, current, ieee_value(0.0_dp, ieee_quiet_nan), 0.1_dp, 0.1_dp, check, status)
301 2 : IF (status == 0) ERROR STOP 'NaN chemical potential accepted'
302 2 : changed = current
303 2 : changed(1, 2, 1) = changed(1, 2, 1) + 0.1_dp
304 2 : CALL hall_response(energies, changed, 0.0_dp, 0.1_dp, 0.1_dp, check, status)
305 2 : IF (status == 0) ERROR STOP 'Non-anti-Hermitian current accepted'
306 :
307 : ! Two-band lattice model: compare every point against analytic Berry curvature,
308 : ! then integrate its two opposite Chern phases and a trivial phase.
309 2 : masses = [-1.0_dp, 1.0_dp, 3.0_dp]
310 2 : expected = [-1.0_dp, 1.0_dp, 0.0_dp]
311 2 : s2 = 0.0_dp
312 2 : s2(1, 1) = 1.0_dp
313 2 : s2(2, 2) = 1.0_dp
314 2 : ds2 = 0.0_dp
315 2 : d2 = 0.0_dp
316 2 : gamma = 1.e-3_dp
317 8 : DO imass = 1, 3
318 6 : mass = masses(imass)
319 6 : integrated = 0.0_dp
320 6 : curvature = 0.0_dp
321 6 : error = 0.0_dp
322 6 : gauge_error = 0.0_dp
323 294 : DO iy = 0, mesh - 1
324 288 : ky = 2*pi*REAL(iy, dp)/mesh
325 14118 : DO ix = 0, mesh - 1
326 13824 : kx = 2*pi*REAL(ix, dp)/mesh
327 13824 : h2(1, 1) = mass + COS(kx) + COS(ky)
328 13824 : h2(2, 2) = -h2(1, 1)
329 13824 : h2(1, 2) = SIN(kx) - iu*SIN(ky)
330 13824 : h2(2, 1) = CONJG(h2(1, 2))
331 13824 : dh2 = 0.0_dp
332 13824 : dh2(1, 1, 1) = -SIN(kx)
333 13824 : dh2(2, 2, 1) = SIN(kx)
334 13824 : dh2(1, 2, 1) = COS(kx)
335 13824 : dh2(2, 1, 1) = COS(kx)
336 13824 : dh2(1, 1, 2) = -SIN(ky)
337 13824 : dh2(2, 2, 2) = SIN(ky)
338 13824 : dh2(1, 2, 2) = -iu*COS(ky)
339 13824 : dh2(2, 1, 2) = iu*COS(ky)
340 13824 : CALL bloch_current(h2, s2, dh2, ds2, d2, e2, j2, status)
341 13824 : IF (status /= 0) ERROR STOP 'Chern-model Bloch current'
342 13824 : CALL hall_response(e2, j2, 0.0_dp, 0.01_dp, gamma, check, status)
343 13824 : IF (status /= 0) ERROR STOP 'Chern-model Hall kernel'
344 13824 : radius = SQRT(SIN(kx)**2 + SIN(ky)**2 + (mass + COS(kx) + COS(ky))**2)
345 13824 : omega = (COS(kx) + COS(ky) + mass*COS(kx)*COS(ky))/(2*radius**3)
346 13824 : analytic = -omega*4*radius**2/(4*radius**2 + gamma**2)*TANH(radius/0.02_dp)
347 13824 : error = MAX(error, ABS(check(1, 2) - analytic))
348 13824 : CALL ao_gauge_reference(kx, ky, h2, dh2, gamma, reference)
349 179712 : gauge_error = MAX(gauge_error, MAXVAL(ABS(check - reference)))
350 13824 : integrated = integrated + check(1, 2)*2*pi/mesh**2
351 14112 : curvature = curvature + omega*2*pi/mesh**2
352 : END DO
353 : END DO
354 6 : IF (error > 1.e-12_dp) ERROR STOP 'Analytic Berry-curvature response mismatch'
355 6 : IF (ABS(integrated - expected(imass)) > 1.e-6_dp) ERROR STOP 'Chern Hall limit'
356 6 : IF (ABS(curvature + expected(imass)) > 1.e-10_dp) ERROR STOP 'Independent Chern integral'
357 6 : IF (gauge_error > 1.e-12_dp) ERROR STOP 'Nonorthogonal AO Hall covariance'
358 6 : PRINT '(A,F5.1,3ES18.10)', 'Hall Chern model mass, sigma, Chern, residual: ', mass, integrated, curvature, error
359 8 : PRINT '(A,ES14.6)', 'Hall k-dependent nonorthogonal AO residual: ', gauge_error
360 : END DO
361 2 : PRINT *, 'Hall sign, Liouvillian, gauge, time reversal, degeneracy and Chern tests passed.'
362 : CONTAINS
363 : ! **************************************************************************************************
364 : !> \brief Transform the lattice model to a k-dependent nonorthogonal AO frame.
365 : !> \param kx First wave-vector component
366 : !> \param ky Second wave-vector component
367 : !> \param h Orthogonal Hamiltonian
368 : !> \param dh Its Cartesian derivatives
369 : !> \param gamma Relaxation energy
370 : !> \param hall Hall response in the transformed AO frame
371 : ! **************************************************************************************************
372 13824 : SUBROUTINE ao_gauge_reference(kx, ky, h, dh, gamma, hall)
373 : REAL(dp), INTENT(IN) :: kx, ky
374 : COMPLEX(dp), INTENT(IN) :: h(2, 2), dh(2, 2, 3)
375 : REAL(dp), INTENT(IN) :: gamma
376 : REAL(dp), INTENT(OUT) :: hall(3, 3)
377 :
378 : COMPLEX(dp) :: basis(2, 2), db(2, 2, 3), dhx(2, 2, 3), &
379 : dsx(2, 2, 3), dx(2, 2, 3), hx(2, 2), &
380 : jx(2, 2, 3), sx(2, 2)
381 : INTEGER :: a, status
382 : REAL(dp) :: e(2)
383 :
384 13824 : basis = 0.0_dp
385 13824 : basis(1, 1) = 1.1_dp + 0.1_dp*COS(kx)
386 13824 : basis(2, 2) = 0.9_dp + 0.1_dp*SIN(ky)
387 13824 : basis(1, 2) = 0.13_dp*EXP(iu*(kx + ky))
388 13824 : basis(2, 1) = 0.07_dp*EXP(iu*(2*kx - ky))
389 13824 : db = 0.0_dp
390 13824 : db(1, 1, 1) = -0.1_dp*SIN(kx)
391 13824 : db(2, 2, 2) = 0.1_dp*COS(ky)
392 41472 : db(1, 2, 1:2) = iu*basis(1, 2)
393 13824 : db(2, 1, 1) = 2*iu*basis(2, 1)
394 13824 : db(2, 1, 2) = -iu*basis(2, 1)
395 400896 : hx = MATMUL(CONJG(TRANSPOSE(basis)), MATMUL(h, basis))
396 207360 : sx = MATMUL(CONJG(TRANSPOSE(basis)), basis)
397 55296 : DO a = 1, 3
398 : dhx(:, :, a) = MATMUL(CONJG(TRANSPOSE(basis)), MATMUL(dh(:, :, a), basis)) + &
399 : MATMUL(CONJG(TRANSPOSE(db(:, :, a))), MATMUL(h, basis)) + &
400 4022784 : MATMUL(CONJG(TRANSPOSE(basis)), MATMUL(h, db(:, :, a)))
401 : dsx(:, :, a) = MATMUL(CONJG(TRANSPOSE(db(:, :, a))), basis) + &
402 1700352 : MATMUL(CONJG(TRANSPOSE(basis)), db(:, :, a))
403 1133568 : dx(:, :, a) = -iu*MATMUL(CONJG(TRANSPOSE(db(:, :, a))), basis) + REAL(a, dp)*sx
404 : END DO
405 13824 : CALL bloch_current(hx, sx, dhx, dsx, dx, e, jx, status)
406 13824 : IF (status /= 0) ERROR STOP 'Nonorthogonal AO Chern model'
407 13824 : CALL hall_response(e, jx, 0.0_dp, 0.01_dp, gamma, hall, status)
408 13824 : IF (status /= 0) ERROR STOP 'Nonorthogonal AO Hall kernel'
409 2 : END SUBROUTINE ao_gauge_reference
410 :
411 : ! **************************************************************************************************
412 : !> \brief Invert the full complex Liouvillian in a nondiagonal four-level frame.
413 : !> \param e Eigenvalues
414 : !> \param currents Anti-Hermitian currents in the diagonal frame
415 : !> \param mu Chemical potential
416 : !> \param kt Temperature
417 : !> \param gamma Relaxation rate
418 : !> \param sigma Complete real DC tensor without physical prefactors
419 : ! **************************************************************************************************
420 2 : SUBROUTINE liouville_reference(e, currents, mu, kt, gamma, sigma)
421 : REAL(dp), INTENT(IN) :: e(4)
422 : COMPLEX(dp), INTENT(IN) :: currents(4, 4, 3)
423 : REAL(dp), INTENT(IN) :: mu, kt, gamma
424 : REAL(dp), INTENT(OUT) :: sigma(3, 3)
425 :
426 : COMPLEX(dp) :: basis(4, 4), drive(4, 4), h(4, 4), &
427 : liouvillian(16, 16), response(4, 4), &
428 : rhs(16, 3), trace, v(4, 4, 3)
429 : INTEGER :: a, b, column, i, info, j, k, l, &
430 : pivot(16), row
431 :
432 10 : DO j = 1, 4
433 42 : DO i = 1, 4
434 40 : basis(i, j) = EXP(iu*2*pi*REAL((i - 1)*(j - 1), dp)/4)/2
435 : END DO
436 : END DO
437 2 : h = 0.0_dp
438 10 : DO i = 1, 4
439 10 : h(i, i) = e(i)
440 : END DO
441 338 : h = MATMUL(basis, MATMUL(h, CONJG(TRANSPOSE(basis))))
442 8 : DO a = 1, 3
443 1254 : v(:, :, a) = MATMUL(basis, MATMUL(iu*currents(:, :, a), CONJG(TRANSPOSE(basis))))
444 30 : DO j = 1, 4
445 126 : DO i = 1, 4
446 120 : drive(i, j) = fermi_slope(e(i), e(j), mu, kt)*iu*currents(i, j, a)
447 : END DO
448 : END DO
449 1014 : drive = MATMUL(basis, MATMUL(drive, CONJG(TRANSPOSE(basis))))
450 8 : rhs(:, a) = RESHAPE(drive, [16])
451 : END DO
452 2 : liouvillian = 0.0_dp
453 10 : DO j = 1, 4
454 42 : DO i = 1, 4
455 32 : row = i + 4*(j - 1)
456 32 : liouvillian(row, row) = gamma
457 168 : DO l = 1, 4
458 672 : DO k = 1, 4
459 512 : column = k + 4*(l - 1)
460 512 : IF (j == l) liouvillian(row, column) = liouvillian(row, column) + iu*h(i, k)
461 640 : IF (i == k) liouvillian(row, column) = liouvillian(row, column) - iu*h(l, j)
462 : END DO
463 : END DO
464 : END DO
465 : END DO
466 2 : CALL zgesv(16, 3, liouvillian, 16, pivot, rhs, 16, info)
467 2 : IF (info /= 0) ERROR STOP 'Liouvillian solve failed'
468 8 : DO b = 1, 3
469 6 : response = RESHAPE(rhs(:, b), [4, 4])
470 26 : DO a = 1, 3
471 378 : trace = SUM(v(:, :, a)*TRANSPOSE(response))
472 24 : sigma(a, b) = -REAL(trace, dp)
473 : END DO
474 : END DO
475 2 : END SUBROUTINE liouville_reference
476 : END SUBROUTINE check_hall_response
477 :
478 : ! **************************************************************************************************
479 : !> \brief Periodic QWZ phases, physical windows, dual gauges, nonorthogonal bases and spinor TR.
480 :
481 : END MODULE kubo_test_cases
482 :
483 2 : PROGRAM kubo_unittest
484 2 : USE kubo_test_cases, ONLY: check_projected_current, check_bloch_current, check_hall_response
485 : IMPLICIT NONE
486 2 : CALL check_projected_current()
487 2 : CALL check_bloch_current()
488 2 : CALL check_hall_response()
489 2 : END PROGRAM kubo_unittest
|
