Skip to content

Fix time-reversal symmetry in plane-wave basis Hartree matrix elements v_h and dump_hartree - #46

Open
jaemolihm wants to merge 1 commit into
AbInitioQHub:developfrom
jaemolihm:fix/vh-hartree-time-reversal
Open

Fix time-reversal symmetry in plane-wave basis Hartree matrix elements v_h and dump_hartree#46
jaemolihm wants to merge 1 commit into
AbInitioQHub:developfrom
jaemolihm:fix/vh-hartree-time-reversal

Conversation

@jaemolihm

@jaemolihm jaemolihm commented Jul 18, 2026

Copy link
Copy Markdown

Problem

coqui.post_proc.dump_hartree (plane-wave Hartree matrix elements VH_skij) was wrong for systems with time-reversal symmetry but no symmorphic inversion (inversion about the grid origin), such as silicon with inversion center at (1/8, 1/8, 1/8).

Fix

Updated src/hamiltonian/v_h.hpp

  • before: conj(FFT[c @ -S·g])(r) = u₋ₛₖ(-r) → density n_{Sk}(-r) (inverted about the origin)
  • after:  conj(FFT[c @ +S·g])(r) = u₋ₛₖ(r)  → density n_{Sk}(r) (correct)

Same as

if (kp_trev(ik))
S_tsk_full(it, is, ik, nda::ellipsis{}) = nda::conj(S_tsk(it, is, ik_ibz, nda::ellipsis{}));

Test

New fixture qe_si333_sym (diamond Si, 3×3×3, force_symmorphic = .true.) + two check_Hartree sections.

  • J: from the THC ERIs (hf.evaluate)
  • J2: from the plane-wave path (hamilt::Vhartreev_h, the code under test)
section ‖J − J2‖ before ‖J − J2‖ after
si333_sym (full ρ) 0.098 (FAIL) 2e-7
si333_sym_diag (diagonal ρ) 0.098 (FAIL) 2e-7

Existing LiH/GaAs fixtures do not cover this because LiH inversion center at origin; GaAs does not use symmetry.

@jaemolihm jaemolihm changed the title fix(hartree): correct time-reversal density reconstruction in v_h Fix time-reversal symmetry in plane-wave basis Hartree matrix elements v_h and dump_hartree Jul 18, 2026
hamilt::v_h passed kp_trev(k) to transform_k2g, which applies the sign
flip sg=-1 and places the plane-wave coefficients at -S*g. Combined with
the separate real-space conjugation of the orbital, this double-applies
the sign and produces the per-k density inverted about the FFT grid
origin, n_{Sk}(-r) instead of n_{Sk}(r).

For crystals whose inversion center is not at the grid origin (e.g.
diamond Si, inversion center at the bond midpoint), the time-reversed
half of the BZ is then misregistered relative to the non-TR half. The
total charge (G=0) is preserved, but the G!=0 structure partially
cancels, halving the Hartree matrix elements (e.g. Si Gamma band 3:
2.65 eV vs QE 5.30 eV). Only the plane-wave v_h path is affected; the
THC J[Dm] path used for HF/GW energies is correct.

Fix: pass false to transform_k2g in both v_h overloads. Time reversal is
already applied via the real-space conjugation of the orbital (full-
density path) and is unnecessary for |psi|^2 (diagonal-density path).
This matches the idiom used elsewhere (unfold_bz.cpp): spatial rotation
via transform_k2g(false, ...) and time reversal via a separate conj.

Test: new fixture qe_si333_sym (diamond Si, 3x3x3, force_symmorphic ->
27 kpts reduce to 4 IBZ with 13 time-reversal pairs and no inversion op)
with two check_Hartree sections (full and diagonal density). Before fix
||J - J2|| = 0.098 (FAIL); after fix O(1e-7). Existing LiH/GaAs fixtures
do not catch this: LiH's inversion center is at the origin (masks the
misregistration) and the GaAs fixture has no TR pairs.
@jaemolihm
jaemolihm force-pushed the fix/vh-hartree-time-reversal branch from 47983c2 to 397ec30 Compare July 18, 2026 20:48
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant