kubo_unittest.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
MODULE kubo_test_cases
   IMPLICIT NONE
   PRIVATE
   PUBLIC :: check_projected_current, check_bloch_current, check_hall_response
CONTAINS
 
! **************************************************************************************************
!> \brief Validate projected currents and their AO-metric covariance.
! **************************************************************************************************
   SUBROUTINE check_projected_current()
      USE kinds, ONLY: dp
      USE kubo_projected, ONLY: projected_current
      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), &
         s(4, 4), s2(4, 4), x(4, 4, 3), x2(4, 4, 3)
      INTEGER                                            :: a, c, i, status
      REAL(dp)                                           :: e(4), e2(4)
 
      h = 0.0_dp
      s = 0.0_dp
      x = 0.0_dp
      DO i = 1, 4
         s(i, i) = 1.0_dp
      END DO
      ! Kramers-degenerate spinor model with imaginary SOC matrix elements.
      h(2, 2) = 0.8_dp
      h(4, 4) = 0.8_dp
      h(1, 2) = cmplx(0.0_dp, 0.3_dp, dp)
      h(2, 1) = conjg(h(1, 2))
      h(3, 4) = conjg(h(1, 2))
      h(4, 3) = h(1, 2)
      x(1, 2, 1) = 1.0_dp
      x(2, 1, 1) = 1.0_dp
      x(3, 4, 1) = 1.0_dp
      x(4, 3, 1) = 1.0_dp
      x(1, 1, 2) = -0.5_dp
      x(2, 2, 2) = 0.5_dp
      x(3, 3, 2) = -0.5_dp
      x(4, 4, 2) = 0.5_dp
      x(:, :, 3) = 2*x(:, :, 1) - x(:, :, 2)
      CALL projected_current(h, s, x, e, j, status)
      IF (status /= 0) error stop 'Spinor current failed'
      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'
      DO a = 1, 3
         IF (maxval(abs(j(:, :, a) + conjg(transpose(j(:, :, a))))) > 1.e-12_dp) error stop 'Current is not anti-Hermitian'
         DO c = 1, 3
            norm(a, c) = sum(conjg(j(:, :, a))*j(:, :, c))
         END DO
      END DO
      IF (abs(norm(1, 1) - 4.0_dp) > 1.e-12_dp) error stop 'Analytic current strength failed'
      ! Complex nonorthogonal basis mixing spin and orbitals, including a degenerate subspace.
      b = s
      b(1, 1) = 1.7_dp
      b(2, 2) = 0.8_dp
      b(3, 1) = cmplx(0.2_dp, 0.4_dp, dp)
      b(1, 4) = cmplx(-0.3_dp, 0.1_dp, dp)
      b(2, 3) = 0.4_dp
      h2 = matmul(conjg(transpose(b)), matmul(h, b))
      s2 = matmul(conjg(transpose(b)), b)
      DO a = 1, 3
         x2(:, :, a) = matmul(conjg(transpose(b)), matmul(x(:, :, a), b)) + real(a, dp)*s2
      END DO
      CALL projected_current(h2, s2, x2, e2, j2, status)
      IF (status /= 0 .OR. maxval(abs(e - e2)) > 1.e-12_dp) error stop 'Nonorthogonal current spectrum failed'
      DO a = 1, 3
         DO c = 1, 3
            norm2(a, c) = sum(conjg(j2(:, :, a))*j2(:, :, c))
         END DO
      END DO
      IF (maxval(abs(norm - norm2)) > 1.e-11_dp) error stop 'Basis/origin/gauge invariance failed'
      s(1, 1) = -1.0_dp
      CALL projected_current(h, s, x, e, j, status)
      IF (status /= -2) error stop 'Indefinite overlap accepted'
      s(1, 1) = 1.0_dp
      x(1, 2, 1) = 2.0_dp
      CALL projected_current(h, s, x, e, j, status)
      IF (status /= -1) error stop 'Non-Hermitian position accepted'
      print *, 'Projected spinor current, metric and gauge tests passed.'
   END SUBROUTINE check_projected_current
 
! **************************************************************************************************
!> \brief Validate Bloch currents, degenerate blocks and k-dependent AO gauges.
! **************************************************************************************************
   SUBROUTINE check_bloch_current()
      USE kinds, ONLY: dp
      USE kubo_bloch, ONLY: bloch_current, &
                            bloch_current_from_frame, &
                            fermi_slope
      USE kubo_projected, ONLY: projected_current
      COMPLEX(dp), PARAMETER                             :: iu = (0.0_dp, 1.0_dp)
 
      COMPLEX(dp) :: b(2, 2), d(2, 2, 3), d2(2, 2, 3), db(2, 2), dh(2, 2, 3), dh2(2, 2, 3), &
         ds(2, 2, 3), ds2(2, 2, 3), expected(2, 2), frame(2, 2), h(2, 2), h2(2, 2), j(2, 2, 3), &
         j2(2, 2, 3), s(2, 2), s2(2, 2)
      INTEGER                                            :: i, ik, status
      REAL(dp)                                           :: e(2), e2(2), k, residuals(2), root, &
                                                            shift, speed, strength, strength2
 
      s = 0.0_dp
      s(1, 1) = 1.0_dp
      s(2, 2) = 1.0_dp
      DO ik = -7, 7
         k = real(ik, dp)/5
         h = k*k*s
         h(1, 1) = h(1, 1) + 0.3_dp
         h(2, 2) = h(2, 2) - 0.3_dp
         h(1, 2) = 0.4_dp*sin(k) - iu*0.2_dp*cos(k)
         h(2, 1) = conjg(h(1, 2))
         dh = 0.0_dp
         dh(:, :, 1) = 2*k*s
         dh(1, 2, 1) = 0.4_dp*cos(k) + iu*0.2_dp*sin(k)
         dh(2, 1, 1) = conjg(dh(1, 2, 1))
         ds = 0.0_dp
         d = 0.0_dp
         CALL bloch_current(h, s, dh, ds, d, e, j, status, frame)
         IF (status /= 0) error stop 'Bloch current failed'
         CALL bloch_current_from_frame(h, s, dh, ds, d, e, frame, j2, 1.e-9_dp, residuals, status)
         IF (status /= 0 .OR. maxval(abs(j2 - j)) > 1.e-12_dp) error stop 'Reused eigenframe failed'
         CALL bloch_current_from_frame(h, s, dh, ds, d, e + 0.1_dp, frame, j2, 1.e-9_dp, residuals, status)
         IF (status /= -4) error stop 'Wrong frame eigenvalues accepted'
         CALL bloch_current_from_frame(h, s, dh, ds, d, e, 1.1_dp*frame, j2, 1.e-9_dp, residuals, status)
         IF (status /= -4) error stop 'Nonorthonormal frame accepted'
         root = sqrt(0.09_dp + 0.16_dp*sin(k)**2 + 0.04_dp*cos(k)**2)
         IF (maxval(abs(e - [k*k - root, k*k + root])) > 1.e-12_dp) error stop 'Analytic bands failed'
         speed = 0.12_dp*sin(k)*cos(k)/root
         IF (abs(iu*j(1, 1, 1) - (2*k - speed)) > 1.e-12_dp) error stop 'Lower band velocity failed'
         IF (abs(iu*j(2, 2, 1) - (2*k + speed)) > 1.e-12_dp) error stop 'Upper band velocity failed'
         IF (maxval(abs(j(:, :, 1) + conjg(transpose(j(:, :, 1))))) > 1.e-12_dp) THEN
            error stop 'Bloch current not anti-Hermitian'
         END IF
 
         ! A k-dependent nonunitary AO change requires both overlap and connection derivatives.
         db = 0.0_dp
         db(1, 1) = 0.13_dp
         db(2, 1) = cmplx(0.17_dp, 0.23_dp, dp)
         db(1, 2) = cmplx(-0.09_dp, 0.11_dp, dp)
         b = s + k*db
         h2 = matmul(conjg(transpose(b)), matmul(h, b))
         s2 = matmul(conjg(transpose(b)), b)
         dh2 = 0.0_dp
         ds2 = 0.0_dp
         d2 = 0.0_dp
         dh2(:, :, 1) = matmul(conjg(transpose(b)), matmul(dh(:, :, 1), b)) + &
                        matmul(conjg(transpose(db)), matmul(h, b)) + matmul(conjg(transpose(b)), matmul(h, db))
         ds2(:, :, 1) = matmul(conjg(transpose(db)), b) + matmul(conjg(transpose(b)), db)
         d2(:, :, 1) = -iu*matmul(conjg(transpose(db)), b) + 7*s2
         CALL bloch_current(h2, s2, dh2, ds2, d2, e2, j2, status)
         IF (status /= 0 .OR. maxval(abs(e - e2)) > 1.e-12_dp) error stop 'Bloch basis covariance failed'
         DO i = 1, 2
            IF (abs(j(i, i, 1) - j2(i, i, 1)) > 1.e-11_dp) error stop 'Band velocity basis dependence'
         END DO
         IF (maxval(abs(abs(j) - abs(j2))) > 1.e-11_dp) error stop 'Interband velocity gauge dependence'
         ! An energy-zero shift must include the k derivative of the metric.
         shift = 4.2_dp
         CALL bloch_current(h2 + shift*s2, s2, dh2 + shift*ds2, ds2, d2, e2, j2, status)
         IF (status /= 0 .OR. maxval(abs(e + shift - e2)) > 1.e-11_dp) error stop 'Energy shift failed'
         IF (maxval(abs(abs(j) - abs(j2))) > 1.e-10_dp) error stop 'Velocity energy-zero dependence'
         d2(:, :, 1) = 0.0_dp
         CALL bloch_current(h2, s2, dh2, ds2, d2, e2, j2, status)
         IF (status /= -3) error stop 'Inconsistent connection accepted'
      END DO
 
      dh = 0.0_dp
      ds = 0.0_dp
      d = 0.0_dp
      d(1, 2, 1) = iu
      d(2, 1, 1) = -iu
      CALL bloch_current(h, s, dh, ds, d, e, j, status)
      IF (status /= 0) error stop 'Finite Bloch limit failed'
      CALL projected_current(h, s, d, e2, j2, status)
      IF (status /= 0 .OR. maxval(abs(j - j2)) > 1.e-12_dp) error stop 'Finite commutator limit failed'
 
      ! Every matrix element within an exactly degenerate block has the same Fermi weight.
      h = 0.0_dp
      dh = 0.0_dp
      d = 0.0_dp
      dh(:, :, 1) = s
      dh(1, 2, 1) = iu
      dh(2, 1, 1) = -iu
      CALL bloch_current(h, s, dh, ds, d, e, j, status)
      IF (status /= 0) error stop 'Degenerate Bloch block failed'
      strength = sum(abs(j(:, :, 1))**2)*(-fermi_slope(0.0_dp, 0.0_dp, 0.0_dp, 0.1_dp))
      IF (abs(strength - 10.0_dp) > 1.e-12_dp) error stop 'Intraband Drude weight failed'
      b = reshape([iu, (1.0_dp, 0.0_dp), (1.0_dp, 0.0_dp), iu], [2, 2])/sqrt(2.0_dp)
      expected = matmul(conjg(transpose(b)), matmul(dh(:, :, 1), b))
      dh(:, :, 1) = expected
      CALL bloch_current(h, s, dh, ds, d, e, j, status)
      strength2 = sum(abs(j(:, :, 1))**2)*(-fermi_slope(0.0_dp, 0.0_dp, 0.0_dp, 0.1_dp))
      IF (status /= 0 .OR. abs(strength - strength2) > 1.e-12_dp) error stop 'Degenerate response gauge failed'
      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
         error stop 'Fermi near-degeneracy limit failed'
      END IF
      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
         error stop 'Fermi divided difference failed'
      END IF
      print *, 'Bloch metric/connection, band velocity, gauge, degeneracy and Drude tests passed.'
   END SUBROUTINE check_bloch_current
 
