[wien2k] Charge self-consistent DFT+DMFT driver

[harrisonlabollita]

Summary

This PR provides a Wien2k Python driver that lets modest's DftDriver run charge self-consistent DFT+DMFT calculations. This allows the TRIQS user to take control of the DFT+DMFT similar to the QE and VASP workflows entirely from their pythons script. We include a worked SrVO3 example to parallel the VASP example.

Wien2k is a chain of Fortran programs launched by the x script rather than a
single executable, and its own DMFT support inverts the control flow: run_lapw calls out to a user python script from inside its cycle, whereas we want DftDriver to be the caller. Thus, we write the SCF loop in python in place of run_lapw.

run_update_stage writes the DMFT density matrices to case.qdmft, runs
lapw2 -qdmft, then lcore and mixer for the new density, lapw0 and lapw1 for
the potential and eigenvectors that go with it, and finally rebuilds the
projectors so the HDF5 archive is current when it returns.

Contents

python/triqs_dftkit/wien2k/driver.py	the driver: SCF loop, run_initial_stage, run_update_stage, run_dft_only, scf_converged
doc/examples/wien2k_csc_svo/	SrVO3 CSC example: struct, indmftpr, script, README

Notes for review

• Scope is serial, non-magnetic (SP=0, SO=0). lapw2 -qdmft has no parallel
  version. The impurity solver still uses all ranks.
• The SCF is restart-aware: a converged case.scf is reused and a partial one
  continued, since modest calls run_initial_stage unconditionally and keeps no
  DFT state across a resumed run. force_scf=True redoes it.

Testing

A SrVO3 csc run to convergence using the shipped example; A companion tutorial will be added on the ModEST side which will link the documentation between the two; run_dft_only() against a plain run_lapw SCF for comparison.

[the-hampel]

Review of the Wien2k CSC driver

Attribution: this review was produced by Claude (Opus 5) running in Claude Code, driven by @the-hampel. It is a static review plus local experiments — I have no Wien2k installation here, so anything that depends on Wien2k's own conventions is explicitly marked verify. Everything else was checked by running code.

Method. Replayed dmftproj.f's case.indmftpr reader against the shipped example; round-tripped _write_qdmft against the in-tree fixture test/python/wien2k/SrVO3.oubwin; and exercised the full SCF loop, run_update_stage and all three restart branches against stub x / dmftproj scripts (using the seam the class docstring advertises). Cross-checked against triqs_dftkit.vasp.driver, triqs_dftkit.qe.driver, triqs_modest.dft_driver.DftDriver, modest/c++/triqs_modest/chargedensity.hpp and dft_tools' SumkDFT.calc_density_correction(dm_type='wien2k').

Overall this is solid work and the best-documented driver in the tree — the "why" comments (_scf_tags record pairing, _save_old, the three qdmft conventions) are genuinely worth having. Two blockers, both concrete and both cheap to fix.

Verified correct — no need to re-check these

• case.qdmft format matches dft_tools' calc_density_correction(dm_type='wien2k') byte-for-byte: mu in Ry, beta in Ry⁻¹, per-k nn + nn rows of Re Im, the blank throwaway record, correner raw in eV. Round-tripped with a token reader; all values exact.
• The N_k convention is right. modest::charge_density_correction (chargedensity.hpp) returns the full band-basis density matrix, zero-padded to max_k N_nu(0,k), with index 0 = that k-point's own lowest in-window band. So writing the full matrix (not the deviation, unlike VASP/QE), unweighted, sliced :nn, is correct. The fixture has non-uniform windows so this is actually exercised.
• N_k is not mutated. VASP and QE both clobber the caller's array in place; nice that this one doesn't, and the comment saying why is worth keeping.
• _scf_tags record pairing — verified against all three :ENE variants (**INFO**** / *WARNING** / **********) and an :ENE-without-:DIS first cycle. Pairing by record rather than zipping two lists is the right call.
• A single window line in indmftpr is correct for this dmftproj (dmftproj.f, proj_mode = 0).
• The converter always overwrites dft_input, so _regenerate_projectors really does refresh the archive on every update.
• mpi.report(..., stderr=True) is supported by both TRIQS mpi backends.
• Program ordering in run_update_stage matches the docstring exactly (traced): lapw2 -qdmft, lcore, concat, mixer, concat scfm, lapw0, lapw1, lapw2 -almd, dmftproj, convert.

---

1. Blockers

1.1 SrVO3.indmftpr puts the correlated t2g shell on Sr, not V

Replaying dmftproj.f's reader against the shipped file:

sort 1 = Sr (ATOM 1, MULT=1):  lsort=[0,0,2,0]  lnreps=[0,0,2,0]  correps=01  SO=0   <-- CORRELATED
sort 2 = V  (ATOM 2, MULT=1):  lsort=[1,0,0,0]  lnreps=[0,0,0,0]                     <-- only s
sort 3 = O  (ATOM 3, MULT=3):  lsort=[0,1,0,0]  lnreps=[0,0,0,0]

The struct order is Sr, V, O(x3), and nmult = 1 1 3 confirms the mapping. The README says "we project onto the V-t2g shell" and the script builds a three-orbital Kanamori with U = 4.5, J = 0.65 for vanadium. Because Sr-d t2g is also three orbitals, the run completes and produces silently wrong physics — no error anywhere.

Fix: swap the sort-1 and sort-2 blocks (Sr becomes 0 0 0 0 / 0 0 0 0; V gets 0 0 2 0 / 0 0 2 0 / 01 / 0). fortran/dmftproj/SRC_templates/case.indmftpr in this repo already has the correct shape to copy from.

1.2 .restart is written and never removed

Traced through two restart cycles with the stub Wien2k:

restarted update stage   -> .restart present AFTER the mixer that consumed it: True
another update stage     -> .restart present: True

_mixing_flush_pending is one-shot, but the file is not. After a single flush, every subsequent mixer call in the run — and every future run in that directory — sees .restart. If mixer honours it, which the whole design assumes, the CSC loop can never build mixing state. _prepare_fresh_scf doesn't clear a stale one either, so a crash poisons the directory permanently.

Fix: os.remove('.restart') immediately after the mixer call in _run_mixer, and clear it in _prepare_fresh_scf. Better still, drop .restart in favour of removing case.broyd* — unambiguous, what run_lapw does, and already what _prepare_fresh_scf does.

Verify: who actually consumes .restart in SRC_mixer? If nothing does, _flush_mixing_history is a no-op that emits a reassuring warning while the inherited Broyden history is still used — worse than not having it, because the warning tells the user the problem was handled.

---

2. Should fix

2.1 Complex cases are silently unsupported

_cmplx gates in1c / in2c naming, _check_inputs and _set_in2_mode, but -c is never passed to x lapw1, x lapw2, lapw2 -almd or lapw2 -qdmft, and x does not auto-detect. Either thread the flag through _run_x or reject self._cmplx up front in _check_inputs — the current state implies support that isn't there.

2.2 Restart hole in run_initial_stage

On a resumed CSC run it reuses case.scf and goes straight to lapw2 -almd with no lapw0 / lapw1 (confirmed in the trace). If the interrupted run died between mixer and lapw1 completing — a real window, that's mixer + lapw0 + lapw1 of wall time — case.vector belongs to the pre-update potential while case.clmsum is post-update, so the projectors and H(k) are silently inconsistent with the density on disk.

The code already knows it is resuming (_has_dmft_cycles() is computed right there for the mixing flush); re-running lapw0 + lapw1 in that branch is cheap next to the solver.

2.3 _write_qdmft doesn't cross-check nn against the archive

The only guard is nn > Nk_avg.shape[1], i.e. greater than max_k n_orbitals. If case.oubwin and dft_input/n_orbitals ever disagree at some k with nn smaller than the archive value, the top-left nn x nn sub-block is written with no error and no warning. dft_tools asserts equality in the same place (assert nb == self.n_orbitals[ik, ntoi[bname]], sumk_dft.py); worth mirroring.

2.4 Eint_m_dc multiplicity is not handled

dft_tools' doc/guide/dftdmft_selfcons.rst is explicit: correnerg must be multiplied by the multiplicity of equivalent correlated atoms, and summed over inequivalent ones, before it goes into case.qdmft. _write_qdmft writes it raw. SrVO3 has a single V so the shipped example is fine, but any multi-atom cell gets a wrong total energy with nothing to indicate it. At minimum document this next to the three other convention notes in the _write_qdmft docstring; better, scale by MULT from corr_shells.

2.5 Test coverage

This is the weakest part of the PR: 831 lines of new driver logic, zero tests, CMakeLists.txt untouched.

• The precedent exists and is registered in CTest — Py_vasp_charge_update runs test/python/vasp/charge_update/test_gamma_writer.py, which tests the analogous VASP charge-update writer with no VASP anywhere.
• Roughly 80% of this driver is testable with no Wien2k at all. I covered it in ~150 lines of throwaway script while writing this review: _write_qdmft (format, :nn slicing on the non-uniform window, spin average, no-mutation, excluded-k rejection), _read_oubwin, _scf_tags / _converged / _scf_history_on_disk on a synthetic case.scf, _set_in2_mode, _cmplx, the _check_inputs messages — and the entire SCF loop, run_update_stage ordering and all three restart branches via stub x / dmftproj scripts.
• The fixtures are already in-tree. test/python/wien2k/SrVO3.oubwin has genuinely non-uniform windows (5 / 4 / 4 / 3x7 bands), which is precisely the case that makes the :nn slicing worth testing.
• The class docstring advertises the seam — "wienroot and dmftproj_exe are the only coupling to WIEN2k, so a fake WIEN2k can be substituted for testing by pointing them elsewhere" — and then nothing uses it.
• Most telling: both blockers above live in the untested surface, and 1.2 is exactly what a twenty-line stub test catches. That is how I found it.

Genuinely needs a real Wien2k: .restart semantics, .oldin2 names, the testconv criterion, _SCF_PARTS ordering, and whether lapw2 -qdmft consumes the file correctly. Everything else is fair game for CI.

2.6 SP = 1 support was removed

The previous driver.py wrote .qdmftup / .qdmftdn for SP != 0; the new one raises for n_sigma != 2 and has no SP handling. Consistent with the stated scope, but it is a removal as well as an addition — worth being explicit about in the PR description so nobody is surprised.

---

3. Consistency

3.1 run_dft_only(n_iter=N) breaks the master-only + bcast pattern

The only place in the file that does:

if not self._scf_history_on_disk():      # every rank reads case.scf
    self._prepare_fresh_scf()
    mpi.barrier(poll_msec=100)           # conditional collective

Every other decision in the driver is ... if mpi.is_master_node() else None followed by mpi.bcast. If the ranks ever disagree — a lagging shared filesystem is enough — some enter the barrier and some don't, and the job hangs. One-line fix.

3.2 Master-only raises followed by a collective

_check_inputs, _run_x, _run_checked, _read_scfm, _write_qdmft and _run_dmftproj all raise on rank 0 only; the other ranks then block in the next mpi.bcast / mpi.barrier. The VASP and QE drivers have the same shape, so this is not a new problem — but this driver has many more such sites, and a small helper that broadcasts the failure flag would pay for itself.

3.3 _converged treats a missing :DIS as charge-converged

dis_ok = dis is None or dis < self.ccut is fine for the first cycle, but scf_converged() and _ensure_converged_scf use the same predicate to decide whether to reuse an existing SCF. A three-record history whose newest record lost its :DIS — which the driver's own _read_scfm warning calls "meaningless" — reports converged, and the SCF is skipped. On the reuse path, require dis is not None.

3.4 Energy label collision across drivers

VASP and QE print DFT + DMFT Total Energy before the update, as dft_energy + band_energy_correction + Eint_m_dc. This driver prints the same label after the update, read from case.scf. Both are right for their code, but they are different quantities at different points in the cycle and share a label — worth distinguishing in the text.

Nits (click to expand)

• _read_oubwin returns iso, which is really the SO flag (the converter asserts it equals SO); every caller discards it, and it is never checked against the archive.
• _RY_IN_EV is hardcoded as 13.605698; the QE driver uses 13.605693123; the archive carries energy_unit from case.ctqmcout. Only affects read_dft_energy, but it is three values for one constant.
• _prepare_fresh_scf matches '.broyd' in name (substring) rather than a f"{seedname}.broyd" prefix.
• _run_scf raises "did not converge in {max_scf_iter} cycles", but with a seeded history the limit is max_scf_iter additional cycles, so the number in the message is not the number of cycles run.
• _env(): an unset SCRATCH makes the tcsh x script abort with a message the driver cannot explain — cheap to check in _check_inputs. The whitelist also drops MKL/OMP tuning and locale.
• README says mpirun -n 16; worth stating that those 16 ranks idle through the whole lapw chain, so users size the job knowing that.
• Example: dmft_cycle raises NameError on Gimp / Sigma_imp_iw if n_loops == 0. Guard it or document n_dmft_loops >= 1.
• Example reassigns target_density, obe = ...update... whereas the VASP example keeps target_density. The Wien2k version is arguably more correct, since the dmftproj window can shift between updates — just flagging the divergence so it is a decision rather than an accident.

---

Not verifiable without a Wien2k installation

.restart semantics, the .oldin2 / .oldin2a names, the exact testconv energy criterion, the _SCF_PARTS ordering, the _LAPW0_SAVE / _MIXER_SAVE sets, and the fort.225 / fort.77 cleanup. These all read plausibly against run_lapw, but someone with Wien2k on disk should confirm them.

[harrisonlabollita]

Thanks for this review! Everything below is in the ten commits that I just pushed. I've implemented/addressed most of your comments/suggestions and provided additional context for things that I've left off for the moment.

1. Blockers

1.1 — correlated shell on Sr. Fixed in Put the correlated shell on V in the SrVO3 example. The provided case.indmftpr and case.struct files were incompatible.

1.2 — .restart never removed. Fixed, and I took your second suggestion rather than the first: .restart is gone entirely. _flush_mixing_history now arranges for case.broyd* to be deleted immediately before the mixer call that would otherwise have inherited them, which is what run_lapw does when it starts an SCF rather than continuing one.

2. Should fix

2.1 — complex cases. Threaded through rather than rejected. _run_x prepends -c for lapw1 and lapw2 when _cmplx is set, following run_lapw; lapw0, lcore and mixer are excluded since they have no complex variant.

2.2 — restart hole in run_initial_stage. Fixed. The resumed-CSC branch of _ensure_converged_scf now runs _save_old(_LAPW0_SAVE), lapw0 and lapw1 before the projectors are rebuilt, so case.vector always belongs to the density on disk.

2.3 — nn cross-check. I've decided to not add this one. Since this driver is controlling everything, everything should be locked together and we will not hit this case.

2.4 — Eint_m_dc multiplicity. Documented rather than scaled. The _write_qdmft docstring now states the convention, the value must arrive already summed over inequivalent shells, each weighted by its multiplicity, alongside the other three conventions. Here the user is left to get this correct when they implement the dmft_step.

2.5 — tests. Agreed, but I've left this off for the moment. I'd rather add the tests in a follow-up PR than grow this one further.

2.6 — SP = 1 support. I've left this off for the moment, just to get something in. Currently, Wien2k supplies a separate bash script for spin-polarized calculations. I don't think we need to supply a separate driver just extend this, but since magnetic DFT+DMFT calculations typically let the magnetism arise from the DMFT part, we can keep the DFT step unpolarized.

3. Consistency

3.1 — conditional collective in run_dft_only. Fixed; case.scf is read on the master and broadcast like everywhere else. It also reuses that history instead of having _run_scf re-read the file.

3.2 — master-only raises before a collective. Agreed that it's real, but I'd rather not fix it here alone: the shape is shared with the VASP and QE drivers, and a broadcast-the-failure-flag helper in this one file would make this driver the odd one out for whoever reads all three next. Happy to open a separate issue for a shared helper if you agree that's the right home for it.

