Repository navigation
Port dmftproj from Fortran to Python (tracking) #7
Description
Goal
Replace fortran/dmftproj with 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 files triqs_dftkit.wien2k.Converter consumes. 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. No f90wrap, 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:
- against the Fortran dmftproj output on the same input, matrix by matrix (the methodology in [wien2k] Pure-Python case.symqmc generator (first dmftproj port step) #5);
- end to end through the converter golden test [test] Add spin-orbit + spin-polarized Wien2k converter test (CaOs2) #4, which locks
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, .pmat come from Wien2k directly, so the port leaves them alone):
| file | content | Fortran routines | status |
|---|---|---|---|
.symqmc |
correlated-shell symmetry matrices | set_ang_trans, setsym, timeinv, outputqmc |
done, #5 |
.sympar |
partial-projector symmetry matrices | same, applied to all included shells | next |
.oubwin* |
correlated band-window indices per k | outbwin |
next |
.ctqmcout |
the projectors | set_projections, orthogonal, rot_projectmat, outputqmc |
core |
.parproj |
partial (uncorrelated) projectors | orthogonal_wannier, rot_projectmat |
core |
.outband |
projectors on a band path | outband |
optional |
Stages
1. Symmetry and window outputs. .sympar reuses the .symqmc spinor 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 given case.dmftsym and the window.
2. The projectors (.ctqmcout). Read case.almblm (the Wien2k l, m, k, band projectors), 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: read almblm; assemble the raw projectors; orthonormalize; basis-rotate and write. .parproj follows from the same machinery for the uncorrelated states.
3. Charge feedback. density, symmetrize_mat, rot_dens build 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 .outband if band plotting is wanted, move the driver loop (dmftproj.f) into driver.py, and drop fortran/dmftproj from the build once every output is Python-native.
Open
- The basis split (correlated
.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
- [test] Add spin-orbit + spin-polarized Wien2k converter test (CaOs2) #4: SOC + spin-polarized converter golden test (the harness).
- [wien2k] Pure-Python case.symqmc generator (first dmftproj port step) #5:
.symqmcgenerator (full replacement of the symqmc-writing). - [dmftproj] build fromfile basis representation in double precision #6: dmftproj
fromfiledouble-precision fix (found while porting.symqmc).
Activity
krystophny commentedon Jun 16, 2026
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:
case.symqmc.case.oubwin(and the almblm reader).case.ctqmcout, the correlated-shell projector core.case.sympar.case.parproj._dmftprojmodule (deduplication, 2260 to 1699 lines; adds unit tests of the isolated primitives).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).[wien2k] exact double-precision reproduction (with the dmftproj preci…
[wien2k] exact double-precision reproduction (with the dmftproj preci…
[wien2k] port the fromfile mixing basis and proj_mode 1/2 band-index …
[wien2k] port case.outband band-structure projectors (-band mode)
[wien2k] pure-Python dmftproj driver; remove the Fortran executable
krystophny commentedon 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:
CMPLXcast in dmftproj built aCOMPLEX(KIND=8)fromREAL(KIND=8)argswithout
KIND=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.
sympar, [wien2k] pure-Python case.parproj partial projectors #11 parproj. Each reproduces the released dmftproj.
_dmftprojmodule, 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).run_dmftprojdriver, removes the Fortran executable.Accuracy against the precision-fixed dmftproj (#14), all integer fields identical:
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
(scipy
ZHEEV('V','U'), not numpy's divide-and-conquerZHEEVD) and the sameZGEMMtrans 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 behaviorchanges: the converter and its file interface are untouched.
[wien2k] pure-Python dmftproj driver; remove the Fortran executable
[wien2k] pure-Python dmftproj driver; remove the Fortran executable
krystophny commentedon Jun 16, 2026
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).
krystophny commentedon Jun 16, 2026
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.
krystophny commentedon Jun 16, 2026
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.
harrisonlabollita commentedon Jun 23, 2026
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,
Harry
the-hampel commentedon Jun 23, 2026
FYI also see TRIQS/dft_tools#148 (comment) on dft_tools. This should be mentioned here.
krystophny commentedon Jun 23, 2026
@the-hampel thanks for taking a look! Unfortunately I have only access to Wien2k 14 so maybe someone with a newer version can check?
krystophny commentedon Jun 24, 2026
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:
harrisonlabollita commentedon Jun 24, 2026
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.
[wien2k] pure-Python case.oubwin band-window generator
[wien2k] pure-Python case.oubwin band-window generator
[wien2k] pure-Python case.ctqmcout correlated-shell projectors
[wien2k] Pure-Python case.oubwin band-window generator (#8)
[wien2k] pure-Python case.ctqmcout correlated-shell projectors