! **************************************************************************************************
!> \brief Validate Hall signs against Liouvillian and Berry-curvature references.
! **************************************************************************************************
   SUBROUTINE check_hall_response()
      USE ieee_arithmetic, ONLY: ieee_quiet_nan, &
                                 ieee_value
      USE kinds, ONLY: dp
      USE kubo_bloch, ONLY: bloch_current, &
                            fermi_slope
      USE kubo_hall, ONLY: hall_response
 
      IMPLICIT NONE
      COMPLEX(dp), PARAMETER :: iu = (0.0_dp, 1.0_dp)
      REAL(dp), PARAMETER :: pi = 3.1415926535897932384626433832795_dp
      REAL(dp) :: energies(4), hall(3, 3), check(3, 3), reference(3, 3), rotation(3, 3)
      REAL(dp) :: analytic, error, gauge_error, mass, kx, ky, radius, omega, integrated, curvature, gamma
      REAL(dp) :: e2(2), masses(3), expected(3), unit(3)
      COMPLEX(dp) :: current(4, 4, 3), changed(4, 4, 3), u(4, 4)
      COMPLEX(dp) :: h2(2, 2), s2(2, 2), dh2(2, 2, 3), ds2(2, 2, 3), d2(2, 2, 3), j2(2, 2, 3)
      INTEGER :: a, b, i, j, imass, ix, iy, status
      INTEGER, PARAMETER :: mesh = 48
 
      ! A two-level response fixes the sign, occupation factor and relaxation denominator.
      j2 = 0.0_dp
      j2(1, 2, 1) = -iu
      j2(2, 1, 1) = -iu
      j2(1, 2, 2) = -1.0_dp
      j2(2, 1, 2) = 1.0_dp
      CALL hall_response([-0.7_dp, 0.7_dp], j2, 0.0_dp, 0.2_dp, 0.3_dp, hall, status)
      analytic = 2.0_dp*tanh(0.7_dp/0.4_dp)/(1.4_dp**2 + 0.3_dp**2)
      IF (status /= 0 .OR. abs(hall(1, 2) - analytic) > 1.e-13_dp) error stop 'Two-level Hall sign or scale'
      IF (maxval(abs(hall + transpose(hall))) /= 0.0_dp) error stop 'Hall not antisymmetric'
      CALL hall_response([0.0_dp, 0.0_dp], j2, 0.0_dp, 0.2_dp, 0.3_dp, hall, status)
      IF (status /= 0 .OR. maxval(abs(hall)) /= 0.0_dp) error stop 'Degenerate Hall block'
      DO i = 2, 14, 2
         radius = 10.0_dp**(-i)
         CALL hall_response([-radius, radius], j2, 0.0_dp, 0.2_dp, 0.3_dp, hall, status)
         analytic = 2.0_dp*tanh(radius/0.4_dp)/(4*radius**2 + 0.3_dp**2)
         IF (status /= 0 .OR. abs(hall(1, 2)/analytic - 1.0_dp) > 1.e-12_dp) THEN
            error stop 'Near-degenerate Hall limit'
         END IF
      END DO
 
      energies = [-0.9_dp, -0.9_dp, 0.8_dp, 0.8_dp]
      DO a = 1, 3
         DO j = 1, 4
            DO i = 1, 4
               current(i, j, a) = cmplx(sin(real(i + 3*j + a, dp)), cos(real(2*i - j + a, dp)), dp)
            END DO
         END DO
         current(:, :, a) = current(:, :, a) - conjg(transpose(current(:, :, a)))
      END DO
      CALL hall_response(energies, current, 0.1_dp, 0.15_dp, 0.07_dp, hall, status)
      IF (status /= 0 .OR. maxval(abs(hall)) < 0.1_dp) error stop 'Nontrivial four-band response'
      CALL liouville_reference(energies, current, 0.1_dp, 0.15_dp, 0.07_dp, reference)
      error = maxval(abs(hall - 0.5_dp*(reference - transpose(reference))))
      IF (error > 2.e-12_dp) error stop 'Full Liouvillian reference mismatch'
      print '(A,ES14.6)', 'Hall Liouvillian residual: ', error
      DO a = 1, 3
         unit = 0.0_dp
         unit(a) = 1.0_dp
         IF (dot_product(unit, matmul(reference, unit)) < -1.e-12_dp) error stop 'Dissipation sign'
      END DO
 
      u = 0.0_dp
      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)
      u(3:4, 3:4) = reshape([(0.6_dp, 0.0_dp), (0.0_dp, 0.8_dp), (0.0_dp, 0.8_dp), &
                             (0.6_dp, 0.0_dp)], [2, 2])
      DO a = 1, 3
         changed(:, :, a) = matmul(conjg(transpose(u)), matmul(current(:, :, a), u))
      END DO
      CALL hall_response(energies, changed, 0.1_dp, 0.15_dp, 0.07_dp, check, status)
      IF (status /= 0 .OR. maxval(abs(hall - check)) > 2.e-12_dp) error stop 'Degenerate gauge covariance'
      CALL hall_response(energies + 2.3_dp, current, 2.4_dp, 0.15_dp, 0.07_dp, check, status)
      IF (status /= 0 .OR. maxval(abs(hall - check)) > 2.e-12_dp) error stop 'Energy zero covariance'
      changed = current([3, 1, 4, 2], [3, 1, 4, 2], :)
      CALL hall_response(energies([3, 1, 4, 2]), changed, 0.1_dp, 0.15_dp, 0.07_dp, check, status)
      IF (status /= 0 .OR. maxval(abs(hall - check)) > 2.e-12_dp) error stop 'Band permutation covariance'
      CALL hall_response(energies, conjg(current), 0.1_dp, 0.15_dp, 0.07_dp, check, status)
      IF (status /= 0 .OR. maxval(abs(hall + check)) > 2.e-12_dp) error stop 'Time-reversed Hall sign'
      rotation = reshape([0.6_dp, 0.8_dp, 0.0_dp, -0.8_dp, 0.6_dp, 0.0_dp, &
                          0.0_dp, 0.0_dp, -1.0_dp], [3, 3])
      changed = 0.0_dp
      DO a = 1, 3
         DO b = 1, 3
            changed(:, :, a) = changed(:, :, a) + rotation(a, b)*current(:, :, b)
         END DO
      END DO
      CALL hall_response(energies, changed, 0.1_dp, 0.15_dp, 0.07_dp, check, status)
      reference = matmul(rotation, matmul(hall, transpose(rotation)))
      IF (status /= 0 .OR. maxval(abs(reference - check)) > 2.e-12_dp) error stop 'Cartesian tensor covariance'
 
      CALL hall_response(energies, current, 0.0_dp, 0.0_dp, 0.1_dp, check, status)
      IF (status == 0) error stop 'Zero temperature accepted'
      CALL hall_response(energies, current, 0.0_dp, 0.1_dp, 0.0_dp, check, status)
      IF (status == 0) error stop 'Zero relaxation accepted'
      CALL hall_response(energies, current, ieee_value(0.0_dp, ieee_quiet_nan), 0.1_dp, 0.1_dp, check, status)
      IF (status == 0) error stop 'NaN chemical potential accepted'
      changed = current
      changed(1, 2, 1) = changed(1, 2, 1) + 0.1_dp
      CALL hall_response(energies, changed, 0.0_dp, 0.1_dp, 0.1_dp, check, status)
      IF (status == 0) error stop 'Non-anti-Hermitian current accepted'
 
      ! Two-band lattice model: compare every point against analytic Berry curvature,
      ! then integrate its two opposite Chern phases and a trivial phase.
      masses = [-1.0_dp, 1.0_dp, 3.0_dp]
      expected = [-1.0_dp, 1.0_dp, 0.0_dp]
      s2 = 0.0_dp
      s2(1, 1) = 1.0_dp
      s2(2, 2) = 1.0_dp
      ds2 = 0.0_dp
      d2 = 0.0_dp
      gamma = 1.e-3_dp
      DO imass = 1, 3
         mass = masses(imass)
         integrated = 0.0_dp
         curvature = 0.0_dp
         error = 0.0_dp
         gauge_error = 0.0_dp
         DO iy = 0, mesh - 1
            ky = 2*pi*real(iy, dp)/mesh
            DO ix = 0, mesh - 1
               kx = 2*pi*real(ix, dp)/mesh
               h2(1, 1) = mass + cos(kx) + cos(ky)
               h2(2, 2) = -h2(1, 1)
               h2(1, 2) = sin(kx) - iu*sin(ky)
               h2(2, 1) = conjg(h2(1, 2))
               dh2 = 0.0_dp
               dh2(1, 1, 1) = -sin(kx)
               dh2(2, 2, 1) = sin(kx)
               dh2(1, 2, 1) = cos(kx)
               dh2(2, 1, 1) = cos(kx)
               dh2(1, 1, 2) = -sin(ky)
               dh2(2, 2, 2) = sin(ky)
               dh2(1, 2, 2) = -iu*cos(ky)
               dh2(2, 1, 2) = iu*cos(ky)
               CALL bloch_current(h2, s2, dh2, ds2, d2, e2, j2, status)
               IF (status /= 0) error stop 'Chern-model Bloch current'
               CALL hall_response(e2, j2, 0.0_dp, 0.01_dp, gamma, check, status)
               IF (status /= 0) error stop 'Chern-model Hall kernel'
               radius = sqrt(sin(kx)**2 + sin(ky)**2 + (mass + cos(kx) + cos(ky))**2)
               omega = (cos(kx) + cos(ky) + mass*cos(kx)*cos(ky))/(2*radius**3)
               analytic = -omega*4*radius**2/(4*radius**2 + gamma**2)*tanh(radius/0.02_dp)
               error = max(error, abs(check(1, 2) - analytic))
               CALL ao_gauge_reference(kx, ky, h2, dh2, gamma, reference)
               gauge_error = max(gauge_error, maxval(abs(check - reference)))
               integrated = integrated + check(1, 2)*2*pi/mesh**2
               curvature = curvature + omega*2*pi/mesh**2
            END DO
         END DO
         IF (error > 1.e-12_dp) error stop 'Analytic Berry-curvature response mismatch'
         IF (abs(integrated - expected(imass)) > 1.e-6_dp) error stop 'Chern Hall limit'
         IF (abs(curvature + expected(imass)) > 1.e-10_dp) error stop 'Independent Chern integral'
         IF (gauge_error > 1.e-12_dp) error stop 'Nonorthogonal AO Hall covariance'
         print '(A,F5.1,3ES18.10)', 'Hall Chern model mass, sigma, Chern, residual: ', mass, integrated, curvature, error
         print '(A,ES14.6)', 'Hall k-dependent nonorthogonal AO residual: ', gauge_error
      END DO
      print *, 'Hall sign, Liouvillian, gauge, time reversal, degeneracy and Chern tests passed.'
   CONTAINS