3.3 — missing :DIS treated as converged. Fixed, and simpler than a flag: dis_ok now requires dis is not None unconditionally.

3.4 — energy label collision. Fixed. This one now prints DFT + DMFT total energy after the charge update, so it no longer shares a string with the pre-update dft_energy + band_energy_correction + Eint_m_dc that VASP and QE report.

Nits

Taken: the .broyd substring match (folded into 1.2); the _run_scf message, which now reports both the cycles run in that call and the case.scf total; _read_oubwin no longer returns iso (the record is still consumed to advance the reader); SCRATCH is checked in _check_inputs; the env whitelist gained OMP_STACKSIZE, MKL_NUM_THREADS, LANG and LC_ALL; the README now says the Wien2k chain runs on the master rank alone, so the job should be sized for the solver; and dmft_cycle raises on n_loops < 1 instead of hitting a NameError.

Three with a caveat:

• SCRATCH warns, it does not raise. Current x_lapw guards the dereference, so an unset value is harmless on a modern installation and a hard error would break a working setup.
• _RY_IN_EV stays at 13.605698. It's WIEN2k's own constant, and everything converted here is either read from or written for WIEN2k.
• target_density reassignment in the example — thanks for flagging it as a divergence rather than a bug. It is deliberate: the dmftproj window can shift between updates, so re-reading it is the correct thing for this driver even though the VASP example doesn't.

Separately, I trimmed some of the longer explanatory comments while I was reviewing.

[the-hampel]

Thanks — I went through all ten commits and re-ran them against the stub-WIEN2k harness I used for the first review. Everything you claim is there; nothing needed walking back. One ask below, then this is good to go from my side.

Attribution: as before, the verification below was carried out by Claude (Opus 5) running in Claude Code, driven by @the-hampel.

Confirmed fixed

1.1 — correlated shell. Replayed dmftproj.f's reader against the new file: sort 2 (V) now carries 0 0 2 0 / 0 0 2 0 / 01 / 0, Sr takes the 1 0 0 0 block, O keeps its p, window unchanged. Consistent with the README and the script. One small correction to your note: only case.indmftpr changed — case.struct was already right, the two blocks were simply attached to the wrong sorts.

1.2 — mixing history. .restart is gone from the file entirely. I scope-tested the new _remove_broyden_files:

before mixer : SrVO3.broyd1, SrVO3.broyd2, SrVO3.broydENE, keepme.broyd1, other.broyd1
after  mixer : keepme.broyd1, other.broyd1

Prefix rather than substring, so the unrelated files survive — that folds in the .broyd nit properly. Also verified it is one-shot (a second _run_mixer with no flush pending deletes nothing), that case.clmsum is untouched, and that the no-broyd-files case does not raise. Taking the second option was the right call: it removes the unverifiable .restart assumption from the design altogether.

2.1 — complex cases. Traced with a non-empty case.in1c:

X lapw1 -c   |   X lapw2 -c   |   X lapw2 -c -almd   |   X lapw0   |   X lcore   |   X mixer

Correct set of programs, and _in2_file switches to case.in2c so _set_in2_mode edits the right file. The -c prepend sits ahead of the _IN2_MODE_FLAGS check but does not disturb it. One thing I cannot check without a WIEN2k here: the flag order in x lapw2 -c -qdmft. x's argument loop should be order-independent, but worth a glance at your x_lapw since run_lapw writes it the other way round.

2.2 — restart hole. The resumed-CSC trace is now lapw0, lapw1, lapw2 -almd, dmftproj, so case.vector belongs to the density on disk before the projectors are rebuilt. The mixing flush stays correctly deferred: case.broyd* is still present when run_initial_stage returns and is removed at the next mixer call.

3.1 / 3.3 / 3.4. All correct. run_dft_only reads case.scf on the master and broadcasts, and feeding that history straight into _run_scf also removes the second read — nice. Passing history=[] still behaves because the guard is is None, not falsiness. dis_ok is now unconditional, and the energy string no longer collides with the VASP/QE one.

