From 0464cc52efea82ded965877be75a3e11b1cdd4e7 Mon Sep 17 00:00:00 2001 From: Surefire618 Date: Wed, 12 Aug 2026 14:33:21 +0200 Subject: [PATCH] Keep the full complex velocity in kappa_offdiag get_kappa_offdiag discarded the imaginary part of the off-diagonal group velocity before squaring it: buf_vel(:, i, j) = real(cv0, r8) justified by "I can take the real part since at the end we sum over both modes and the imaginary components disappear". The premise holds -- sum_ij Im v_ij is zero because v_ij = conj(v_ji) -- but the quantity accumulated afterwards is quadratic: buf_velsq(:, :, i, j) += lo_outerproduct(v1, v1) and sum_ij Im v^a Im v^b does not vanish. Dropping it leaves Re v^a Re v^b alone, which is not invariant under the arbitrary per-mode phases zheev returns, so kappa_offdiag depends on the eigenvector basis. Measured on three materials, a random phase gauge moved it by 27-59 % while leaving kappa_sma bit-identical. Keeping cv0 complex and taking the real part of the outer product gives Re v^a Re v^b + Im v^a Im v^b, the gauge invariant bilinear (lo_complex_outerproduct conjugates its first argument). kappa_offdiag is then bit-identical across gauges; kappa_sma, frequencies, linewidths, lifetimes, heat capacities and group velocities are unchanged. kappa_offdiag itself moves: 0.00481941 -> 0.00830067 (KI, Pm-3m), 0.05503408 -> 0.06626228 (Rb2O, Fm-3m), 0.07273678/0.14298927 -> 0.16748196/0.37433862 (KPTe2, R-3m) W/mK at 300 K on a 6^3 mesh. --- src/thermal_conductivity/kappa.f90 | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/thermal_conductivity/kappa.f90 b/src/thermal_conductivity/kappa.f90 index 71747e38..e262b76a 100644 --- a/src/thermal_conductivity/kappa.f90 +++ b/src/thermal_conductivity/kappa.f90 @@ -95,7 +95,7 @@ subroutine get_kappa_offdiag(dr, qp, uc, fc, temperature, classical, mem, mw, ka real(r8), dimension(3, 3), intent(out) :: kappa_offdiag !> The off diagonal group velocity - real(r8), dimension(:, :, :), allocatable :: buf_vel + complex(r8), dimension(:, :, :), allocatable :: buf_vel !> The off diagonal group velocity, squared real(r8), dimension(:, :, :, :), allocatable :: buf_velsq !> The qpoint @@ -116,7 +116,7 @@ subroutine get_kappa_offdiag(dr, qp, uc, fc, temperature, classical, mem, mw, ka complex(r8), dimension(:, :, :), allocatable :: buf_grad_dynmat complex(r8), dimension(:, :), allocatable :: kronegv, buf_egv, buf_egw, buf_cm0, buf_cm1, buf_cm2 complex(r8), dimension(3) :: cv0 - real(r8), dimension(3) :: v0, v1 + complex(r8), dimension(3) :: v0, v1 integer :: a1, a2, ia, ib, ic, ix, iy, iz, k, iop, i, ii, j, jj ! Some buffers @@ -207,9 +207,12 @@ subroutine get_kappa_offdiag(dr, qp, uc, fc, temperature, classical, mem, mw, ka cv0 = buf_cm2(ii, :) ! remove tiny numbers. cv0 = lo_chop(cv0, 1E-10/(lo_groupvel_Hartreebohr_to_ms/1000)) - ! I can take the real part since at the end we sum over - ! both modes and the imaginary components disappear. - buf_vel(:, i, j) = real(cv0, r8) + ! Keep the full complex velocity. The imaginary parts do cancel + ! in sum_ij v_ij, but the quantity accumulated below is quadratic + ! and sum_ij Im v^a Im v^b does not vanish. Discarding it makes + ! kappa_offdiag depend on the arbitrary phases of the eigenvectors + ! returned by zheev. + buf_vel(:, i, j) = cv0 end do end do @@ -226,7 +229,10 @@ subroutine get_kappa_offdiag(dr, qp, uc, fc, temperature, classical, mem, mw, ka do k = 1, qp%ip(iq)%n_full_point iop = qp%ip(iq)%operation_full_point(k) v1 = matmul(uc%sym%op(abs(iop))%m, v0) - buf_velsq(:, :, i, j) = buf_velsq(:, :, i, j) + lo_outerproduct(v1, v1) + ! lo_outerproduct conjugates its first argument, so the real part + ! here is Re v^a Re v^b + Im v^a Im v^b, the gauge invariant form. + buf_velsq(:, :, i, j) = buf_velsq(:, :, i, j) + & + real(lo_outerproduct(v1, v1), r8) end do end do end do