! **************************************************************************************************
!> \brief Transform the lattice model to a k-dependent nonorthogonal AO frame.
!> \param kx First wave-vector component
!> \param ky Second wave-vector component
!> \param h Orthogonal Hamiltonian
!> \param dh Its Cartesian derivatives
!> \param gamma Relaxation energy
!> \param hall Hall response in the transformed AO frame
! **************************************************************************************************
      SUBROUTINE ao_gauge_reference(kx, ky, h, dh, gamma, hall)
      REAL(dp), INTENT(IN)                               :: kx, ky
      COMPLEX(dp), INTENT(IN)                            :: h(2, 2), dh(2, 2, 3)
      REAL(dp), INTENT(IN)                               :: gamma
      REAL(dp), INTENT(OUT)                              :: hall(3, 3)
 
      COMPLEX(dp)                                        :: basis(2, 2), db(2, 2, 3), dhx(2, 2, 3), &
                                                            dsx(2, 2, 3), dx(2, 2, 3), hx(2, 2), &
                                                            jx(2, 2, 3), sx(2, 2)
      INTEGER                                            :: a, status
      REAL(dp)                                           :: e(2)
 
         basis = 0.0_dp
         basis(1, 1) = 1.1_dp + 0.1_dp*cos(kx)
         basis(2, 2) = 0.9_dp + 0.1_dp*sin(ky)
         basis(1, 2) = 0.13_dp*exp(iu*(kx + ky))
         basis(2, 1) = 0.07_dp*exp(iu*(2*kx - ky))
         db = 0.0_dp
         db(1, 1, 1) = -0.1_dp*sin(kx)
         db(2, 2, 2) = 0.1_dp*cos(ky)
         db(1, 2, 1:2) = iu*basis(1, 2)
         db(2, 1, 1) = 2*iu*basis(2, 1)
         db(2, 1, 2) = -iu*basis(2, 1)
         hx = matmul(conjg(transpose(basis)), matmul(h, basis))
         sx = matmul(conjg(transpose(basis)), basis)
         DO a = 1, 3
            dhx(:, :, a) = matmul(conjg(transpose(basis)), matmul(dh(:, :, a), basis)) + &
                           matmul(conjg(transpose(db(:, :, a))), matmul(h, basis)) + &
                           matmul(conjg(transpose(basis)), matmul(h, db(:, :, a)))
            dsx(:, :, a) = matmul(conjg(transpose(db(:, :, a))), basis) + &
                           matmul(conjg(transpose(basis)), db(:, :, a))
            dx(:, :, a) = -iu*matmul(conjg(transpose(db(:, :, a))), basis) + real(a, dp)*sx
         END DO
         CALL bloch_current(hx, sx, dhx, dsx, dx, e, jx, status)
         IF (status /= 0) error stop 'Nonorthogonal AO Chern model'
         CALL hall_response(e, jx, 0.0_dp, 0.01_dp, gamma, hall, status)
         IF (status /= 0) error stop 'Nonorthogonal AO Hall kernel'
      END SUBROUTINE ao_gauge_reference
 
! **************************************************************************************************
!> \brief Invert the full complex Liouvillian in a nondiagonal four-level frame.
!> \param e Eigenvalues
!> \param currents Anti-Hermitian currents in the diagonal frame
!> \param mu Chemical potential
!> \param kt Temperature
!> \param gamma Relaxation rate
!> \param sigma Complete real DC tensor without physical prefactors
! **************************************************************************************************
      SUBROUTINE liouville_reference(e, currents, mu, kt, gamma, sigma)
      REAL(dp), INTENT(IN)                               :: e(4)
      COMPLEX(dp), INTENT(IN)                            :: currents(4, 4, 3)
      REAL(dp), INTENT(IN)                               :: mu, kt, gamma
      REAL(dp), INTENT(OUT)                              :: sigma(3, 3)
 
      COMPLEX(dp)                                        :: basis(4, 4), drive(4, 4), h(4, 4), &
                                                            liouvillian(16, 16), response(4, 4), &
                                                            rhs(16, 3), trace, v(4, 4, 3)
      INTEGER                                            :: a, b, column, i, info, j, k, l, &
                                                            pivot(16), row
 
         DO j = 1, 4
            DO i = 1, 4
               basis(i, j) = exp(iu*2*pi*real((i - 1)*(j - 1), dp)/4)/2
            END DO
         END DO
         h = 0.0_dp
         DO i = 1, 4
            h(i, i) = e(i)
         END DO
         h = matmul(basis, matmul(h, conjg(transpose(basis))))
         DO a = 1, 3
            v(:, :, a) = matmul(basis, matmul(iu*currents(:, :, a), conjg(transpose(basis))))
            DO j = 1, 4
               DO i = 1, 4
                  drive(i, j) = fermi_slope(e(i), e(j), mu, kt)*iu*currents(i, j, a)
               END DO
            END DO
            drive = matmul(basis, matmul(drive, conjg(transpose(basis))))
            rhs(:, a) = reshape(drive, [16])
         END DO
         liouvillian = 0.0_dp
         DO j = 1, 4
            DO i = 1, 4
               row = i + 4*(j - 1)
               liouvillian(row, row) = gamma
               DO l = 1, 4
                  DO k = 1, 4
                     column = k + 4*(l - 1)
                     IF (j == l) liouvillian(row, column) = liouvillian(row, column) + iu*h(i, k)
                     IF (i == k) liouvillian(row, column) = liouvillian(row, column) - iu*h(l, j)
                  END DO
               END DO
            END DO
         END DO
         CALL zgesv(16, 3, liouvillian, 16, pivot, rhs, 16, info)
         IF (info /= 0) error stop 'Liouvillian solve failed'
         DO b = 1, 3
            response = reshape(rhs(:, b), [4, 4])
            DO a = 1, 3
               trace = sum(v(:, :, a)*transpose(response))
               sigma(a, b) = -real(trace, dp)
            END DO
         END DO
      END SUBROUTINE liouville_reference
   END SUBROUTINE check_hall_response
 
! **************************************************************************************************
!> \brief Periodic QWZ phases, physical windows, dual gauges, nonorthogonal bases and spinor TR.
 
END MODULE kubo_test_cases
 
PROGRAM kubo_unittest
   USE kubo_test_cases, ONLY: check_projected_current, check_bloch_current, check_hall_response
   IMPLICIT NONE
   CALL check_projected_current()
   CALL check_bloch_current()
   CALL check_hall_response()
END PROGRAM kubo_unittest