diff --git a/doc/examples/wien2k_csc_svo/README.md b/doc/examples/wien2k_csc_svo/README.md new file mode 100644 index 0000000..27e1539 --- /dev/null +++ b/doc/examples/wien2k_csc_svo/README.md @@ -0,0 +1,85 @@ +# 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 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 +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..4d0ed1d --- /dev/null +++ b/doc/examples/wien2k_csc_svo/SrVO3.indmftpr @@ -0,0 +1,15 @@ +3 +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 +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..33459f4 --- /dev/null +++ b/doc/examples/wien2k_csc_svo/wien2k_modest_csc.py @@ -0,0 +1,170 @@ +# ============================================================================ +# 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']. + """ + 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): + 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() 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..4fd1d9a 100644 --- a/python/triqs_dftkit/wien2k/driver.py +++ b/python/triqs_dftkit/wien2k/driver.py @@ -1,48 +1,879 @@ +""" +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. + +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. TODO: SOC and SP support will be added in the future +""" +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 +# 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. +_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') + +# 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. +_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') + +# DMFT label marker in the : LABEL style of run_lapw +_DMFT_MARKER = ':QDMFT: CHARGE DENSITY UPDATED FROM DMFT OCCUPATIONS' + 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'} + # 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}, " + 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 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 + # 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')}" + "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 _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 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(): + self._remove_broyden_files() + 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}; 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): + """ + 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: 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 + 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, stop_at_dmft=False): + """ + 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. + + 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. + """ + history, dis = [], None + if not os.path.isfile(path): + 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'): + history.append((float(line.split()[-1]), dis)) + dis = None + elif line.startswith(':DIS'): + match = _DIS_RE.search(line) + dis = (float(match.group(1).replace('D', 'E').replace('d', 'e')) + if match else float(line.split()[-1])) + return history + + 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") + history = self._scf_tags(path) + if not history: + raise DFTWorkflowError(f"no :ENE line in {path}") + 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 ene, dis + + 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. + """ + 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") + # 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 = history[-1][0] * _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 + 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. + 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. + """ + 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_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. + + 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 not None and 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') + 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 + # 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): + """ + 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 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): + """ + 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)") + 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") + # 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: + 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") + 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. + + 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) + + 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, history=history) + + def _scf_history_on_disk(self): + """ + (ene, dis) pairs recovered from an existing case.scf, for restart. + """ + return self._scf_tags(self._f('scf'), stop_at_dmft=True) + + # --------------------------------------------------------- dmft interface + + def _read_oubwin(self): + """ + Read case.oubwin, written by dmftproj. + + 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): + 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)) + next(reader) # the SO flag + + 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 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 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* + 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. + + ``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() + 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: + # 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) + 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`` (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) + + # 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_mixer() + self._append_to_scf(('m',)) + + self._save_old(_LAPW0_SAVE) + self._run_x('lapw0') + self._run_x('lapw1') + self._regenerate_projectors() - 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 + if mpi.is_master_node() and os.path.isfile('fort.77'): + os.remove('fort.77') + 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