Repository navigation
Port dmftproj from Fortran to Python (tracking) #7
Description
Activity
Strategy and stacked PRs
Port dmftproj to Python one
case.*output at a time, each PR proven against the Fortran dmftproj output matrix by matrix. Tail-backward, so the Fortran shrinks from the output end and every step leaves a working converter. Generators were written self-contained first (so each could be validated in isolation), then the shared machinery was extracted once. The final PR swaps in a pure-Python driver and removesfortran/dmftproj.Review and merge in order:
- [test] Add spin-orbit + spin-polarized Wien2k converter test (CaOs2) #4: SOC + spin-polarized golden test (CaOs2). The characterization harness.
- [wien2k] Pure-Python case.symqmc generator (first dmftproj port step) #5:
case.symqmc. - [wien2k] pure-Python case.oubwin band-window generator #8:
case.oubwin(and the almblm reader). - [wien2k] pure-Python case.ctqmcout correlated-shell projectors #9:
case.ctqmcout, the correlated-shell projector core. - [wien2k] pure-Python case.sympar partial-shell symmetry matrices #10:
case.sympar. - [wien2k] pure-Python case.parproj partial projectors #11:
case.parproj. - [wien2k] complete the port to machine precision: shared module, fromfile, proj_mode, outband #12: extract the shared
_dmftprojmodule (deduplication, 2260 to 1699 lines; adds unit tests of the isolated primitives). - [wien2k] pure-Python dmftproj driver; remove the Fortran executable #13: pure-Python
run_dmftprojdriver, removes the Fortran executable.
Every PR ships a test against the committed dmftproj reference; #12 adds isolated unit tests; #13 runs driver to converter to
h5diffend to end. The tree is test-covered at each step, including after the refactor.On numerical accuracy
Agreement is 1e-6 to 1e-7, not bit-exact. dmftproj builds the cubic/fromfile basis transform through a single-precision
CMPLXcast (set_ang_trans.f), sosymqmc/ctqmcout/sympar/parprojcarry ~1e-7 noise; the Python is double precision and reproduces the values to that floor.oubwin(integer band windows) is byte-exact. Going tighter would mean mimicking the Fortran single-precision truncation at each site bit for bit; if that is wanted, we would still need to redo parts of the basis-transform path to match the cast exactly (or fix the precision in Fortran first, cf. the fromfile-precision fix discussed under dft_tools #148).- added 5 commits that reference this issue
on Jun 16, 2026 Port complete: machine precision on every output
dmftproj is fully ported to Python and reproduces the Fortran to machine
precision. The Fortran precision fix is kept separate and stacked below.Stacked PRs, review/merge order:
- [dmftproj] eliminate single-precision loss in the projector construction #14: Fortran precision fix (separate, reproducible). Every single-precision
CMPLXcast in dmftproj built aCOMPLEX(KIND=8)fromREAL(KIND=8)args
withoutKIND=8(the spinrot phase in outputqmc.f, the Loewdin eigenvalue,
the k-weights), a ~1e-7 loss; plus the shipped cubic templates stored at full
double precision. Supersedes the fromfile-only [dmftproj] build fromfile basis representation in double precision #6. - [test] Add spin-orbit + spin-polarized Wien2k converter test (CaOs2) #4: SOC + spin-polarized golden test (CaOs2 harness).
- [wien2k] Pure-Python case.symqmc generator (first dmftproj port step) #5 symqmc, [wien2k] pure-Python case.oubwin band-window generator #8 oubwin (+ almblm reader), [wien2k] pure-Python case.ctqmcout correlated-shell projectors #9 ctqmcout (projector core), [wien2k] pure-Python case.sympar partial-shell symmetry matrices #10
sympar, [wien2k] pure-Python case.parproj partial projectors #11 parproj. Each reproduces the released dmftproj. - [wien2k] complete the port to machine precision: shared module, fromfile, proj_mode, outband #12: complete the port to machine precision: extract the shared
_dmftproj
module, switch to exact harmonics and the same LAPACK/BLAS routines dmftproj
uses, and port the remaining functionality the converter exposes downstream:
the fromfile spin-mixing basis, the proj_mode 1/2 band-index window, and
case.outband(band mode, feedingconvert_bands_input). - [wien2k] pure-Python dmftproj driver; remove the Fortran executable #13: pure-Python
run_dmftprojdriver, removes the Fortran executable.
Accuracy against the precision-fixed dmftproj (#14), all integer fields identical:
output accuracy output accuracy oubwin byte-identical sympar 1.9e-14 symqmc 1.5e-14 parproj 1.5e-14 ctqmcout 1.5e-14 fromfile 1.9e-14 outband 2.8e-16 proj_mode 1.5e-14 Two findings worth recording
The Loewdin overlap must be full rank. ctqmcout/parproj/outband
orthonormalize the correlated projectors with O^(-1/2), O = D D^H. For a narrow
energy window with more correlated spin-orbitals than bands (CaOs2: 20 vs 17) O
is rank-deficient; its near-null eigenvectors are non-unique and differ across
LAPACK builds, so O^(-1/2) amplifies the last-bit difference. dmftproj is
byte-identical to itself across BLAS, but no cross-language port matches it there
without changing the (preserved) algorithm. A physically sensible window (enough
bands) makes O non-singular and the result exact. The tests use full-rank
windows; the narrow-window result is a numerical property of that degenerate
case, not the port.Match the Fortran's exact operations. The projector D is bit-identical to the
Fortran; matching ctqmcout then required calling the identical LAPACK routine
(scipyZHEEV('V','U'), not numpy's divide-and-conquerZHEEVD) and the same
ZGEMMtrans flags rather than numpy conj-transpose copies.The outband fixture was generated by a Wien2k band run (
lapw1/lapwso/lapw2 -band) plusdmftproj -bandon a full-rank window. No downstream behavior
changes: the converter and its file interface are untouched.- [dmftproj] eliminate single-precision loss in the projector construction #14: Fortran precision fix (separate, reproducible). Every single-precision
Addendum: closing a verification gap. The generators were validated only on the spin-orbit case (CaOs2); the non-spin-orbit path (ifSO=0, the common DMFT case e.g. SrVO3) was emitting SO-sized 2(2l+1) matrices instead of (2l+1) and is now fixed and verified in #15 (ctqmcout 1.4e-15, symqmc 2.1e-14, sympar/parproj machine precision, against a precision-fixed dmftproj run without -so). New review order inserts #15 between #12 and the capstone #13; all PRs still target unstable. Full suite is 12/12 (numpy generators + TRIQS converter/driver tests).
Coverage pass before the capstone surfaced a real bug and closed the remaining gaps. All at machine precision against precision-fixed dmftproj.
Odd-l parity bug (PR #16). dmat took the rotation parity from det(krotm), but krotm in case.dmftsym is the proper part of the operation (det +1 even for improper rotations); the parity is iprop. The Wigner-D factor is (-1)^l, so for even l the sign is +1 either way and the d-only fixtures never exposed it. For odd l the symmetry matrices came out with the wrong overall sign. Fixed in symqmc, sympar and the mixing spinor representation. New p (l=1) and f (l=3) tests pin the sign; p-shell symqmc differed by 2.0 before, 1.2e-14 after.
f-shell rank-deficiency. Os has no f weight near E_F, so the f projector overlap is rank-deficient on the standard bands (~12 near-null eigenvalues, cond 8e9). Recomputing the almblm with a raised lapw1/lapwso band cutoff puts a high-energy window on bands with real l=3 character; the overlap is then full rank (cond 1.3) and ctqmcout matches at 2.6e-14. Same resolution as the narrow-window SO case: use a full-rank window.
End-to-end at machine precision (PR #13). Two new driver -> converter -> h5diff tests: non-SO partial (convert_dft_input + convert_parproj_input) and wide-window SO, both 1e-12. The narrow-window SO driver test stays at 1e-6, the documented rank-deficiency floor.
Full wien2k suite 15/15 in the TRIQS 3.3.1 container.
Correction to the previous note: the narrow-window SO driver test is removed, not kept. That window is rank-deficient (overlap condition number infinite: more correlated spin-orbitals than bands), so a Python-vs-Fortran h5diff on it cannot exceed ~1e-6 and is not a meaningful machine-precision target. The driver is instead verified end to end on full-rank windows only (wien2k_socfull_convert SO and wien2k_noso_convert non-SO, both 1e-12). The narrow physical case stays locked exactly by wien2k_soc_convert, which runs the converter on the Fortran output (no Python regeneration, so no rank-deficiency). No 1e-6 tolerance remains anywhere.
Hi @krystophny, thank you for opening this issue and the related PRs so far. It's very well-documented and easy to follow. It will be nice to have a pure Python version of dmftproj!
In this port, can I also ask about how you are handling spin-orbit coupling? I learned last year that there is some bug in Wien2k 14 that dmftproj accommodated for, but since this bug has been fixed in Wien2k >14, there's now a bug in dmftproj. Can you comment on this issue?
The conventions used for spin-orbit coupling are important for when we do the symmetrization over the entire Brillouin zone to get the k-summation correct. A colleague of mine did some benchmarking with the Rutgers Wien2k+eDMFT implementation and found some discrepancies when considering the spin-orbit case.
Best,
HarryFYI also see TRIQS/dft_tools#148 (comment) on dft_tools. This should be mentioned here.
Reacted by Harry LaBollita@the-hampel thanks for taking a look! Unfortunately I have only access to Wien2k 14 so maybe someone with a newer version can check?
Follow-up on the spin-orbit / Wien2k-version question: we now have access to Wien2k 11.1 through 24.1 (14.2, 16.1, 19.1, 21.1, 24.1) and can test the port against any of them ourselves, so no external check is needed.
Multi-version results:
- The pure-Python port reproduces the Fortran dmftproj for every version tested; the gauge-invariant projector output (orbital occupations, local-Hamiltonian eigenvalues) is version-independent. The only version-dependent item is the 19.x almblm formatting change (dft_tools#127), which carries no data loss for cases without local orbitals in the affected l-channels.
- On the SOC symmetry question: the dmftproj symmetrization error (cf. dft_tools#148) does reproduce, but only for non-centrosymmetric spin-polarized + SO cells. Centrosymmetric cases (CaOs2, Sr2MgOsO6) are unaffected, which is why it stayed hidden in earlier checks. It is version-independent (identical drift across 16.1/19.1/21.1), i.e. a dmftproj convention issue, not a Wien2k regression. Localization so far: the intra-atomic symmetrization (including time reversal) is correct; the inter-atomic symmetry mapping is where it breaks. The originally proposed ephase sign flip does not fix it. Still under investigation.
Hi @krystophny, thanks for the update. It's nice to know that the "bug" is localized. Do you happen to have a set of notes with these SOC conventions that are being used in dmftproj. I'd imagine the Python implementation will likely become more readable.
Goal
Replace
fortran/dmftprojwith Python, one output at a time, leaving the converter and its file interface unchanged.dmftproj is the last Fortran component in the Wien2k path. It reads the Wien2k output (
case.almblm,case.struct,case.dmftsym,case.indmftpr, k-list) and writes the filestriqs_dftkit.wien2k.Converterconsumes. Each of those files is an independent target: generate it in Python, prove it reproduces the Fortran output, and the converter neither knows nor cares which side produced it. Nof90wrap, no C bindings, no change to how anything is wired.Method
Each step is one PR that adds a pure-Python generator for one
case.*file and is proven two ways:convert_dft_inputso a ported stage cannot move the HDF5.Order is tail-backward: the self-contained symmetry and bookkeeping outputs first, then the projector core, so the Fortran shrinks from the output end and every intermediate state is a working converter.
Targets
Files the converter reads from dmftproj (
.struct,.outputs,.pmatcome from Wien2k directly, so the port leaves them alone):.symqmcset_ang_trans,setsym,timeinv,outputqmc.sympar.oubwin*outbwin.ctqmcoutset_projections,orthogonal,rot_projectmat,outputqmc.parprojorthogonal_wannier,rot_projectmat.outbandoutbandStages
1. Symmetry and window outputs.
.symparreuses the.symqmcspinor machinery over the full shell set; small..oubwin*is index bookkeeping over the energy window and k-list, no projector math. Both are self-contained givencase.dmftsymand the window.2. The projectors (
.ctqmcout). Readcase.almblm(the Wien2kl, m, k, bandprojectors), restrict to the correlated energy window, orthonormalize within that window (Lowdin), and rotate to the chosen angular basis. This is the bulk of dmftproj and will land as several PRs: readalmblm; assemble the raw projectors; orthonormalize; basis-rotate and write..parprojfollows from the same machinery for the uncorrelated states.3. Charge feedback.
density,symmetrize_mat,rot_densbuild the local density matrix and the charge-density correction for the self-consistency loop (Driver.write_charge_correction). A one-shot conversion does not need them; charge-self-consistent DFT+DMFT does.4. Retire the Fortran. Port
.outbandif band plotting is wanted, move the driver loop (dmftproj.f) intodriver.py, and dropfortran/dmftprojfrom the build once every output is Python-native.Open
.symqmcvs partial.sympar) suggests.symqmcand.symparshare one symmetry module; [wien2k] Pure-Python case.symqmc generator (first dmftproj port step) #5 is structured for that..ctqmcoutorthonormalization fixes the numerical contract for everything downstream; it deserves its own design discussion before implementation.In flight
.symqmcgenerator (full replacement of the symqmc-writing).fromfiledouble-precision fix (found while porting.symqmc).