Nits all landed as described: _read_oubwin's single caller updated, env whitelist, the SCRATCH warning, the _run_scf message, the n_loops guard and the README note.

Your pushback — all accepted

2.3, 3.2, 2.5 and 2.6 are all fine by me.

On 2.3 specifically, the new _write_qdmft paragraph is accurate — I traced it, and _regenerate_projectors really is the only producer of both case.oubwin and the archive, so they cannot drift. Since that invariant is now what makes the missing check safe, it is worth keeping that paragraph verbatim if the file is ever trimmed again.

On 3.2, agreed, and your reasoning is better than mine: a broadcast-the-failure helper in one of three drivers would make this one the odd file out. A separate issue for a shared helper is the right home — please do open it.

One ask: 707abf2 trimmed four comments that were load bearing

Most of that commit is fair — but four of the removals took out rationale rather than verbosity, and I would like them back:

• _scf_tags lost the paragraph explaining why the (ene, dis) pairing is done per record rather than by zipping two independent lists. That is the subtlest logic in the file, and the comment is precisely what stops someone later "simplifying" it to zip(enes, dises) — which silently disables the charge convergence test rather than failing. I reproduced that failure mode while writing the first review.
• _run_x and _run_checked both lost "Success is checked twice, because WIEN2k signals failure two ways…". That paragraph is the only thing justifying _error_text alongside the return-code check; without it the .error handling reads as redundant and is an obvious candidate for deletion.
• The module docstring lost "lapw2 -qdmft cannot be run in parallel at all", replaced by the TODO. That sentence is the reason the driver is serial and the single most useful line for anyone sizing a job. Please keep it alongside the TODO rather than instead of it.
• _LAPW0_SAVE now reads # case.clmsum_old is mixer's previous-iteration density with the predicate cut off — "on unit 10, so the second set is required for mixing to work at all" is gone and the sentence dangles. That one looks like an accident rather than a decision.

The rest of the trim I have no issue with.

Minor

• _write_qdmft's docstring now says lapw2 "divides it by _RY_IN_EV". lapw2 divides by its own hardcoded constant, which happens to equal ours; the previous literal stated the Fortran fact, whereas the symbol implies a coupling that does not exist.
• Three lines picked up trailing whitespace in the reworked comments (the _DEFAULT_ENV_VARS block, the _RY_IN_EV line, and the _converged docstring).

Restore those four comments and open the 3.2 issue and I am happy to approve.

[the-hampel]

Approving. Both blockers are fixed, the four comments are back, and #24 is a better-scoped home for 3.2 than what I asked for.

Attribution: as with the earlier reviews on this PR, the verification below was carried out by Claude (Opus 5) running in Claude Code, driven by @the-hampel.

Re-ran the stub-WIEN2k harness against the final state after ac643ce to confirm the comment-only commit changed no behaviour:

• SCF converges in 9 cycles; no .restart anywhere in the tree
• case.broyd* survives normal operation and is removed only by the deferred flush, scoped by prefix so unrelated *.broyd* files are left alone
• resumed-CSC trace is lapw0, lapw1, lapw2 -almd, dmftproj, so case.vector matches the density on disk before the projectors are rebuilt
• -c reaches lapw1 and lapw2 only, including lapw2 -c -almd, and _in2_file switches to case.in2c
• no trailing whitespace left in the file

On the restored comments: _scf_tags came back condensed rather than verbatim. It keeps the instruction, the :ENE/:DIS length asymmetry that makes zipping wrong, and the record-closing rule, which is enough that nobody simplifies it by accident. Fine as it stands.

Two things carried forward by agreement, neither blocking:

• Tests. Deferred to a follow-up PR. Worth tracking alongside #24 so it does not get lost — the driver has no coverage today, and the parser work on #23 is the case in point for how much a small stub-driven test finds.
• Flag order. x lapw2 -c -qdmft is the one thing neither of us can check from here — run_lapw writes it the other way round, and while x's argument loop should be order-independent, it is worth one look at your x_lapw before this sees a real run.

Nice work on the turnaround.