From a042d0958b0e698668285df3057313acce28a19f Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 27 Aug 2026 07:14:46 -0400 Subject: [PATCH 01/14] [wien2k] Add CSC DFT+DMFT driver A driver that orchestrates the chain of Wien2k programs so modest's DftDriver can drive charge self-consistency for Wien2k calculations from python. WIEN2k is a chain of programs launched by the x script rather than a single executable, and its own qdmft support inverts the control flow DftDriver needs. The driver therefore owns the SCF loop in place of run_lapw: lapw0 -> lapw1 -> lapw2 -> lcore -> mixer, tracking convergence from case.scfm. run_update_stage writes the DMFT density matrices to case.qdmft, runs lapw2 -qdmft (serial only), 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 on return. Also handles an x-script side effect the driver has to work around: a stale .oldin2 silently corrupting case.in2. Adds run_dft_only() and scf_converged() for plain DFT runs and for stepping the cycle one iteration at a time. Scope: serial, non-magnetic (SP=0, SO=0). --- python/triqs_dftkit/wien2k/__init__.py | 6 +- python/triqs_dftkit/wien2k/driver.py | 782 +++++++++++++++++++++++-- 2 files changed, 742 insertions(+), 46 deletions(-) diff --git a/python/triqs_dftkit/wien2k/__init__.py b/python/triqs_dftkit/wien2k/__init__.py index bf222a9..eb645ce 100644 --- a/python/triqs_dftkit/wien2k/__init__.py +++ b/python/triqs_dftkit/wien2k/__init__.py @@ -1,8 +1,8 @@ """ -Wien2k converter for DFT+DMFT calculations +Wien2k converter and driver for DFT+DMFT calculations """ from .converter import Converter -from .driver import Driver +from .driver import Driver, DFTWorkflowError -__all__ = ['Converter', 'Driver'] +__all__ = ['Converter', 'Driver', 'DFTWorkflowError'] diff --git a/python/triqs_dftkit/wien2k/driver.py b/python/triqs_dftkit/wien2k/driver.py index 5e6dca7..b2bd6fb 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -1,48 +1,744 @@ +""" +WIEN2k driver for TRIQS+DFT workflow automation. + +This module provides a driver class for automating WIEN2k calculations in the +context of charge self-consistent (CSC) DMFT calculations with TRIQS/modest. + +Unlike VASP (one persistent forked process) or Quantum ESPRESSO (one executable +per step), WIEN2k is a *chain* of Fortran programs -- lapw0, lapw1, lapw2, lcore, +mixer -- each launched through the ``x`` tcsh script, which writes the ``.def`` +file of unit -> filename assignments that the program reads. The SCF loop that +normally drives them lives in the ``run_lapw`` tcsh script. + +This driver replaces ``run_lapw``: it owns the loop in Python so that modest's +DftDriver can inject a DMFT charge density correction between iterations. +WIEN2k's own QDMFT support cannot be reused because it inverts the control flow +-- ``run_lapw`` calls out to a user python script from inside its cycle, whereas +DftDriver requires python to be the caller. + +Scope: serial, non-magnetic (SP=0, SO=0). ``lapw2 -qdmft`` cannot be run in +parallel at all: the entire DMFT code path is inside ``#ifndef Parallel`` in +SRC_lapw2/qdmft.F and is therefore compiled out of ``lapw2_mpi``. Independently, +its k-point counter is process-local while the density matrix array is indexed +globally, so k-point parallelism is broken too. +""" +import os, re, shutil, subprocess +from datetime import datetime + +import numpy as np import triqs.utility.mpi as mpi -from h5 import HDFArchive + from .converter import Converter +from ..converter_tools import ConverterTools + + +class DFTWorkflowError(Exception): + """Exception raised for errors in the DFT workflow.""" + pass + + +# Default environment variables to preserve for subprocess execution. WIENROOT +# and SCRATCH are WIEN2k specific: x interpolates $WIENROOT into every .def file +# it writes (e.g. the xc_funcs.h path for lapw0) and tcsh aborts outright on an +# undefined variable, so the environment cannot simply be stripped to the +# defaults the other drivers use. +_DEFAULT_ENV_VARS = ['PATH', 'LD_LIBRARY_PATH', 'SHELL', 'PWD', 'HOME', 'OMP_NUM_THREADS', + 'OMPI_MCA_btl_vader_single_copy_mechanism', 'WIENROOT', 'SCRATCH'] + +# WIEN2k works in Rydberg, modest in eV. +_RY_IN_EV = 13.605698 + +# Flags that make x rewrite the first five characters of case.in2 in place and +# restore the original from .oldin2 afterwards. +_IN2_MODE_FLAGS = ('-almd', '-qdmft', '-fermi', '-qtl', '-alm', '-efg') + +# Files saved to _old before lapw0 and before mixer, as run_lapw:471-474 +# and run_lapw:991-995 do. case.clmsum_old is mixer's previous-iteration density +# on unit 10, so the second set is required for mixing to work at all. +_LAPW0_SAVE = ('vsp', 'vns', 'r2v') +_MIXER_SAVE = ('clmsum', 'vrespsum', 'tausum') + +# Partial scf files concatenated into case.scf each cycle, in run_lapw's order +# for the non-spin-polarised non-HF branch. Missing ones are skipped. +_SCF_PARTS = ('0', '1', 'so', '2', '1s', '2s', 'c') + +# ':DIS : CHARGE DISTANCE ( for atom spin)'. The +# value inside the parentheses is the per-atom maximum, which is what testconv +# tests against the -cc limit; the trailing number is the cell total. +_DIS_RE = re.compile(r'\(\s*([-+0-9.EDed]+)\s+for atom') + class Driver(object): + """ + Driver orchestrating the WIEN2k program chain for CSC DFT+DMFT. + + Satisfies modest's DftDriver contract: it exposes ``seedname``, + ``run_initial_stage`` and ``run_update_stage``, and leaves + ``.h5`` fully re-converted before either hook returns. + + Attributes + ---------- + seedname : str + WIEN2k case name. All WIEN2k files are ``.`` and the + archive modest reads is ``.h5``. + wienroot : str + WIEN2k installation directory; ``/x`` launches every program. + dmftproj_exe : str + The dmftproj executable. Looked up on $PATH by default, matching how + run_lapw invokes it. + max_scf_iter : int + Iteration cap for the initial SCF loop. + ecut, ccut : float + Energy (Ry) and charge convergence limits, as run_lapw's -ec and -cc. + verbosity : int + 0 -- warnings only; 1 -- also per-cycle convergence and restart messages + (default); 2 -- also one line per WIEN2k program launched. Output is + written on the master rank only. Settable at any time as an attribute. + + Notes + ----- + ``wienroot`` and ``dmftproj_exe`` are the only coupling to WIEN2k, so a fake + WIEN2k can be substituted for testing by pointing them elsewhere. + """ + + def __init__(self, seedname, wienroot=None, dmftproj_exe='dmftproj', + max_scf_iter=100, ecut=1e-4, ccut=1e-4, verbosity=1): + self.seedname = seedname + self.wienroot = wienroot if wienroot is not None else os.getenv('WIENROOT') + if not self.wienroot: + raise DFTWorkflowError( + "WIENROOT is not set and no wienroot= was given; cannot locate the WIEN2k 'x' script.") + self.x_exe = os.path.join(self.wienroot, 'x') + if not os.path.isfile(self.x_exe): + raise DFTWorkflowError(f"No 'x' script at {self.x_exe} (is wienroot= correct?)") + self.dmftproj_exe = dmftproj_exe + self.max_scf_iter = max_scf_iter + self.ecut = ecut + self.ccut = ccut + self.verbosity = verbosity + self.fortran_to_replace = {'D': 'E'} + + def __repr__(self): + return (f"Wien2kDriver(seedname={self.seedname}, wienroot={self.wienroot}, " + f"dmftproj_exe={self.dmftproj_exe})") + + __str__ = __repr__ + + # ------------------------------------------------------------------ files + + def _report(self, message, level=1): + """ + Print progress on the master rank, if verbosity allows. + + Deliberately a per-driver setting rather than TRIQS's module-level + ``mpi.Verbosity_Level_Report_Max``, since turning that down to quieten a + few hundred lines of DFT bookkeeping would also quieten the solver. + """ + if level <= self.verbosity: + mpi.report(message) + + @staticmethod + def _warn(message): + """ + Report a problem that does not stop the run. + + Not gated on verbosity, and sent to stderr: these mean the calculation is + probably wrong rather than merely noisy, so silencing progress output + should not silence them too. + """ + mpi.report(f"WARNING: {message}", stderr=True) + + def _f(self, ext): + """Path of the WIEN2k file ``.``.""" + return f"{self.seedname}.{ext}" + + @property + def _cmplx(self): + """'c' when this is a complex (no inversion symmetry) case, else ''.""" + in1c = self._f('in1c') + return 'c' if os.path.isfile(in1c) and os.path.getsize(in1c) > 0 else '' + + def _in2_file(self): + return self._f('in2' + self._cmplx) + + def _set_in2_mode(self, mode): + """ + Overwrite the first five characters of case.in2 line 1 with ``mode``. + + x does this itself for -almd/-qdmft, but a crashed run can leave the + wrong keyword behind, which would silently turn a plain density pass + into a projector pass. + """ + path = self._in2_file() + with open(path) as fh: + lines = fh.readlines() + if not lines: + raise DFTWorkflowError(f"{path} is empty") + lines[0] = f"{mode:<5}" + lines[0].rstrip('\n')[5:] + '\n' + with open(path, 'w') as fh: + fh.writelines(lines) + + # -------------------------------------------------------------- execution + + def _env(self): + """Environment for WIEN2k subprocesses, with WIENROOT forced to ours.""" + env = {} + for name in _DEFAULT_ENV_VARS: + value = os.getenv(name) + if value: + env[name] = value + env['WIENROOT'] = self.wienroot + return env + + def _run_x(self, program, *flags): + """ + Run ``x `` and verify it succeeded. Master node only. + + Success is checked twice, because WIEN2k signals failure two ways: x + exits 9 when the program returns non-zero, and every program writes a + message into ``.error`` on entry and truncates it again just + before a successful exit. The error file therefore reports failures + that never reach the exit status, including untrapped runtime errors. + """ + if not mpi.is_master_node(): + return 0 + + # x backs case.in2 up to .oldin2 only when that file is absent, and + # -qdmft seds .oldin2 rather than case.in2 -- so a .oldin2 left behind by + # an earlier crash is silently used as the source and then restored over + # the real input. x warns about this but does not clean it up. + if any(flag in _IN2_MODE_FLAGS for flag in flags): + for stale in ('.oldin2', '.oldin2a'): + if os.path.isfile(stale): + os.remove(stale) + + return self._run_checked([self.x_exe, program, *flags], f'{program}.error', + f"x {program} {' '.join(flags)}".rstrip()) + + def _run_checked(self, command, error_file, label): + """ + Run a WIEN2k command and verify it succeeded. + + Success is checked twice, because WIEN2k signals failure two ways: a + non-zero exit status, and a message left in the .error file. Every program + writes that message on entry and truncates it again just before a + successful exit, so it catches failures that never reach the exit status, + including untrapped runtime errors. + """ + self._report(f"[{datetime.now()}] running {' '.join(command)}", level=2) + result = subprocess.run(command, capture_output=True, text=True, env=self._env()) + tail = f"{result.stdout[-2000:]}{result.stderr[-2000:]}" + if result.returncode != 0: + raise DFTWorkflowError(f"{label} failed with code {result.returncode}\n" + f"{self._error_text(error_file)}{tail}") + message = self._error_text(error_file) + if message: + raise DFTWorkflowError(f"{label} exited cleanly but left {error_file}:\n{message}") + return 0 + + @staticmethod + def _error_text(error_file): + """Contents of a WIEN2k .error file; empty string means success.""" + if not os.path.isfile(error_file) or os.path.getsize(error_file) == 0: + return '' + with open(error_file) as fh: + return fh.read().strip() + + def _run_dmftproj(self): + """Run dmftproj to build the projectors the converter reads.""" + if not mpi.is_master_node(): + return 0 + if not os.path.isfile(self._f('indmftpr')): + raise DFTWorkflowError( + f"{self._f('indmftpr')} not found; copy and edit " + f"{os.path.join(self.wienroot, 'SRC_templates', 'case.indmftpr')}") + return self._run_checked([self.dmftproj_exe], 'dmftproj.error', 'dmftproj') + + # ------------------------------------------------------------- scf output + + def _append_to_scf(self, parts): + """ + Append partial scf files to case.scf, as run_lapw does. + + Called twice per cycle: once for the programs before mixer, once for + case.scfm after it. + """ + if not mpi.is_master_node(): + return + with open(self._f('scf'), 'a') as out: + for part in parts: + path = self._f('scf' + part) + if os.path.isfile(path): + with open(path) as fh: + out.write(fh.read()) + + def _save_old(self, exts): + """ + Copy case. to case._old for each ext that exists. + + run_lapw does this at two points in every cycle and both are load + bearing. Most importantly mixer reads case.clmsum_old on unit 10 as the + previous-iteration density (SRC_mixer/mixer.F:881-895): without the copy + it mixes against a missing or stale density, which emits no :DIS line and + sends the SCF diverging until the linearisation energies go bad and + select.f aborts with "no energy limits found". + """ + if not mpi.is_master_node(): + return + for ext in exts: + path = self._f(ext) + if os.path.isfile(path): + shutil.copyfile(path, self._f(ext + '_old')) + + @staticmethod + def _scf_tags(path): + """ + Parse the (:ENE, :DIS) values out of a WIEN2k scf file. + + :ENE is written in three variants (**INFO****, *WARNING**, **********), so + the energy is the last whitespace token rather than a fixed column. :DIS + carries two numbers; the one inside the parentheses is the per-atom + maximum, which is what testconv compares against the -cc limit, and the + trailing one is the cell total. + """ + energies, distances = [], [] + if not os.path.isfile(path): + return energies, distances + with open(path) as fh: + for line in fh: + if line.startswith(':ENE'): + energies.append(float(line.split()[-1])) + elif line.startswith(':DIS'): + match = _DIS_RE.search(line) + distances.append( + float(match.group(1).replace('D', 'E').replace('d', 'e')) + if match else float(line.split()[-1])) + return energies, distances + + def _read_scfm(self): + """ + Return ``(ene, dis)`` from case.scfm: the total energy in Ry and the + per-atom charge distance. + + :ENE is written in three variants (**INFO****, *WARNING**, **********), + so the value is taken as the last whitespace token rather than by column. + """ + path = self._f('scfm') + if not os.path.isfile(path): + raise DFTWorkflowError(f"{path} was not written; mixer did not run") + energies, distances = self._scf_tags(path) + if not energies: + raise DFTWorkflowError(f"no :ENE line in {path}") + if not distances: + # mixer omits :DIS when it has no previous density to compare + # against, which means case.clmsum_old was missing -- the mixing is + # then meaningless even though mixer exits cleanly. + self._warn(f"no :DIS line in {path}; mixer had no previous density, " + "so the charge convergence test is being skipped") + return energies[-1], (distances[-1] if distances else None) + + def read_dft_energy(self): + """ + Total energy in eV, from the last :ENE line in case.scf. + + After a charge update this **already includes** the interaction energy: + lapw2 -qdmft folds ``correner`` into :SUM (ETOT = ETOT + correner) and + subtracts the DFT in-window band energy itself. Callers must therefore + not add ``Eint_m_dc`` again, and there is no separate band energy + correction to compute -- unlike the VASP and QE drivers. + """ + energy = None + if mpi.is_master_node(): + path = self._f('scf') + if not os.path.isfile(path): + raise DFTWorkflowError(f"{path} does not exist; no SCF has run") + energies, _ = self._scf_tags(path) + if not energies: + raise DFTWorkflowError(f"no :ENE line in {path}") + energy = energies[-1] * _RY_IN_EV + return mpi.bcast(energy) + + # -------------------------------------------------------------- scf cycle + + def _check_inputs(self, need_projectors=True): + """ + Verify the inputs are in place before launching anything. + + ``case.indmftpr`` is only consumed at the very end, by dmftproj, so + without an up-front check a missing one would not surface until a full + SCF had already run. + """ + if not mpi.is_master_node(): + return + # x takes the case name from the directory it runs in, while this driver + # takes it from seedname. If the two disagree, x reads and writes a + # different set of case.* files than the driver looks at. + cwd_case = os.path.basename(os.getcwd()) + if self.seedname != cwd_case: + self._warn(f"seedname is '{self.seedname}' but the working directory is " + f"'{cwd_case}'; the x script derives the case name from the " + "directory, so these must normally match") + + required = ['struct', 'in0', 'in1' + self._cmplx, 'in2' + self._cmplx, 'inm', 'inc'] + if need_projectors: + required.append('indmftpr') + missing = [self._f(ext) for ext in required if not os.path.isfile(self._f(ext))] + if missing: + hint = 'run init_lapw' + if self._f('indmftpr') in missing: + hint += ("; case.indmftpr is not made by init_lapw -- generate it with " + "init_dmftpr, or copy and edit " + f"{os.path.join(self.wienroot, 'SRC_templates', 'case.indmftpr')}") + raise DFTWorkflowError(f"missing WIEN2k input file(s): {', '.join(missing)}; {hint}") + if not os.path.isfile(self._f('clmsum')): + if os.path.isfile(self._f('clmsum_old')): + self._warn(f"{self._f('clmsum')} missing, recovering from clmsum_old") + with open(self._f('clmsum_old'), 'rb') as src, open(self._f('clmsum'), 'wb') as dst: + dst.write(src.read()) + else: + raise DFTWorkflowError( + f"no {self._f('clmsum')} or {self._f('clmsum_old')}, which lapw0 needs; run dstart") + + def _scf_iteration(self): + """ + One SCF cycle: lapw0 -> lapw1 -> lapw2 -> lcore -> mixer. + + Returns ``(ene, dis)`` read from case.scfm and broadcast, so that every + rank reaches the same convergence verdict and stays in lockstep. + """ + self._save_old(_LAPW0_SAVE) + for program in ('lapw0', 'lapw1', 'lapw2', 'lcore'): + self._run_x(program) + self._append_to_scf(_SCF_PARTS) + self._save_old(_MIXER_SAVE) + self._run_x('mixer') + self._append_to_scf(('m',)) + + values = self._read_scfm() if mpi.is_master_node() else None + return mpi.bcast(values) + + def _converged(self, history): + """ + Apply testconv's criterion to the (ene, dis) history. + + The energy test is the mean of the last two |dE| over three iterations, + not a single difference, and needs three points before it can fire. + """ + if len(history) < 3: + return False + e3, e2, e1 = (entry[0] for entry in history[-3:]) + ene_ok = 0.5 * (abs(e1 - e3) + abs(e1 - e2)) < self.ecut + dis = history[-1][1] + dis_ok = dis is None or dis < self.ccut + return ene_ok and dis_ok + + def _prepare_fresh_scf(self): + """ + Housekeeping before the first cycle of a new SCF. + + Only ever for a genuine fresh start. Dropping case.broyd* part way + through would discard the mixing history and stall convergence, which + matters when the cycle is being stepped one iteration at a time. + """ + if not mpi.is_master_node(): + return + self._set_in2_mode('TOT') + for name in os.listdir('.'): + if '.broyd' in name: + os.remove(name) + + def _run_scf(self, n_iter=None, history=None): + """ + Iterate the SCF cycle. + + ``n_iter=None`` runs to convergence, up to ``max_scf_iter``, and raises if + it is never reached. ``n_iter=N`` runs exactly N cycles and returns + whatever state they reached, without requiring convergence. + + ``history`` seeds the convergence test; when omitted it is recovered from + case.scf. That matters for a stepped run: the energy criterion needs + three points, which a history rebuilt from scratch on every call would + never accumulate. + """ + if history is None: + history = self._scf_history_on_disk() if mpi.is_master_node() else None + history = mpi.bcast(history) + history = list(history) + done = len(history) + limit = self.max_scf_iter if n_iter is None else n_iter + + for step in range(1, limit + 1): + history.append(self._scf_iteration()) + self._report(f" cycle {done + step}: :ENE = {history[-1][0]:.8f} Ry " + f":DIS = {history[-1][1]}") + if n_iter is None and self._converged(history): + self._report(f"SCF converged after {done + step} cycles") + return history + + if n_iter is None: + raise DFTWorkflowError( + f"SCF did not converge in {self.max_scf_iter} cycles " + f"(ecut={self.ecut} Ry, ccut={self.ccut})") + return history + + def _ensure_converged_scf(self, force_scf=False): + """ + Bring the SCF to convergence, reusing whatever case.scf already holds. + + Already converged -> nothing to run. Partially converged -> continue from + there, keeping the Broyden history, since restarting would throw away both + the completed cycles and the mixing state. Cold start, or force_scf -> + begin afresh. + """ + history = self._scf_history_on_disk() if mpi.is_master_node() else None + history = mpi.bcast(history) + + if not force_scf and self._converged(history): + self._report(f"case.scf already holds a converged SCF ({len(history)} cycles); " + "reusing it (force_scf=True to redo)") + return history + + if force_scf or not history: + self._prepare_fresh_scf() + mpi.barrier(poll_msec=100) + return self._run_scf(history=[]) + + self._report(f"continuing the SCF from {len(history)} cycles already in case.scf") + return self._run_scf(history=history) + + def scf_converged(self): + """ + Whether the SCF history in case.scf already meets the ecut/ccut criteria. + + Exposed so the cycle can be driven a step at a time:: + + while not driver.scf_converged(): + driver.run_dft_only(n_iter=1) + """ + history = self._scf_history_on_disk() if mpi.is_master_node() else None + return self._converged(mpi.bcast(history)) + + def run_dft_only(self, n_iter=None, force_scf=False): + """ + Run the DFT SCF cycle and stop: no projectors, no dmftproj, no HDF5. + + Not part of the DftDriver contract. This is for plain WIEN2k runs and for + driving the cycle under external control; case.indmftpr is not needed. + + Parameters + ---------- + n_iter : int, optional + Run exactly this many cycles and return, converged or not. Pass 1 to + advance a single step. The default, None, iterates to convergence and + raises if max_scf_iter is exhausted. + force_scf : bool, optional + Start over from scratch, discarding the accumulated case.scf history + and the Broyden files, instead of reusing a converged result or + continuing a partial one. Ignored when n_iter is given, since a fixed + number of cycles was asked for explicitly. + + Returns + ------- + list of (float, float) + The (:ENE in Ry, :DIS) history, including cycles already on disk. + """ + self._check_inputs(need_projectors=False) + + if n_iter is None: + return self._ensure_converged_scf(force_scf) + + # Stepping: prepare only on a genuine cold start, so that repeated + # single-step calls keep their mixing history and their :ENE record. + if not self._scf_history_on_disk(): + self._prepare_fresh_scf() + mpi.barrier(poll_msec=100) + return self._run_scf(n_iter=n_iter) + + def _scf_history_on_disk(self): + """(ene, dis) pairs recovered from an existing case.scf, for restart.""" + energies, distances = self._scf_tags(self._f('scf')) + # :DIS may be absent from older runs; pad so the pairs line up. + distances += [None] * (len(energies) - len(distances)) + return list(zip(energies, distances)) + + # --------------------------------------------------------- dmft interface + + def _read_oubwin(self): + """ + Read case.oubwin, written by dmftproj. + + Returns ``(iso, windows)`` where windows is a list of + ``(included, nb_bot, nb_top, weight)``, one per k-point, with 1-based + inclusive band indices. This file -- not dft_input/n_orbitals -- is what + lapw2 -qdmft cross-checks case.qdmft against, so it is the authority on + the per-k window. + """ + path = self._f('oubwin') + if not os.path.isfile(path): + raise DFTWorkflowError(f"{path} not found; dmftproj must run before the charge update") + reader = ConverterTools.read_fortran_file(self, path, self.fortran_to_replace) + try: + n_k = int(next(reader)) + iso = int(next(reader)) + windows = [] + for _ in range(n_k): + included = int(next(reader)) + if included == 1: + nb_bot, nb_top = int(next(reader)), int(next(reader)) + weight = next(reader) + windows.append((included, nb_bot, nb_top, weight)) + else: + windows.append((included, None, None, None)) + except StopIteration: + raise DFTWorkflowError(f"wien2k: reading file {path} failed!") + return iso, windows + + def _write_qdmft(self, N_k, Eint_m_dc, mu=0.0, beta=0.0): + """ + Write the DMFT band occupation matrices to case.qdmft for lapw2 -qdmft. + + The layout is fixed by the reader in SRC_lapw2/qdmft.F::readdata_qdmft:: + + mu (read, then unused) + beta (read, then unused) + per k-point: + nn must equal nb_top - nb_bot + 1 + nn records of 2*nn reals: Re Im, one record per matrix *row* + one throwaway record + correner in eV; lapw2 divides it by 13.605698 + + Three conventions differ from the VASP and QE writers: + + * the **full** occupation matrix is written, not the deviation from the + Kohn-Sham density -- lapw2 excises the in-window DFT bands itself; + * the matrix must be **unweighted**, because lapw2 multiplies it by the + k-point weight it reads from case.oubwin; + * ``correner`` is written in eV with no conversion. + + The blank line after each matrix is mandatory: the reader issues a + ``READ(32,*)`` with an empty io-list, which consumes one whole record. + """ + if not mpi.is_master_node(): + return + + _, windows = self._read_oubwin() + if len(windows) != N_k.shape[0]: + raise DFTWorkflowError( + f"case.oubwin has {len(windows)} k-points but N_k has {N_k.shape[0]}") + + excluded = [ik for ik, w in enumerate(windows) if w[0] != 1] + if excluded: + # lapw2 loops over every k-point when reading case.qdmft but only + # sets nn for included ones, so an excluded k-point is compared + # against uninitialised memory. + raise DFTWorkflowError( + f"case.oubwin marks k-point(s) {excluded} as not included; lapw2 -qdmft " + "requires every k-point to be inside the correlated window " + "(widen the energy window in case.indmftpr)") + + n_sigma = N_k.shape[1] + if n_sigma != 2: + raise DFTWorkflowError( + f"expected 2 spin channels for a non-magnetic case, got {n_sigma}; " + "spin-polarised and spin-orbit cases are not supported yet") + # Both channels are computed independently even when SP=0, so average + # them. Build a new array: modest hands over the caller's N_k uncopied. + Nk_avg = 0.5 * (N_k[:, 0] + N_k[:, 1]) + + for ik, (_, nb_bot, nb_top, _w) in enumerate(windows): + nn = nb_top - nb_bot + 1 + if nn > Nk_avg.shape[1]: + raise DFTWorkflowError( + f"k-point {ik}: case.oubwin window is {nn} bands but N_k only " + f"has {Nk_avg.shape[1]}") + + with open(self._f('qdmft'), 'w') as fh: + fh.write("%.14f\n" % mu) + fh.write("%.14f\n" % beta) + for ik, (_, nb_bot, nb_top, _w) in enumerate(windows): + nn = nb_top - nb_bot + 1 + fh.write("%s\n" % nn) + block = Nk_avg[ik, :nn, :nn] + for row in range(nn): + fh.write(''.join(f"{block[row, col].real:.14f} {block[row, col].imag:.14f} " + for col in range(nn))) + fh.write("\n") + fh.write("\n") # the mandatory throwaway record + fh.write("%.16f\n" % np.real(Eint_m_dc)) + + def _regenerate_projectors(self): + """ + Rebuild the projectors and reconvert the archive. + """ + self._run_x('lapw2', '-almd') + self._run_dmftproj() + if mpi.is_master_node(): + Converter(filename=self.seedname).convert_dft_input() + # lapw2 -almd writes one fort.225 record per band, l and m for every + # atom and k-point through a unit with no .def entry. + if os.path.isfile('fort.225'): + os.remove('fort.225') + mpi.barrier(poll_msec=100) + + # ---------------------------------------------------------------- the API + + def run_initial_stage(self, force_scf=False, **kwargs): + """ + Converge the DFT SCF cycle, then build projectors and convert to HDF5. + + Restart-aware: a converged case.scf is reused and a partial one continued + rather than recomputed, because modest calls this hook unconditionally and + keeps no DFT state of its own across a resumed run. Pass + ``force_scf=True`` to redo the SCF from scratch regardless. + """ + self._check_inputs() + self._ensure_converged_scf(force_scf) + self._regenerate_projectors() + return 0 + + def run_update_stage(self, N_k, Eint_m_dc, mu=0.0, beta=0.0, **kwargs): + """ + Apply the DMFT charge density correction and reconverge one cycle. + + Runs lapw2 -qdmft (which rebuilds the valence density from ``N_k``), + lcore and mixer to obtain the new density, then lapw0 and lapw1 to get + the potential and eigenvectors that go with it, and finally rebuilds the + projectors so the archive is current when this returns. + + Safe to call repeatedly with a fresh ``N_k``, which modest's CSC loop + does several times per DMFT iteration. + + ``mu`` and ``beta`` are accepted only to fill their slots in case.qdmft; + lapw2 reads both and uses neither. + """ + self._write_qdmft(N_k, Eint_m_dc, mu=mu, beta=beta) + mpi.barrier(poll_msec=100) + + # Serial by necessity: the DMFT path is compiled out of lapw2_mpi. + self._run_x('lapw2', '-qdmft') + self._run_x('lcore') + self._append_to_scf(_SCF_PARTS) + self._save_old(_MIXER_SAVE) + self._run_x('mixer') + self._append_to_scf(('m',)) + + self._save_old(_LAPW0_SAVE) + self._run_x('lapw0') + self._run_x('lapw1') + self._regenerate_projectors() + + if mpi.is_master_node() and os.path.isfile('fort.77'): + os.remove('fort.77') + self._report(f"DFT + DMFT Total Energy: {self.read_dft_energy()} eV") + mpi.barrier(poll_msec=100) + return 0 + + def kill(self): + """ + No-op teardown. - def __init__(self, seedname): self.seedname = seedname - - def run_initial_stage(self): - Converter(self.seedname).convert_dft_input() - return - - def run_update_stage(self, N_k, Ecorr, **kwargs): - beta = kwargs.pop('beta'); mu = kwargs.pop('mu') - self.write_charge_correction(N_k, Ecorr, beta, mu) - return - - def write_charge_correction(self, N_k, Ecorr, beta, mu): - energy_unit = 13.605698 # eV to Ry - if not mpi.is_master_node(): return - - n_k = N_k.shape[0] - with HDFArchive(f"{self.seedname}.h5", 'r') as ar: - n_bands_per_k = ar['dft_input']['n_orbitals'] - SO = ar['dft_input']['SO'] - SP = ar['dft_input']['SP'] - - def write_spin_block(file, Nksp, isp): - file.write("%.14f\n" % (mu / energy_unit) ) - file.write("%.14f\n" % (beta * energy_unit) ) - for ik in range(n_k): - file.write("%s\n" % n_bands_per_k[ik,isp]) - mat = Nksp[ik] - for inu in range(n_bands_per_k[ik,isp]): - for imu in range(n_bands_per_k[ik,isp]): - file.write(f"{mat[inu,imu].real:.14f} {mat[inu,imu].imag:.14f} ") - file.write("\n") - file.write("\n") - file.write("%.16f\n" % Ecorr.real) - - if SP == 0: - Nk_avg = 0.5 * (N_k[:, 0] + N_k[:, 1]) - with open(f"{self.seedname}.qdmft", "w") as f: write_spin_block(f, Nk_avg, 0) - else: - spin_to_data_idx = [0,1] if SO == 0 else [0,0] - with open(f"{self.seedname}.qdmftup", "w") as fup, open(f"{self.seedname}.qdmftdn", "w") as fdn: - for f, isp in zip([fup, fdn], spin_to_data_idx): write_spin_block(f, N_k[:, isp], isp) - return + WIEN2k runs as short-lived subprocesses, so there is nothing to stop. + Defined because CSC drivers are torn down with ``driver.kill()`` in a + finally block, and the VASP driver does have a process to terminate. + """ + return None From 2efbaef270aeb8d733c508db0430ce4be88a8ab0 Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 27 Aug 2026 12:07:05 -0400 Subject: [PATCH 02/14] [wien2k] csc DFT+DMFT example for SrVO3 --- doc/examples/wien2k_csc_svo/README.md | 81 +++++++ doc/examples/wien2k_csc_svo/SrVO3.indmftpr | 15 ++ doc/examples/wien2k_csc_svo/SrVO3.struct | 217 ++++++++++++++++++ .../wien2k_csc_svo/wien2k_modest_csc.py | 166 ++++++++++++++ 4 files changed, 479 insertions(+) create mode 100644 doc/examples/wien2k_csc_svo/README.md create mode 100644 doc/examples/wien2k_csc_svo/SrVO3.indmftpr create mode 100644 doc/examples/wien2k_csc_svo/SrVO3.struct create mode 100644 doc/examples/wien2k_csc_svo/wien2k_modest_csc.py diff --git a/doc/examples/wien2k_csc_svo/README.md b/doc/examples/wien2k_csc_svo/README.md new file mode 100644 index 0000000..010cc08 --- /dev/null +++ b/doc/examples/wien2k_csc_svo/README.md @@ -0,0 +1,81 @@ +# Charge self-consistent DFT+DMFT for SrVO3 with Wien2k + +A charge self-consistent DFT+DMFT calculation using Wien2k and ModEST. + +Here, we project onto the V-t2g shell in a wide energy window (-10 eV to +10 eV) around the Fermi energy. + +## What is here + +| file | role | +|---|---| +| `SrVO3.struct` | the Wien2k structure file | +| `SrVO3.indmftpr` | the dmftproj projector definition | +| `wien2k_modest_csc.py` | the CSC DFT+DMFT script | + +## Prerequisites + +* Wien2k, with `WIENROOT` set and `$WIENROOT/x` on disk +* `dmftproj` on `$PATH` (built and installed by this package) +* `triqs_modest`, and an impurity solver — the script uses `triqs_cthyb` + +## Running it + +Wien2k takes the case name from the directory name, so the directory must be +called `SrVO3` to match `seedname` in the script. + +```bash +mkdir SrVO3 && cd SrVO3 +cp /path/to/wien2k_csc_svo/SrVO3.struct . + +# generate the Wien2k inputs (in0, in1, in2, inm, inc, klist, the starting +# density, ...). Accept the defaults unless you know you want otherwise. +init_lapw -b -numk 1000 -rkmax 7.0 + +# the dmftproj input: either copy the one provided +cp /path/to/wien2k_csc_svo/SrVO3.indmftpr . +# ...or generate your own interactively +init_dmftpr + +cp /path/to/wien2k_csc_svo/wien2k_modest_csc.py . +mpirun -n 16 python wien2k_modest_csc.py +``` + +The driver converges the DFT SCF cycle itself on the first call, so there is no +separate `run_lapw` step. An already-converged `case.scf` is reused; pass +`force_scf=True` to `one_body_elements_from_dft()` to redo it. + +The run is checkpointed every iteration into a `svo_csc_beta_U_J.ckpt` +directory, so it can be killed and restarted: the DMFT state is restored from the +checkpoint and the Wien2k state from the files on disk. Restarting is just +re-running the same command. + +To check the DFT part alone first, without any DMFT: + +```python +from triqs_dftkit.wien2k import Driver +Driver("SrVO3").run_dft_only() # or run_dft_only(n_iter=1) for one cycle +``` + +## The projection window + +Line 15 of `SrVO3.indmftpr` is the correlated energy window relative to E_F, in +Rydberg. + +Every k-point must have bands inside this window: `lapw2 -qdmft` indexes its +density matrices per k-point but only sets their dimension for k-points that +dmftproj marked as included, so an excluded k-point would be read against +uninitialised memory. The driver checks `case.oubwin` and refuses up front +rather than letting that happen, so if it reports excluded k-points, widen the +window. + +## Notes + +`lapw2 -qdmft` runs serially regardless of how the rest of the calculation is +parallelised. The impurity solver still uses all available ranks. + +The driver takes a `verbosity` argument: `1` (the default) prints one line per SCF +cycle, `2` adds one line per Wien2k program launched, and `0` leaves only warnings. +Warnings ignore the setting and go to stderr, so quietening the progress output +cannot hide a problem. + +Only non-magnetic calculations (`SP=0`, `SO=0`) are supported at present. diff --git a/doc/examples/wien2k_csc_svo/SrVO3.indmftpr b/doc/examples/wien2k_csc_svo/SrVO3.indmftpr new file mode 100644 index 0000000..57228d8 --- /dev/null +++ b/doc/examples/wien2k_csc_svo/SrVO3.indmftpr @@ -0,0 +1,15 @@ +3 +1 1 3 +3 +cubic +0 0 2 0 +0 0 2 0 +01 +0 +cubic +1 0 0 0 +0 0 0 0 +cubic +0 1 0 0 +0 0 0 0 +-0.73499 0.73499 diff --git a/doc/examples/wien2k_csc_svo/SrVO3.struct b/doc/examples/wien2k_csc_svo/SrVO3.struct new file mode 100644 index 0000000..457b9fc --- /dev/null +++ b/doc/examples/wien2k_csc_svo/SrVO3.struct @@ -0,0 +1,217 @@ +Title +P LATTICE,NONEQUIV.ATOMS: 3221_Pm-3m +MODE OF CALC=RELA unit=bohr + 7.260500 7.260500 7.260500 90.000000 90.000000 90.000000 +ATOM 1: X=0.50000000 Y=0.50000000 Z=0.50000000 + MULT= 1 ISPLIT= 2 +Sr NPT= 781 R0=0.00001000 RMT= 2.50000 Z: 38.0 +LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000 + 0.0000000 1.0000000 0.0000000 + 0.0000000 0.0000000 1.0000000 +ATOM 2: X=0.00000000 Y=0.00000000 Z=0.00000000 + MULT= 1 ISPLIT= 2 +V NPT= 781 R0=0.00005000 RMT= 1.91 Z: 23.0 +LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000 + 0.0000000 1.0000000 0.0000000 + 0.0000000 0.0000000 1.0000000 +ATOM -3: X=0.50000000 Y=0.00000000 Z=0.00000000 + MULT= 3 ISPLIT=-2 + -3: X=0.00000000 Y=0.50000000 Z=0.00000000 + -3: X=0.00000000 Y=0.00000000 Z=0.50000000 +O NPT= 781 R0=0.00010000 RMT= 1.70 Z: 8.0 +LOCAL ROT MATRIX: 0.0000000 0.0000000 1.0000000 + 0.0000000 1.0000000 0.0000000 + -1.0000000 0.0000000 0.0000000 + 48 NUMBER OF SYMMETRY OPERATIONS +-1 0 0 0.00000000 + 0-1 0 0.00000000 + 0 0-1 0.00000000 + 1 +-1 0 0 0.00000000 + 0-1 0 0.00000000 + 0 0 1 0.00000000 + 2 +-1 0 0 0.00000000 + 0 0-1 0.00000000 + 0-1 0 0.00000000 + 3 +-1 0 0 0.00000000 + 0 0 1 0.00000000 + 0-1 0 0.00000000 + 4 +-1 0 0 0.00000000 + 0 0-1 0.00000000 + 0 1 0 0.00000000 + 5 +-1 0 0 0.00000000 + 0 0 1 0.00000000 + 0 1 0 0.00000000 + 6 +-1 0 0 0.00000000 + 0 1 0 0.00000000 + 0 0-1 0.00000000 + 7 +-1 0 0 0.00000000 + 0 1 0 0.00000000 + 0 0 1 0.00000000 + 8 + 0-1 0 0.00000000 +-1 0 0 0.00000000 + 0 0-1 0.00000000 + 9 + 0-1 0 0.00000000 +-1 0 0 0.00000000 + 0 0 1 0.00000000 + 10 + 0 0-1 0.00000000 +-1 0 0 0.00000000 + 0-1 0 0.00000000 + 11 + 0 0 1 0.00000000 +-1 0 0 0.00000000 + 0-1 0 0.00000000 + 12 + 0 0-1 0.00000000 +-1 0 0 0.00000000 + 0 1 0 0.00000000 + 13 + 0 0 1 0.00000000 +-1 0 0 0.00000000 + 0 1 0 0.00000000 + 14 + 0 1 0 0.00000000 +-1 0 0 0.00000000 + 0 0-1 0.00000000 + 15 + 0 1 0 0.00000000 +-1 0 0 0.00000000 + 0 0 1 0.00000000 + 16 + 0-1 0 0.00000000 + 0 0-1 0.00000000 +-1 0 0 0.00000000 + 17 + 0-1 0 0.00000000 + 0 0 1 0.00000000 +-1 0 0 0.00000000 + 18 + 0 0-1 0.00000000 + 0-1 0 0.00000000 +-1 0 0 0.00000000 + 19 + 0 0 1 0.00000000 + 0-1 0 0.00000000 +-1 0 0 0.00000000 + 20 + 0 0-1 0.00000000 + 0 1 0 0.00000000 +-1 0 0 0.00000000 + 21 + 0 0 1 0.00000000 + 0 1 0 0.00000000 +-1 0 0 0.00000000 + 22 + 0 1 0 0.00000000 + 0 0-1 0.00000000 +-1 0 0 0.00000000 + 23 + 0 1 0 0.00000000 + 0 0 1 0.00000000 +-1 0 0 0.00000000 + 24 + 0-1 0 0.00000000 + 0 0-1 0.00000000 + 1 0 0 0.00000000 + 25 + 0-1 0 0.00000000 + 0 0 1 0.00000000 + 1 0 0 0.00000000 + 26 + 0 0-1 0.00000000 + 0-1 0 0.00000000 + 1 0 0 0.00000000 + 27 + 0 0 1 0.00000000 + 0-1 0 0.00000000 + 1 0 0 0.00000000 + 28 + 0 0-1 0.00000000 + 0 1 0 0.00000000 + 1 0 0 0.00000000 + 29 + 0 0 1 0.00000000 + 0 1 0 0.00000000 + 1 0 0 0.00000000 + 30 + 0 1 0 0.00000000 + 0 0-1 0.00000000 + 1 0 0 0.00000000 + 31 + 0 1 0 0.00000000 + 0 0 1 0.00000000 + 1 0 0 0.00000000 + 32 + 0-1 0 0.00000000 + 1 0 0 0.00000000 + 0 0-1 0.00000000 + 33 + 0-1 0 0.00000000 + 1 0 0 0.00000000 + 0 0 1 0.00000000 + 34 + 0 0-1 0.00000000 + 1 0 0 0.00000000 + 0-1 0 0.00000000 + 35 + 0 0 1 0.00000000 + 1 0 0 0.00000000 + 0-1 0 0.00000000 + 36 + 0 0-1 0.00000000 + 1 0 0 0.00000000 + 0 1 0 0.00000000 + 37 + 0 0 1 0.00000000 + 1 0 0 0.00000000 + 0 1 0 0.00000000 + 38 + 0 1 0 0.00000000 + 1 0 0 0.00000000 + 0 0-1 0.00000000 + 39 + 0 1 0 0.00000000 + 1 0 0 0.00000000 + 0 0 1 0.00000000 + 40 + 1 0 0 0.00000000 + 0-1 0 0.00000000 + 0 0-1 0.00000000 + 41 + 1 0 0 0.00000000 + 0-1 0 0.00000000 + 0 0 1 0.00000000 + 42 + 1 0 0 0.00000000 + 0 0-1 0.00000000 + 0-1 0 0.00000000 + 43 + 1 0 0 0.00000000 + 0 0 1 0.00000000 + 0-1 0 0.00000000 + 44 + 1 0 0 0.00000000 + 0 0-1 0.00000000 + 0 1 0 0.00000000 + 45 + 1 0 0 0.00000000 + 0 0 1 0.00000000 + 0 1 0 0.00000000 + 46 + 1 0 0 0.00000000 + 0 1 0 0.00000000 + 0 0-1 0.00000000 + 47 + 1 0 0 0.00000000 + 0 1 0 0.00000000 + 0 0 1 0.00000000 + 48 diff --git a/doc/examples/wien2k_csc_svo/wien2k_modest_csc.py b/doc/examples/wien2k_csc_svo/wien2k_modest_csc.py new file mode 100644 index 0000000..726848c --- /dev/null +++ b/doc/examples/wien2k_csc_svo/wien2k_modest_csc.py @@ -0,0 +1,166 @@ +# ============================================================================ +# Charge self-consistent (CSC) DFT+DMFT for SrVO3 with Wien2k + TRIQS/modest +# ============================================================================ + +import numpy as np + +import triqs.utility.mpi as mpi + +from triqs.gfs import MeshImFreq + +from triqs_cthyb import solve_generic, TailFitParams + +import triqs_modest as tm +from triqs_modest.utils import Checkpointer, IterationData +from triqs_modest.dft_driver import DftDriver + +from triqs_dftkit.wien2k import Driver as Wien2kDriver + + +# --- physical and run parameters --------------------------------------------- +seedname = "SrVO3" +U, J = 4.5, 0.65 # eV +Up = U - 2 * J # rotationally invariant Kanamori +beta = 20.0 # 1/eV +n_iw = 1000 + +n_csc_loops = 50 # outer DFT charge updates +n_dmft_loops = 1 # DMFT iterations per charge update +n_dmft_loops_first = 5 # ... except the first, where Sigma is converged + # against the unmodified DFT charge density +dc_method = "cHeld" # Held's DC, the usual choice for a t2g-only shell + +mesh = MeshImFreq(beta=beta, statistic="Fermion", n_iw=n_iw) + +# --- solver setup ------------------------------------------------------------ +solver_params = dict( + n_tau=10 * n_iw, + length_cycle=500, + n_cycles=int(7e6 / mpi.size), + n_warmup_cycles=int(1e4), +) +tail_fit_params = TailFitParams(fit_min_w=20, fit_max_w=24, fit_max_moment=4) + + +def dmft_cycle(obe, target_density, Sigma_imp_dyn, Sigma_imp_hf, Sigma_dc, n_loops): + """Run `n_loops` DMFT iterations at fixed H(k), then prepare the DFT feedback. + + The run-level constants (embedding `E`, `h_int`, `dc`, `deg_blocks`, mesh and + solver parameters) are read from module scope; everything that changes from + one CSC iteration to the next is passed in explicitly. + + After the last impurity solve the double counting is recomputed from the new + Gimp and mu is re-found, so that the returned N_k is consistent with the + self-energy and DC that go with it. + + Returns + ------- + (IterationData, float, numpy.ndarray) + The iteration data (mu, self-energies, DC, Gimp/Gloc/Delta), the + impurity interaction energy minus the DC energy in eV, and the + band-basis charge density correction N_k[k, sigma, nu, nu']. + """ + epsilon_d = E.extract(tm.impurity_levels(obe))[0] + + for n in range(n_loops): + Sigma_imp_hf_m_dc = [hf - sig_dc for (hf, sig_dc) in zip(Sigma_imp_hf, Sigma_dc)] + + Sigma_C = E.embed([Sigma_imp_dyn], [Sigma_imp_hf_m_dc]) # Embed self-energy + + mu = tm.find_chemical_potential(target_density, obe, *Sigma_C, verbosity=False) # Find chemical potential + + Gloc = E.extract(tm.gloc(obe, mu, *Sigma_C))[0] # Compute Gloc + + ed = [(eps - mu * np.eye(eps.shape[0]) - sig_dc).real for (eps, sig_dc) in zip(epsilon_d, Sigma_dc)] + + Delta = tm.symmetrize(tm.hybridization(ed, Gloc, Sigma_imp_dyn, Sigma_imp_hf), deg_blocks) # Compute Delta + + res = solve_generic(Delta, ed, h_int, postprocess=tail_fit_params, **solver_params) # Solve the impurity problem + + Sigma_imp_dyn = tm.symmetrize(res.Sigma_dynamic, deg_blocks) # Update Sigma_imp + Sigma_imp_hf = tm.symmetrize(res.Sigma_HartreeFock, deg_blocks) + Sigma_imp_iw = tm.symmetrize(res.Sigma_iw, deg_blocks) + Gimp = tm.symmetrize(res.G_iw, deg_blocks) + + mpi.report(f" [dmft {n + 1}/{n_loops}] mu = {mu:.6f} " + f"n_loc = {Gloc.total_density().real:.6f} " + f"n_imp = {Gimp.total_density().real:.6f}") + + # --- double counting and the charge density correction --- + Sigma_dc = dc.dc_self_energy(Gimp) + Eint = 0.5 * np.real((Sigma_imp_iw * Gimp).total_density()) + Eint_m_dc = Eint - dc.dc_energy(Gimp) + + Sigma_C = E.embed([Sigma_imp_dyn], [[hf - sig_dc for (hf, sig_dc) in zip(Sigma_imp_hf, Sigma_dc)]]) + mu = tm.find_chemical_potential(target_density, obe, *Sigma_C, verbosity=False) + N_k = tm.charge_density_correction(obe, mu, *Sigma_C) + + it_data = IterationData(mu=mu, + Sigma_imp_list=[Sigma_imp_dyn], + Sigma_hartree_list=[Sigma_imp_hf], + Sigma_dc_list=[Sigma_dc], + Gimp_list=[Gimp], Gloc_list=[Gloc], Delta_list=[Delta]) + return it_data, Eint_m_dc, N_k + + +# --- DFT driver -------------------------------------------------------------- +# The Wien2k driver owns the lapw0 -> lapw1 -> lapw2 -> lcore -> mixer chain in +# python, replacing run_lapw, so that the charge correction can be injected +# between cycles. run_initial_stage() converges the SCF (reusing a converged +# case.scf if there is one), runs dmftproj and converts to SrVO3.h5. +driver = DftDriver(Wien2kDriver(seedname=seedname, ecut=1e-3, ccut=1e-3)) + +target_density, obe = driver.one_body_elements_from_dft() +mpi.report(f"target_density= {target_density}") +mpi.report(obe) + +# --- embedding: one impurity, three degenerate t2g orbitals ------------------ +E = tm.make_embedding(obe.C_space) +mpi.report(E.description(True)) + +# --- interaction and double counting ---------------------------------------- +h_int = tm.make_kanamori(E.sigma_names, E.imp_decomposition(0), U, Up, J) +dc = tm.DcSolver("NonPolarized", dc_method, U, J) + +# --- DFT-only pass: block structure and the initial DC ----------------------- +mu_dft = tm.find_chemical_potential(target_density, obe, beta, verbosity=False) +Gdft = E.extract(tm.gloc(mesh, obe, mu_dft))[0] +mpi.report(f"mu_dft= {mu_dft:.6f} n_dft= {Gdft.total_density().real:.6f}") + +deg_blocks = tm.analyze_degenerate_blocks(Gdft) +mpi.report(f"degenerate blocks= {deg_blocks}") + +# --- checkpoint: restart if there is something to restart from --------------- +ckpt = Checkpointer(f"svo_csc_beta{beta}_U{U}_J{J}.ckpt") + +if (prev := ckpt.restart()): + Sigma_imp_dyn, Sigma_imp_hf = prev.Sigma_imp_list[0], prev.Sigma_hartree_list[0] + Sigma_dc = prev.Sigma_dc_list[0] + mpi.report(f"restarting from checkpoint at iteration {len(ckpt)}") +else: + Sigma_dc = dc.dc_self_energy(Gdft) + Sigma_imp_dyn, Sigma_imp_hf = E.make_zero_imp_self_energies(mesh)[0] + for ibl in range(len(Sigma_imp_hf)): + Sigma_imp_hf[ibl] += Sigma_dc[ibl] + +# --- DFT + DMFT loop --------------------------------------------------------- +for it in range(len(ckpt), n_csc_loops): + mpi.report(f"\n=== DFT+DMFT iteration {it + 1}/{n_csc_loops} ===") + + n_loops = n_dmft_loops_first if it == 0 else n_dmft_loops + it_data, Eint_m_dc, N_k = dmft_cycle(obe, target_density, Sigma_imp_dyn, Sigma_imp_hf, Sigma_dc, n_loops) + + Sigma_imp_dyn = it_data.Sigma_imp_list[0] + Sigma_imp_hf = it_data.Sigma_hartree_list[0] + Sigma_dc = it_data.Sigma_dc_list[0] + + ckpt.append(it_data, Eint_m_dc=Eint_m_dc) + mpi.report(f"[csc {it + 1}] mu = {it_data.mu:.6f} Eint-Edc = {Eint_m_dc:.6f} eV") + + # Wien2k charge update. Skipped on the last iteration so the calculation + # ends on a DMFT step. lapw2 -qdmft folds Eint-Edc into the total energy + # itself, so nothing further has to be added on this side. + if it < n_csc_loops - 1: + target_density, obe = driver.update_one_body_elements_with_charge_correction(N_k, Eint_m_dc, mu=it_data.mu, beta=beta) + +ckpt.summarize() From 9aacd3d46b26571cbd4baeaa452f4076aae95820 Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Tue, 1 Sep 2026 10:16:37 -0400 Subject: [PATCH 03/14] [wien2k] harden driver including diff restart scenarios --- python/triqs_dftkit/wien2k/driver.py | 197 +++++++++++++++++++-------- 1 file changed, 142 insertions(+), 55 deletions(-) diff --git a/python/triqs_dftkit/wien2k/driver.py b/python/triqs_dftkit/wien2k/driver.py index b2bd6fb..19d341c 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -4,8 +4,7 @@ This module provides a driver class for automating WIEN2k calculations in the context of charge self-consistent (CSC) DMFT calculations with TRIQS/modest. -Unlike VASP (one persistent forked process) or Quantum ESPRESSO (one executable -per step), WIEN2k is a *chain* of Fortran programs -- lapw0, lapw1, lapw2, lcore, +WIEN2k is a *chain* of Fortran programs -- lapw0, lapw1, lapw2, lcore, mixer -- each launched through the ``x`` tcsh script, which writes the ``.def`` file of unit -> filename assignments that the program reads. The SCF loop that normally drives them lives in the ``run_lapw`` tcsh script. @@ -17,10 +16,7 @@ DftDriver requires python to be the caller. Scope: serial, non-magnetic (SP=0, SO=0). ``lapw2 -qdmft`` cannot be run in -parallel at all: the entire DMFT code path is inside ``#ifndef Parallel`` in -SRC_lapw2/qdmft.F and is therefore compiled out of ``lapw2_mpi``. Independently, -its k-point counter is process-local while the density matrix array is indexed -globally, so k-point parallelism is broken too. +parallel at all. """ import os, re, shutil, subprocess from datetime import datetime @@ -45,15 +41,15 @@ class DFTWorkflowError(Exception): _DEFAULT_ENV_VARS = ['PATH', 'LD_LIBRARY_PATH', 'SHELL', 'PWD', 'HOME', 'OMP_NUM_THREADS', 'OMPI_MCA_btl_vader_single_copy_mechanism', 'WIENROOT', 'SCRATCH'] -# WIEN2k works in Rydberg, modest in eV. +# WIEN2k works in Rydberg, TRIQS in eV. _RY_IN_EV = 13.605698 # Flags that make x rewrite the first five characters of case.in2 in place and # restore the original from .oldin2 afterwards. _IN2_MODE_FLAGS = ('-almd', '-qdmft', '-fermi', '-qtl', '-alm', '-efg') -# Files saved to _old before lapw0 and before mixer, as run_lapw:471-474 -# and run_lapw:991-995 do. case.clmsum_old is mixer's previous-iteration density +# Files saved to _old before lapw0 and before mixer as run_lapw. +# case.clmsum_old is mixer's previous-iteration density # on unit 10, so the second set is required for mixing to work at all. _LAPW0_SAVE = ('vsp', 'vns', 'r2v') _MIXER_SAVE = ('clmsum', 'vrespsum', 'tausum') @@ -67,6 +63,9 @@ class DFTWorkflowError(Exception): # tests against the -cc limit; the trailing number is the cell total. _DIS_RE = re.compile(r'\(\s*([-+0-9.EDed]+)\s+for atom') +# DMFT label marker in the : LABEL style of run_lapw +_DMFT_MARKER = ':QDMFT: CHARGE DENSITY UPDATED FROM DMFT OCCUPATIONS' + class Driver(object): """ @@ -117,6 +116,9 @@ def __init__(self, seedname, wienroot=None, dmftproj_exe='dmftproj', self.ccut = ccut self.verbosity = verbosity self.fortran_to_replace = {'D': 'E'} + # Set when this process inherits an SCF it did not run, and cleared by the + # next mixer call. See _flush_mixing_history. + self._mixing_flush_pending = False def __repr__(self): return (f"Wien2kDriver(seedname={self.seedname}, wienroot={self.wienroot}, " @@ -251,11 +253,49 @@ def _run_dmftproj(self): if not os.path.isfile(self._f('indmftpr')): raise DFTWorkflowError( f"{self._f('indmftpr')} not found; copy and edit " - f"{os.path.join(self.wienroot, 'SRC_templates', 'case.indmftpr')}") + f"{os.path.join(self.wienroot, 'SRC_templates', 'case.indmftpr')}" + "or use the command line tool init_dmftpr") return self._run_checked([self.dmftproj_exe], 'dmftproj.error', 'dmftproj') # ------------------------------------------------------------- scf output + def _mark_dmft_cycle(self): + """Record in case.scf that the cycle starting here is a DMFT charge update.""" + if not mpi.is_master_node(): + return + with open(self._f('scf'), 'a') as out: + out.write(f"\n{_DMFT_MARKER}\n\n") + + def _run_mixer(self): + """ + Run mixer, first flushing the mixing history if this process inherited one. + + The flush is deferred to here rather than done where it is decided, so + that ``.restart`` is only ever written immediately before the mixer call + that consumes it. Writing it earlier would leave a file that some later, + unrelated mixer call would pick up and act on. + """ + if self._mixing_flush_pending: + if mpi.is_master_node(): + open('.restart', 'w').close() + self._mixing_flush_pending = False + mpi.barrier(poll_msec=100) + self._run_x('mixer') + + def _flush_mixing_history(self, reason): + """ + Arrange for the next mixer call to restart its mixing, and say so. + + Called when this process is about to mix against a case.broyd* history it + did not build. The Broyden vectors encode a sequence of density updates, + so continuing against a history whose density had a different character -- + a DMFT-corrected one, or one from an interrupted run -- mixes the current + residual into predictions derived from unrelated ones. + """ + self._warn(f"{reason}; flushing the Broyden mixing history via .restart, " + "so the next mixer call restarts its mixing with a reduced DMIX") + self._mixing_flush_pending = True + def _append_to_scf(self, parts): """ Append partial scf files to case.scf, as run_lapw does. @@ -278,10 +318,9 @@ def _save_old(self, exts): run_lapw does this at two points in every cycle and both are load bearing. Most importantly mixer reads case.clmsum_old on unit 10 as the - previous-iteration density (SRC_mixer/mixer.F:881-895): without the copy - it mixes against a missing or stale density, which emits no :DIS line and - sends the SCF diverging until the linearisation energies go bad and - select.f aborts with "no energy limits found". + previous-iteration density: without the copy it mixes against a missing + or stale density, which emits no :DIS line and sends the SCF diverging until + the linearisation energies go bad. """ if not mpi.is_master_node(): return @@ -291,29 +330,42 @@ def _save_old(self, exts): shutil.copyfile(path, self._f(ext + '_old')) @staticmethod - def _scf_tags(path): + def _scf_tags(path, stop_at_dmft=False): """ - Parse the (:ENE, :DIS) values out of a WIEN2k scf file. + Parse a WIEN2k scf file into one ``(ene, dis)`` pair per mixer record. :ENE is written in three variants (**INFO****, *WARNING**, **********), so the energy is the last whitespace token rather than a fixed column. :DIS carries two numbers; the one inside the parentheses is the per-atom maximum, which is what testconv compares against the -cc limit, and the trailing one is the cell total. + + Pairing is done by record rather than by zipping two independent lists. + mixer writes :DIS only inside, i.e. only when it could read case.clmsum_old, + but writes :ENE unconditionally, so a cycle with no previous density contributes + an :ENE with no :DIS. Zipping and padding the shorter list at the end would attribute + every later :DIS to the cycle before its own and leave the newest cycle + with ``dis=None``, which silently disables the charge convergence test. + Within one record :DIS precedes :ENE, so :ENE closes the record. + + ``stop_at_dmft`` stops at the first _DMFT_MARKER, leaving the caller only + the plain-DFT part of the history. """ - energies, distances = [], [] + history, dis = [], None if not os.path.isfile(path): - return energies, distances + return history with open(path) as fh: for line in fh: + if stop_at_dmft and line.startswith(_DMFT_MARKER): + break if line.startswith(':ENE'): - energies.append(float(line.split()[-1])) + history.append((float(line.split()[-1]), dis)) + dis = None elif line.startswith(':DIS'): match = _DIS_RE.search(line) - distances.append( - float(match.group(1).replace('D', 'E').replace('d', 'e')) - if match else float(line.split()[-1])) - return energies, distances + dis = (float(match.group(1).replace('D', 'E').replace('d', 'e')) + if match else float(line.split()[-1])) + return history def _read_scfm(self): """ @@ -326,16 +378,17 @@ def _read_scfm(self): path = self._f('scfm') if not os.path.isfile(path): raise DFTWorkflowError(f"{path} was not written; mixer did not run") - energies, distances = self._scf_tags(path) - if not energies: + history = self._scf_tags(path) + if not history: raise DFTWorkflowError(f"no :ENE line in {path}") - if not distances: + ene, dis = history[-1] + if dis is None: # mixer omits :DIS when it has no previous density to compare # against, which means case.clmsum_old was missing -- the mixing is # then meaningless even though mixer exits cleanly. self._warn(f"no :DIS line in {path}; mixer had no previous density, " "so the charge convergence test is being skipped") - return energies[-1], (distances[-1] if distances else None) + return ene, dis def read_dft_energy(self): """ @@ -345,17 +398,19 @@ def read_dft_energy(self): lapw2 -qdmft folds ``correner`` into :SUM (ETOT = ETOT + correner) and subtracts the DFT in-window band energy itself. Callers must therefore not add ``Eint_m_dc`` again, and there is no separate band energy - correction to compute -- unlike the VASP and QE drivers. + correction to compute. """ energy = None if mpi.is_master_node(): path = self._f('scf') if not os.path.isfile(path): raise DFTWorkflowError(f"{path} does not exist; no SCF has run") - energies, _ = self._scf_tags(path) - if not energies: + # Deliberately not stop_at_dmft: this must report the current total + # energy, which after a charge update is the DMFT one. + history = self._scf_tags(path) + if not history: raise DFTWorkflowError(f"no :ENE line in {path}") - energy = energies[-1] * _RY_IN_EV + energy = history[-1][0] * _RY_IN_EV return mpi.bcast(energy) # -------------------------------------------------------------- scf cycle @@ -406,12 +461,15 @@ def _scf_iteration(self): Returns ``(ene, dis)`` read from case.scfm and broadcast, so that every rank reaches the same convergence verdict and stays in lockstep. """ + if mpi.is_master_node(): + self._set_in2_mode('TOT') + self._save_old(_LAPW0_SAVE) for program in ('lapw0', 'lapw1', 'lapw2', 'lcore'): self._run_x(program) self._append_to_scf(_SCF_PARTS) self._save_old(_MIXER_SAVE) - self._run_x('mixer') + self._run_mixer() self._append_to_scf(('m',)) values = self._read_scfm() if mpi.is_master_node() else None @@ -446,6 +504,14 @@ def _prepare_fresh_scf(self): for name in os.listdir('.'): if '.broyd' in name: os.remove(name) + # Truncate case.scf too. Without this the abandoned cycles stay on disk + # and _scf_history_on_disk splices them onto the new ones, so a later + # restart applies the convergence test across the seam between two + # unrelated runs -- and the caller was promised the history was discarded. + if os.path.isfile(self._f('scf')): + self._report(f"discarding the {len(self._scf_history_on_disk())} cycles " + f"already in {self._f('scf')}") + open(self._f('scf'), 'w').close() def _run_scf(self, n_iter=None, history=None): """ @@ -496,6 +562,15 @@ def _ensure_converged_scf(self, force_scf=False): if not force_scf and self._converged(history): self._report(f"case.scf already holds a converged SCF ({len(history)} cycles); " "reusing it (force_scf=True to redo)") + resumed = self._has_dmft_cycles() if mpi.is_master_node() else None + if mpi.bcast(resumed): + # A CSC run being resumed: case.clmsum is the DMFT-corrected + # density from the interrupted run, and the mixing history behind + # it belongs to that run. Nothing here reruns the SCF -- the + # density is kept as it stands -- but the next charge update will + # mix against that inherited history. + self._flush_mixing_history( + "case.scf records DMFT charge updates from an earlier run") return history if force_scf or not history: @@ -504,8 +579,18 @@ def _ensure_converged_scf(self, force_scf=False): return self._run_scf(history=[]) self._report(f"continuing the SCF from {len(history)} cycles already in case.scf") + self._flush_mixing_history( + f"resuming an SCF from {len(history)} cycles this process did not run") return self._run_scf(history=history) + def _has_dmft_cycles(self): + """Whether case.scf records a DMFT charge update from an earlier run.""" + path = self._f('scf') + if not os.path.isfile(path): + return False + with open(path) as fh: + return any(line.startswith(_DMFT_MARKER) for line in fh) + def scf_converged(self): """ Whether the SCF history in case.scf already meets the ecut/ccut criteria. @@ -555,11 +640,14 @@ def run_dft_only(self, n_iter=None, force_scf=False): return self._run_scf(n_iter=n_iter) def _scf_history_on_disk(self): - """(ene, dis) pairs recovered from an existing case.scf, for restart.""" - energies, distances = self._scf_tags(self._f('scf')) - # :DIS may be absent from older runs; pad so the pairs line up. - distances += [None] * (len(energies) - len(distances)) - return list(zip(energies, distances)) + """ + (ene, dis) pairs recovered from an existing case.scf, for restart. + + Stops at the first DMFT marker: the cycles after it were produced by + run_update_stage against a DMFT-corrected density, and the DFT ecut/ccut + limits say nothing useful about them. + """ + return self._scf_tags(self._f('scf'), stop_at_dmft=True) # --------------------------------------------------------- dmft interface @@ -599,8 +687,8 @@ def _write_qdmft(self, N_k, Eint_m_dc, mu=0.0, beta=0.0): The layout is fixed by the reader in SRC_lapw2/qdmft.F::readdata_qdmft:: - mu (read, then unused) - beta (read, then unused) + mu in Ry (read, then unused) + beta in Ry^-1 (read, then unused) per k-point: nn must equal nb_top - nb_bot + 1 nn records of 2*nn reals: Re Im, one record per matrix *row* @@ -653,8 +741,11 @@ def _write_qdmft(self, N_k, Eint_m_dc, mu=0.0, beta=0.0): f"has {Nk_avg.shape[1]}") with open(self._f('qdmft'), 'w') as fh: - fh.write("%.14f\n" % mu) - fh.write("%.14f\n" % beta) + # Converted to Ry / Ry^-1 for the file even though lapw2 discards + # both, so the values match what the format documents and what + # legacy dft_tools wrote. + fh.write("%.14f\n" % (mu / _RY_IN_EV)) + fh.write("%.14f\n" % (beta * _RY_IN_EV)) for ik, (_, nb_bot, nb_top, _w) in enumerate(windows): nn = nb_top - nb_bot + 1 fh.write("%s\n" % nn) @@ -708,8 +799,9 @@ def run_update_stage(self, N_k, Eint_m_dc, mu=0.0, beta=0.0, **kwargs): Safe to call repeatedly with a fresh ``N_k``, which modest's CSC loop does several times per DMFT iteration. - ``mu`` and ``beta`` are accepted only to fill their slots in case.qdmft; - lapw2 reads both and uses neither. + ``mu`` (eV) and ``beta`` (1/eV) are accepted only to fill their slots in + case.qdmft, converted to Ry there. lapw2 reads both and uses neither, so + omitting them changes nothing. """ self._write_qdmft(N_k, Eint_m_dc, mu=mu, beta=beta) mpi.barrier(poll_msec=100) @@ -717,9 +809,10 @@ def run_update_stage(self, N_k, Eint_m_dc, mu=0.0, beta=0.0, **kwargs): # Serial by necessity: the DMFT path is compiled out of lapw2_mpi. self._run_x('lapw2', '-qdmft') self._run_x('lcore') + self._mark_dmft_cycle() self._append_to_scf(_SCF_PARTS) self._save_old(_MIXER_SAVE) - self._run_x('mixer') + self._run_mixer() self._append_to_scf(('m',)) self._save_old(_LAPW0_SAVE) @@ -729,16 +822,10 @@ def run_update_stage(self, N_k, Eint_m_dc, mu=0.0, beta=0.0, **kwargs): if mpi.is_master_node() and os.path.isfile('fort.77'): os.remove('fort.77') - self._report(f"DFT + DMFT Total Energy: {self.read_dft_energy()} eV") + # Four decimals, not the full float repr: a single charge update carries + # the solver's statistical error through N_k, so consecutive updates + # scatter by a few meV even after the density has settled. Printing 17 + # significant digits advertises a precision the number does not have. + self._report(f"DFT + DMFT Total Energy: {self.read_dft_energy():.4f} eV") mpi.barrier(poll_msec=100) return 0 - - def kill(self): - """ - No-op teardown. - - WIEN2k runs as short-lived subprocesses, so there is nothing to stop. - Defined because CSC drivers are torn down with ``driver.kill()`` in a - finally block, and the VASP driver does have a process to terminate. - """ - return None From 37bb53b932783ab3b90aea3849989fc050d23ec4 Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 10 Sep 2026 12:48:36 -0400 Subject: [PATCH 04/14] [wien2k] Put the correlated shell on V in the SrVO3 example --- doc/examples/wien2k_csc_svo/SrVO3.indmftpr | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/examples/wien2k_csc_svo/SrVO3.indmftpr b/doc/examples/wien2k_csc_svo/SrVO3.indmftpr index 57228d8..4d0ed1d 100644 --- a/doc/examples/wien2k_csc_svo/SrVO3.indmftpr +++ b/doc/examples/wien2k_csc_svo/SrVO3.indmftpr @@ -2,14 +2,14 @@ 1 1 3 3 cubic +1 0 0 0 +0 0 0 0 +cubic 0 0 2 0 0 0 2 0 01 0 cubic -1 0 0 0 -0 0 0 0 -cubic 0 1 0 0 0 0 0 0 -0.73499 0.73499 From 0707fcb5688acf7eaf4df1a2f4dbadf4960ead1c Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 10 Sep 2026 13:54:58 -0400 Subject: [PATCH 05/14] [wien2k] Remove case.broyd* to discard the mixing history --- python/triqs_dftkit/wien2k/driver.py | 29 +++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/python/triqs_dftkit/wien2k/driver.py b/python/triqs_dftkit/wien2k/driver.py index 19d341c..22ff491 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -266,18 +266,31 @@ def _mark_dmft_cycle(self): with open(self._f('scf'), 'a') as out: out.write(f"\n{_DMFT_MARKER}\n\n") + def _remove_broyden_files(self): + """ + Delete case.broyd*, which is how run_lapw discards a mixing history. + + run_lapw removes them whenever it starts an SCF rather than continuing + one, unless -NI is given. Deleting the files is unambiguous: mixer + rebuilds the history from scratch because there is nothing to read. + """ + prefix = self._f('broyd') + for name in os.listdir('.'): + if name.startswith(prefix): + os.remove(name) + def _run_mixer(self): """ Run mixer, first flushing the mixing history if this process inherited one. The flush is deferred to here rather than done where it is decided, so - that ``.restart`` is only ever written immediately before the mixer call - that consumes it. Writing it earlier would leave a file that some later, - unrelated mixer call would pick up and act on. + that case.broyd* is only ever removed immediately before the mixer call + that would otherwise have read it, leaving the history in place for every + other mixer call in the run. """ if self._mixing_flush_pending: if mpi.is_master_node(): - open('.restart', 'w').close() + self._remove_broyden_files() self._mixing_flush_pending = False mpi.barrier(poll_msec=100) self._run_x('mixer') @@ -292,8 +305,8 @@ def _flush_mixing_history(self, reason): a DMFT-corrected one, or one from an interrupted run -- mixes the current residual into predictions derived from unrelated ones. """ - self._warn(f"{reason}; flushing the Broyden mixing history via .restart, " - "so the next mixer call restarts its mixing with a reduced DMIX") + self._warn(f"{reason}; removing case.broyd* as run_lapw does, so the next " + "mixer call rebuilds its mixing history from scratch") self._mixing_flush_pending = True def _append_to_scf(self, parts): @@ -501,9 +514,7 @@ def _prepare_fresh_scf(self): if not mpi.is_master_node(): return self._set_in2_mode('TOT') - for name in os.listdir('.'): - if '.broyd' in name: - os.remove(name) + self._remove_broyden_files() # Truncate case.scf too. Without this the abandoned cycles stay on disk # and _scf_history_on_disk splices them onto the new ones, so a later # restart applies the convergence test across the seam between two From 3744972b036dee4b65879fb8c67fb2b72c7a2cbd Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 10 Sep 2026 13:56:41 -0400 Subject: [PATCH 06/14] [wien2k] Pass -c to lapw1 and lapw2 for complex cases --- python/triqs_dftkit/wien2k/driver.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/python/triqs_dftkit/wien2k/driver.py b/python/triqs_dftkit/wien2k/driver.py index 22ff491..c45314d 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -48,6 +48,11 @@ class DFTWorkflowError(Exception): # restore the original from .oldin2 afterwards. _IN2_MODE_FLAGS = ('-almd', '-qdmft', '-fermi', '-qtl', '-alm', '-efg') +# Programs that run_lapw calls with -c on a complex (no inversion symmetry) case, +# which makes x launch lapw1c/lapw2c and read case.in1c/case.in2c. lapw0, lcore +# and mixer have no complex variant and reject the flag. +_CMPLX_PROGRAMS = ('lapw1', 'lapw2') + # Files saved to _old before lapw0 and before mixer as run_lapw. # case.clmsum_old is mixer's previous-iteration density # on unit 10, so the second set is required for mixing to work at all. @@ -205,6 +210,10 @@ def _run_x(self, program, *flags): if not mpi.is_master_node(): return 0 + # x does not detect a complex case on its own; run_lapw passes -c. + if self._cmplx and program in _CMPLX_PROGRAMS: + flags = ('-c',) + flags + # x backs case.in2 up to .oldin2 only when that file is absent, and # -qdmft seds .oldin2 rather than case.in2 -- so a .oldin2 left behind by # an earlier crash is silently used as the source and then restored over From 929859a296790c9c0dc3bfe805b079fe4710fb26 Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 10 Sep 2026 13:56:42 -0400 Subject: [PATCH 07/14] [wien2k] Redo lapw0 and lapw1 when resuming a CSC run --- python/triqs_dftkit/wien2k/driver.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/python/triqs_dftkit/wien2k/driver.py b/python/triqs_dftkit/wien2k/driver.py index c45314d..ce59a7a 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -591,6 +591,18 @@ def _ensure_converged_scf(self, force_scf=False): # mix against that inherited history. self._flush_mixing_history( "case.scf records DMFT charge updates from an earlier run") + # run_update_stage ends with mixer, lapw0, lapw1, so an + # interrupted one can leave case.vector belonging to the + # pre-update potential while case.clmsum is already the + # post-update density. The projectors built from that vector + # would then be inconsistent with the density on disk, silently. + # Redoing lapw0 and lapw1 is cheap next to one solver call. + self._report("redoing lapw0 and lapw1 so the eigenvectors match " + "the density left by the interrupted run") + self._save_old(_LAPW0_SAVE) + self._run_x('lapw0') + self._run_x('lapw1') + mpi.barrier(poll_msec=100) return history if force_scf or not history: From c032eceab20442eaef3763073e514d8d89249796 Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 10 Sep 2026 13:56:42 -0400 Subject: [PATCH 08/14] [wien2k] Require a :DIS record before reusing an SCF --- python/triqs_dftkit/wien2k/driver.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/python/triqs_dftkit/wien2k/driver.py b/python/triqs_dftkit/wien2k/driver.py index ce59a7a..db0c549 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -502,14 +502,19 @@ def _converged(self, history): Apply testconv's criterion to the (ene, dis) history. The energy test is the mean of the last two |dE| over three iterations, - not a single difference, and needs three points before it can fire. + not a single difference, and needs three points. + + A newest record with no :DIS counts as not converged. mixer omits :DIS + only when it had no case.clmsum_old to compare against, which makes that + cycle's mixing meaningless (see _read_scfm) -- so there is no case in + which the absence of the line should let the SCF stop or be reused. """ if len(history) < 3: return False e3, e2, e1 = (entry[0] for entry in history[-3:]) ene_ok = 0.5 * (abs(e1 - e3) + abs(e1 - e2)) < self.ecut dis = history[-1][1] - dis_ok = dis is None or dis < self.ccut + dis_ok = dis is not None and dis < self.ccut return ene_ok and dis_ok def _prepare_fresh_scf(self): From 1a4ae8a73d489f7640f36345a86ae479e72296ef Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 10 Sep 2026 13:56:42 -0400 Subject: [PATCH 09/14] [wien2k] Read case.scf on the master rank in run_dft_only --- python/triqs_dftkit/wien2k/driver.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/triqs_dftkit/wien2k/driver.py b/python/triqs_dftkit/wien2k/driver.py index db0c549..f109779 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -669,12 +669,12 @@ def run_dft_only(self, n_iter=None, force_scf=False): if n_iter is None: return self._ensure_converged_scf(force_scf) - # Stepping: prepare only on a genuine cold start, so that repeated - # single-step calls keep their mixing history and their :ENE record. - if not self._scf_history_on_disk(): + history = self._scf_history_on_disk() if mpi.is_master_node() else None + history = mpi.bcast(history) + if not history: self._prepare_fresh_scf() mpi.barrier(poll_msec=100) - return self._run_scf(n_iter=n_iter) + return self._run_scf(n_iter=n_iter, history=history) def _scf_history_on_disk(self): """ From a9babf74e521b4db818ae519df2566601758b0c1 Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 10 Sep 2026 14:05:24 -0400 Subject: [PATCH 10/14] [wien2k] Address review nits in the driver --- python/triqs_dftkit/wien2k/driver.py | 52 +++++++++++++++++++--------- 1 file changed, 36 insertions(+), 16 deletions(-) diff --git a/python/triqs_dftkit/wien2k/driver.py b/python/triqs_dftkit/wien2k/driver.py index f109779..14cbed0 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -39,6 +39,7 @@ class DFTWorkflowError(Exception): # undefined variable, so the environment cannot simply be stripped to the # defaults the other drivers use. _DEFAULT_ENV_VARS = ['PATH', 'LD_LIBRARY_PATH', 'SHELL', 'PWD', 'HOME', 'OMP_NUM_THREADS', + 'OMP_STACKSIZE', 'MKL_NUM_THREADS', 'LANG', 'LC_ALL', 'OMPI_MCA_btl_vader_single_copy_mechanism', 'WIENROOT', 'SCRATCH'] # WIEN2k works in Rydberg, TRIQS in eV. @@ -447,6 +448,12 @@ def _check_inputs(self, need_projectors=True): """ if not mpi.is_master_node(): return + if not os.getenv('SCRATCH'): + self._warn("SCRATCH is not set; x is a tcsh script and older WIEN2k " + "versions dereference $SCRATCH unguarded, which aborts with " + "an error this driver cannot attribute. siteconfig normally " + "sets it to ./") + # x takes the case name from the directory it runs in, while this driver # takes it from seedname. If the two disagree, x reads and writes a # different set of case.* files than the driver looks at. @@ -568,8 +575,10 @@ def _run_scf(self, n_iter=None, history=None): if n_iter is None: raise DFTWorkflowError( - f"SCF did not converge in {self.max_scf_iter} cycles " - f"(ecut={self.ecut} Ry, ccut={self.ccut})") + f"SCF did not converge in the {limit} cycles run here " + f"({done + limit} in case.scf in total; ecut={self.ecut} Ry, " + f"ccut={self.ccut}); raise max_scf_iter and rerun to continue from " + "where this stopped") return history def _ensure_converged_scf(self, force_scf=False): @@ -692,11 +701,14 @@ def _read_oubwin(self): """ Read case.oubwin, written by dmftproj. - Returns ``(iso, windows)`` where windows is a list of - ``(included, nb_bot, nb_top, weight)``, one per k-point, with 1-based - inclusive band indices. This file -- not dft_input/n_orbitals -- is what - lapw2 -qdmft cross-checks case.qdmft against, so it is the authority on - the per-k window. + Returns the windows as a list of ``(included, nb_bot, nb_top, weight)``, + one per k-point, with 1-based inclusive band indices. This file -- not + dft_input/n_orbitals -- is what lapw2 -qdmft cross-checks case.qdmft + against, so it is the authority on the per-k window. + + The SO flag on the second record is read to advance past it and then + dropped: the converter already asserts it against case.ctqmcout, and + nothing here has any use for it. """ path = self._f('oubwin') if not os.path.isfile(path): @@ -704,7 +716,8 @@ def _read_oubwin(self): reader = ConverterTools.read_fortran_file(self, path, self.fortran_to_replace) try: n_k = int(next(reader)) - iso = int(next(reader)) + next(reader) # the SO flag + windows = [] for _ in range(n_k): included = int(next(reader)) @@ -716,7 +729,7 @@ def _read_oubwin(self): windows.append((included, None, None, None)) except StopIteration: raise DFTWorkflowError(f"wien2k: reading file {path} failed!") - return iso, windows + return windows def _write_qdmft(self, N_k, Eint_m_dc, mu=0.0, beta=0.0): """ @@ -730,7 +743,7 @@ def _write_qdmft(self, N_k, Eint_m_dc, mu=0.0, beta=0.0): nn must equal nb_top - nb_bot + 1 nn records of 2*nn reals: Re Im, one record per matrix *row* one throwaway record - correner in eV; lapw2 divides it by 13.605698 + correner in eV; lapw2 divides it by _RY_IN_EV Three conventions differ from the VASP and QE writers: @@ -742,11 +755,21 @@ def _write_qdmft(self, N_k, Eint_m_dc, mu=0.0, beta=0.0): The blank line after each matrix is mandatory: the reader issues a ``READ(32,*)`` with an empty io-list, which consumes one whole record. + + ``Eint_m_dc`` is written as given and must therefore already be the value + for the whole cell: summed over the inequivalent correlated shells, each + multiplied by the multiplicity of its equivalent atoms. That is the + caller's responsibility. + + The per-k window comes from case.oubwin while ``N_k`` is indexed by the + archive's n_orbitals, so the two must agree. They do by construction: + _regenerate_projectors is the only producer of either, and it runs + dmftproj and the converter together. """ if not mpi.is_master_node(): return - _, windows = self._read_oubwin() + windows = self._read_oubwin() if len(windows) != N_k.shape[0]: raise DFTWorkflowError( f"case.oubwin has {len(windows)} k-points but N_k has {N_k.shape[0]}") @@ -859,10 +882,7 @@ def run_update_stage(self, N_k, Eint_m_dc, mu=0.0, beta=0.0, **kwargs): if mpi.is_master_node() and os.path.isfile('fort.77'): os.remove('fort.77') - # Four decimals, not the full float repr: a single charge update carries - # the solver's statistical error through N_k, so consecutive updates - # scatter by a few meV even after the density has settled. Printing 17 - # significant digits advertises a precision the number does not have. - self._report(f"DFT + DMFT Total Energy: {self.read_dft_energy():.4f} eV") + self._report(f"DFT + DMFT total energy after the charge update: " + f"{self.read_dft_energy():.4f} eV") mpi.barrier(poll_msec=100) return 0 From 707abf276d9b751072c25949aaeeb8e1b8ca441d Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 10 Sep 2026 14:05:24 -0400 Subject: [PATCH 11/14] [wien2k] Trim driver comments --- python/triqs_dftkit/wien2k/driver.py | 35 ++++------------------------ 1 file changed, 4 insertions(+), 31 deletions(-) diff --git a/python/triqs_dftkit/wien2k/driver.py b/python/triqs_dftkit/wien2k/driver.py index 14cbed0..787f926 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -15,8 +15,7 @@ -- ``run_lapw`` calls out to a user python script from inside its cycle, whereas DftDriver requires python to be the caller. -Scope: serial, non-magnetic (SP=0, SO=0). ``lapw2 -qdmft`` cannot be run in -parallel at all. +Scope: serial, non-magnetic (SP=0, SO=0). TODO: SOC and SP support will be added in the future """ import os, re, shutil, subprocess from datetime import datetime @@ -35,14 +34,13 @@ class DFTWorkflowError(Exception): # Default environment variables to preserve for subprocess execution. WIENROOT # and SCRATCH are WIEN2k specific: x interpolates $WIENROOT into every .def file -# it writes (e.g. the xc_funcs.h path for lapw0) and tcsh aborts outright on an -# undefined variable, so the environment cannot simply be stripped to the -# defaults the other drivers use. +# and tcsh aborts outright on an undefined variable, so the environment cannot simply +# be stripped to the defaults the other drivers use. _DEFAULT_ENV_VARS = ['PATH', 'LD_LIBRARY_PATH', 'SHELL', 'PWD', 'HOME', 'OMP_NUM_THREADS', 'OMP_STACKSIZE', 'MKL_NUM_THREADS', 'LANG', 'LC_ALL', 'OMPI_MCA_btl_vader_single_copy_mechanism', 'WIENROOT', 'SCRATCH'] -# WIEN2k works in Rydberg, TRIQS in eV. +# WIEN2k works in Rydberg, TRIQS in eV. _RY_IN_EV = 13.605698 # Flags that make x rewrite the first five characters of case.in2 in place and @@ -56,7 +54,6 @@ class DFTWorkflowError(Exception): # Files saved to _old before lapw0 and before mixer as run_lapw. # case.clmsum_old is mixer's previous-iteration density -# on unit 10, so the second set is required for mixing to work at all. _LAPW0_SAVE = ('vsp', 'vns', 'r2v') _MIXER_SAVE = ('clmsum', 'vrespsum', 'tausum') @@ -201,12 +198,6 @@ def _env(self): def _run_x(self, program, *flags): """ Run ``x `` and verify it succeeded. Master node only. - - Success is checked twice, because WIEN2k signals failure two ways: x - exits 9 when the program returns non-zero, and every program writes a - message into ``.error`` on entry and truncates it again just - before a successful exit. The error file therefore reports failures - that never reach the exit status, including untrapped runtime errors. """ if not mpi.is_master_node(): return 0 @@ -230,12 +221,6 @@ def _run_x(self, program, *flags): def _run_checked(self, command, error_file, label): """ Run a WIEN2k command and verify it succeeded. - - Success is checked twice, because WIEN2k signals failure two ways: a - non-zero exit status, and a message left in the .error file. Every program - writes that message on entry and truncates it again just before a - successful exit, so it catches failures that never reach the exit status, - including untrapped runtime errors. """ self._report(f"[{datetime.now()}] running {' '.join(command)}", level=2) result = subprocess.run(command, capture_output=True, text=True, env=self._env()) @@ -363,14 +348,6 @@ def _scf_tags(path, stop_at_dmft=False): maximum, which is what testconv compares against the -cc limit, and the trailing one is the cell total. - Pairing is done by record rather than by zipping two independent lists. - mixer writes :DIS only inside, i.e. only when it could read case.clmsum_old, - but writes :ENE unconditionally, so a cycle with no previous density contributes - an :ENE with no :DIS. Zipping and padding the shorter list at the end would attribute - every later :DIS to the cycle before its own and leave the newest cycle - with ``dis=None``, which silently disables the charge convergence test. - Within one record :DIS precedes :ENE, so :ENE closes the record. - ``stop_at_dmft`` stops at the first _DMFT_MARKER, leaving the caller only the plain-DFT part of the history. """ @@ -688,10 +665,6 @@ def run_dft_only(self, n_iter=None, force_scf=False): def _scf_history_on_disk(self): """ (ene, dis) pairs recovered from an existing case.scf, for restart. - - Stops at the first DMFT marker: the cycles after it were produced by - run_update_stage against a DMFT-corrected density, and the DFT ecut/ccut - limits say nothing useful about them. """ return self._scf_tags(self._f('scf'), stop_at_dmft=True) From ae330d2a58f4f0e3303af8ba78f4d67d06aef14d Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 10 Sep 2026 14:06:03 -0400 Subject: [PATCH 12/14] [wien2k] Guard n_loops in the CSC example --- doc/examples/wien2k_csc_svo/wien2k_modest_csc.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/doc/examples/wien2k_csc_svo/wien2k_modest_csc.py b/doc/examples/wien2k_csc_svo/wien2k_modest_csc.py index 726848c..33459f4 100644 --- a/doc/examples/wien2k_csc_svo/wien2k_modest_csc.py +++ b/doc/examples/wien2k_csc_svo/wien2k_modest_csc.py @@ -60,6 +60,10 @@ def dmft_cycle(obe, target_density, Sigma_imp_dyn, Sigma_imp_hf, Sigma_dc, n_loo impurity interaction energy minus the DC energy in eV, and the band-basis charge density correction N_k[k, sigma, nu, nu']. """ + if n_loops < 1: + raise ValueError("n_loops must be at least 1: the double counting, Eint-Edc " + "and N_k below are all built from the impurity solution") + epsilon_d = E.extract(tm.impurity_levels(obe))[0] for n in range(n_loops): From 1f7ea18ae993cfdf0bbb740aebb23faf7b82e26e Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Thu, 10 Sep 2026 14:06:04 -0400 Subject: [PATCH 13/14] [wien2k] Note the master-rank DFT chain in the example README --- doc/examples/wien2k_csc_svo/README.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/doc/examples/wien2k_csc_svo/README.md b/doc/examples/wien2k_csc_svo/README.md index 010cc08..27e1539 100644 --- a/doc/examples/wien2k_csc_svo/README.md +++ b/doc/examples/wien2k_csc_svo/README.md @@ -73,6 +73,10 @@ window. `lapw2 -qdmft` runs serially regardless of how the rest of the calculation is parallelised. The impurity solver still uses all available ranks. +The entire Wien2k chain runs on the master rank alone. Typically, the job size +is determined by the impurity solver, as we expect the DFT part to take a small +wall time. + The driver takes a `verbosity` argument: `1` (the default) prints one line per SCF cycle, `2` adds one line per Wien2k program launched, and `0` leaves only warnings. Warnings ignore the setting and go to stderr, so quietening the progress output From ac643ce712e52bc7198ad6d04891dfd9ea54c385 Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Tue, 15 Sep 2026 05:48:13 -0400 Subject: [PATCH 14/14] [wien2k] Restore the load-bearing driver comments --- python/triqs_dftkit/wien2k/driver.py | 34 +++++++++++++++++++++------- 1 file changed, 26 insertions(+), 8 deletions(-) diff --git a/python/triqs_dftkit/wien2k/driver.py b/python/triqs_dftkit/wien2k/driver.py index 787f926..4fd1d9a 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -15,7 +15,8 @@ -- ``run_lapw`` calls out to a user python script from inside its cycle, whereas DftDriver requires python to be the caller. -Scope: serial, non-magnetic (SP=0, SO=0). TODO: SOC and SP support will be added in the future +Scope: serial, non-magnetic (SP=0, SO=0). ``lapw2 -qdmft`` cannot be run in +parallel at all. TODO: SOC and SP support will be added in the future """ import os, re, shutil, subprocess from datetime import datetime @@ -34,13 +35,13 @@ class DFTWorkflowError(Exception): # Default environment variables to preserve for subprocess execution. WIENROOT # and SCRATCH are WIEN2k specific: x interpolates $WIENROOT into every .def file -# and tcsh aborts outright on an undefined variable, so the environment cannot simply +# and tcsh aborts outright on an undefined variable, so the environment cannot simply # be stripped to the defaults the other drivers use. _DEFAULT_ENV_VARS = ['PATH', 'LD_LIBRARY_PATH', 'SHELL', 'PWD', 'HOME', 'OMP_NUM_THREADS', 'OMP_STACKSIZE', 'MKL_NUM_THREADS', 'LANG', 'LC_ALL', 'OMPI_MCA_btl_vader_single_copy_mechanism', 'WIENROOT', 'SCRATCH'] -# WIEN2k works in Rydberg, TRIQS in eV. +# WIEN2k works in Rydberg, TRIQS in eV. _RY_IN_EV = 13.605698 # Flags that make x rewrite the first five characters of case.in2 in place and @@ -54,6 +55,7 @@ class DFTWorkflowError(Exception): # Files saved to _old before lapw0 and before mixer as run_lapw. # case.clmsum_old is mixer's previous-iteration density +# on unit 10, so the second set is required for mixing to work at all. _LAPW0_SAVE = ('vsp', 'vns', 'r2v') _MIXER_SAVE = ('clmsum', 'vrespsum', 'tausum') @@ -66,7 +68,7 @@ class DFTWorkflowError(Exception): # tests against the -cc limit; the trailing number is the cell total. _DIS_RE = re.compile(r'\(\s*([-+0-9.EDed]+)\s+for atom') -# DMFT label marker in the : LABEL style of run_lapw +# DMFT label marker in the : LABEL style of run_lapw _DMFT_MARKER = ':QDMFT: CHARGE DENSITY UPDATED FROM DMFT OCCUPATIONS' @@ -198,6 +200,12 @@ def _env(self): def _run_x(self, program, *flags): """ Run ``x `` and verify it succeeded. Master node only. + + Success is checked twice, because WIEN2k signals failure two ways: x + exits 9 when the program returns non-zero, and every program writes a + message into ``.error`` on entry and truncates it again just + before a successful exit. The error file therefore reports failures + that never reach the exit status, including untrapped runtime errors. """ if not mpi.is_master_node(): return 0 @@ -221,6 +229,12 @@ def _run_x(self, program, *flags): def _run_checked(self, command, error_file, label): """ Run a WIEN2k command and verify it succeeded. + + Success is checked twice, because WIEN2k signals failure two ways: a + non-zero exit status, and a message left in the .error file. Every program + writes that message on entry and truncates it again just before a + successful exit, so it catches failures that never reach the exit status, + including untrapped runtime errors. """ self._report(f"[{datetime.now()}] running {' '.join(command)}", level=2) result = subprocess.run(command, capture_output=True, text=True, env=self._env()) @@ -326,8 +340,8 @@ def _save_old(self, exts): run_lapw does this at two points in every cycle and both are load bearing. Most importantly mixer reads case.clmsum_old on unit 10 as the - previous-iteration density: without the copy it mixes against a missing - or stale density, which emits no :DIS line and sends the SCF diverging until + previous-iteration density: without the copy it mixes against a missing + or stale density, which emits no :DIS line and sends the SCF diverging until the linearisation energies go bad. """ if not mpi.is_master_node(): @@ -348,6 +362,10 @@ def _scf_tags(path, stop_at_dmft=False): maximum, which is what testconv compares against the -cc limit, and the trailing one is the cell total. + Pair by record, not by zipping: mixer writes :ENE every cycle but :DIS + only when it had a previous density. :DIS precedes :ENE, which closes + the record. + ``stop_at_dmft`` stops at the first _DMFT_MARKER, leaving the caller only the plain-DFT part of the history. """ @@ -486,7 +504,7 @@ def _converged(self, history): Apply testconv's criterion to the (ene, dis) history. The energy test is the mean of the last two |dE| over three iterations, - not a single difference, and needs three points. + not a single difference, and needs three points. A newest record with no :DIS counts as not converged. mixer omits :DIS only when it had no case.clmsum_old to compare against, which makes that @@ -716,7 +734,7 @@ def _write_qdmft(self, N_k, Eint_m_dc, mu=0.0, beta=0.0): nn must equal nb_top - nb_bot + 1 nn records of 2*nn reals: Re Im, one record per matrix *row* one throwaway record - correner in eV; lapw2 divides it by _RY_IN_EV + correner in eV; lapw2 divides it by 13.605698 Three conventions differ from the VASP and QE writers: