From 2dd14bf8a3f9f0679117fa15643b4e9702ad5f54 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 16 Jun 2026 14:42:23 +0200 Subject: [PATCH 1/3] [wien2k] pure-Python case.oubwin band-window generator Port the dmftproj almblm reader and case.oubwin writer to NumPy. Given case.almblm{up,dn} and the energy window in case.indmftpr, select the contiguous band range per k-point and spin and write case.oubwin{up,dn} in the dmftproj format. Support proj_mode 0; band-index modes remain explicitly unsupported. Part of #7. Keep window-open and window-valid flags separate, as set_projections.f does. A window opening only at nbmax remains excluded with zero bounds, including nbmin == nbmax. Preserve this native behavior on both sides. Reuse the existing compressed CaOs2 almblm fixtures, decompressing into a temporary directory. Compare both spin outputs with CaOs2.oubwinup, ignoring only trailing whitespace stripped from the shipped reference. Add focused checks for the final-band edge cases, energy boundaries, Fermi shifts, per-k lowest-band weights and excluded output records. Validation: the generator test passes locally with NumPy, including both reviewer reproducers and regeneration of the 43-92 CaOs2 band window. --- python/triqs_dftkit/wien2k/oubwin.py | 270 +++++++++++++++++++++ test/python/wien2k/CMakeLists.txt | 9 + test/python/wien2k/wien2k_oubwin_python.py | 97 ++++++++ 3 files changed, 376 insertions(+) create mode 100644 python/triqs_dftkit/wien2k/oubwin.py create mode 100644 test/python/wien2k/wien2k_oubwin_python.py diff --git a/python/triqs_dftkit/wien2k/oubwin.py b/python/triqs_dftkit/wien2k/oubwin.py new file mode 100644 index 0000000..b5452c8 --- /dev/null +++ b/python/triqs_dftkit/wien2k/oubwin.py @@ -0,0 +1,270 @@ +################################################################################ +# +# TRIQS: a Toolbox for Research in Interacting Quantum Systems +# +# Copyright (C) 2011 by L. Pourovskii, V. Vildosola, C. Martins, M. Aichhorn +# +# TRIQS is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# TRIQS. If not, see . +# +################################################################################ + +"""Pure-Python generation of the dmftproj band-window file (case.oubwin), +replacing the corresponding output of the dmftproj Fortran executable. + +Given the alm/blm coefficient file (case.almblm, or the spin pair +case.almblmup / case.almblmdn) and the projector definition (case.indmftpr), +this selects, per k-point and spin, the contiguous band range whose Fermi- +shifted Kohn-Sham eigenvalues fall in the energy window, and writes the +result in the case.oubwin format read by the charge self-consistency. + +It reproduces the Fortran path: read elecn/nk/nloat/eferm from the almblm +header (outbwin.f / dmftproj.f), parse the energy window (e_bot, e_top, +proj_mode) from the last line of case.indmftpr, shift every band eigenvalue +by eferm, and apply the proj_mode==0 selection of set_projections.f +(strict lower bound e_bot < E, inclusive upper bound E <= e_top, yielding a +single contiguous [nb_bot, nb_top] per k). The weight written per included +k-point is the tetrahedron weight of the lowest band nbmin. + +Spin-polarization is detected from the input files: case.almblmup / +case.almblmdn present => two spin files (case.oubwinup / case.oubwindn), +otherwise the single case.oubwin. The spin-orbit flag (ifSO) comes from +case.indmftpr; with -so the up and dn windows must coincide and dmftproj +aborts otherwise, a check this generator also performs. +""" + +import os +import numpy as np + + +# --- list-directed token stream over a Fortran free-form text file ----------- + +class _Reader: + """Record-oriented reader mirroring Fortran list-directed and READ('()') + semantics on a free-form almblm file. + + `record()` returns the whitespace tokens of the next physical line (the + common case: every almblm logical record fits on one line). `skip()` + consumes one physical line unconditionally, matching READ('()').""" + + def __init__(self, path): + with open(path) as fh: + self._lines = fh.read().splitlines() + self._i = 0 + + def skip(self): + self._i += 1 + + def record(self): + toks = self._lines[self._i].split() + self._i += 1 + return toks + + +def _to_float(tok): + """Parse a Fortran real token, accepting D/d exponents (1.0D-3) and the + list-directed 4.57E-002 form. Python float already handles E/e and the + embedded sign in the exponent; only D/d needs translation.""" + return float(tok.replace('D', 'E').replace('d', 'e')) + + +# --- inputs ------------------------------------------------------------------ + +def _read_window(indmftpr): + """Energy/band window from the last non-empty line of case.indmftpr. + + Try list-directed `e_bot e_top proj_mode`; on failure fall back to + `e_bot e_top` with proj_mode=0 (dmftproj.f:182-197).""" + lines = [l for l in open(indmftpr).read().splitlines() if l.strip()] + toks = lines[-1].split() + e_bot, e_top = _to_float(toks[0]), _to_float(toks[1]) + proj_mode = int(toks[2]) if len(toks) >= 3 else 0 + return e_bot, e_top, proj_mode + + +def _read_indmftpr_geometry(indmftpr): + """Return (nsort, lmax) from case.indmftpr. nsort is its first value and + lmax its third (dmftproj.f:95). Both are needed to walk the per-sort + almblm header up to the k-loop.""" + raw = [l.split('!')[0].strip() for l in open(indmftpr)] + raw = [l for l in raw if l != ''] + nsort = int(raw[0].split()[0]) + lmax = int(raw[2].split()[0]) + return nsort, lmax + + +def _read_so_flag(indmftpr): + """Spin-orbit flag from case.indmftpr: the single integer line that + follows a correlated sort's irep specification (the value dmftproj stores + as ifSOflag). Returns 0 or 1.""" + raw = [l.split('!')[0].strip() for l in open(indmftpr)] + raw = [l for l in raw if l != ''] + nsort = int(raw[0].split()[0]) + i = 3 + so = 0 + for _ in range(nsort): + basis = raw[i].split()[0] + i += 1 + if basis == 'fromfile': + i += 1 + l_inc = [int(x) for x in raw[i].split()] + i += 1 + ireps = [int(x) for x in raw[i].split()] + i += 1 + correlated = any(v == 2 for v in l_inc) + if any(n > 0 for n in ireps): + i += 1 + if correlated: + so = int(raw[i].split()[0]) + i += 1 + return so + + +# --- almblm parsing ---------------------------------------------------------- + +def _read_almblm_windows(path, nsort, lmax): + """Read case.almblm and return (eferm, [(nbmin, nbmax, eband, weight) per + k-point]) for the FIRST atomic sort. nbmin/nbmax/eband/weight are identical + across sorts at a given k, so the first sort's k-loop suffices for window + selection. + + Follows the dmftproj.f read order field by field: header, then the first + sort's per-l header (u_dot_norm, nLO, nLO overlap lines), then its k-loop + (two skipped banner lines, the idum/nbmin/nbmax line, the per-band + rtetr/eband lines).""" + r = _Reader(path) + r.record() # elecn + nk = int(r.record()[0]) + r.record() # nloat + eferm = _to_float(r.record()[0]) # non-band path: eferm is header line 4 + + # First sort's per-l header: u_dot_norm, nLO, then nLO overlap lines. + for _ in range(lmax + 1): + r.record() # u_dot_norm(l) + nlo = int(r.record()[0]) + for _ in range(nlo): + r.record() # ovl_LO_u, ovl_LO_udot + + kpoints = [] + for _ in range(nk): + r.skip() # "IK = .. jatom=.." banner + r.skip() # 3-int line + idum_nb = r.record() + nbmin, nbmax = int(idum_nb[1]), int(idum_nb[2]) + eband = np.empty(nbmax - nbmin + 1) + weight = None + for off in range(nbmax - nbmin + 1): + toks = r.record() + rtetr = _to_float(toks[0]) + eband[off] = _to_float(toks[1]) + if off == 0: # tetrahedron weight of the lowest band + weight = rtetr + kpoints.append((nbmin, nbmax, eband, weight)) + return eferm, kpoints + + +# --- band-window selection (set_projections.f, proj_mode==0) ----------------- + +def _select_window(nbmin, nbmax, eband, eferm, e_bot, e_top): + """proj_mode==0 contiguous band selection for one k-point. + + Scan absolute band index ib = nbmin..nbmax over the Fermi-shifted energy + E = eband - eferm. The first ib with e_bot < E <= e_top opens the window; + the first ib above it with E > e_top closes it at ib-1 (scan exits). A + window already open on reaching nbmax closes at nbmax. As in Fortran, + opening at nbmax alone does not validate the window. Returns + (included, nb_bot, nb_top), with zero bounds for an invalid window.""" + opened = False + included = False + nb_bot = nb_top = 0 + for off, ib in enumerate(range(nbmin, nbmax + 1)): + e = eband[off] - eferm + if not opened and e_bot < e <= e_top: + opened = True + nb_bot = ib + elif opened and e > e_top: + nb_top = ib - 1 + included = True + break + elif ib == nbmax and e_bot < e <= e_top: + nb_top = ib + included = True + if not included: + nb_bot = nb_top = 0 + return included, nb_bot, nb_top + + +def _windows_for_spin(almblm, indmftpr, nsort, lmax, e_bot, e_top, proj_mode): + """Per-k (included, nb_bot, nb_top, weight) for one spin's almblm file.""" + if proj_mode != 0: + raise NotImplementedError( + 'oubwin generation is implemented for the energy-window ' + 'projection mode (proj_mode==0); band-index modes 1 and 2 are ' + 'not yet covered by a test fixture') + eferm, kpoints = _read_almblm_windows(almblm, nsort, lmax) + out = [] + for nbmin, nbmax, eband, weight in kpoints: + incl, nb_bot, nb_top = _select_window( + nbmin, nbmax, eband, eferm, e_bot, e_top) + out.append((incl, nb_bot, nb_top, weight)) + return out + + +# --- output ------------------------------------------------------------------ + +def _write_oubwin_file(path, ifso, windows): + """Write one case.oubwin spin file (outbwin.f). Integer fields are i6; + the per-k weight is a list-directed real. Not-included k-points emit only + the flag line.""" + with open(path, 'w') as f: + f.write('%6d\n' % len(windows)) + f.write('%6d\n' % (1 if ifso else 0)) + for incl, nb_bot, nb_top, weight in windows: + f.write('%6d\n' % (1 if incl else 0)) + if incl: + f.write('%6d%6d\n' % (nb_bot, nb_top)) + f.write(' %.16f \n' % weight) + + +def write_oubwin(case): + """Read .almblm{up,dn} (or .almblm) and write the matching + .oubwin{up,dn} (or .oubwin) in the dmftproj format. + + Spin-polarization is inferred from the input files. The spin-orbit flag + written into every file comes from case.indmftpr; under -sp+-so the up/dn + windows must coincide (dmftproj aborts otherwise).""" + indmftpr = case + '.indmftpr' + nsort, lmax = _read_indmftpr_geometry(indmftpr) + ifso = _read_so_flag(indmftpr) + e_bot, e_top, proj_mode = _read_window(indmftpr) + + up, dn = case + '.almblmup', case + '.almblmdn' + spin_polarized = os.path.exists(up) and os.path.exists(dn) + + if spin_polarized: + win_up = _windows_for_spin( + up, indmftpr, nsort, lmax, e_bot, e_top, proj_mode) + win_dn = _windows_for_spin( + dn, indmftpr, nsort, lmax, e_bot, e_top, proj_mode) + if ifso: + for u, d in zip(win_up, win_dn): + if u[0] != d[0] or u[1] != d[1] or u[2] != d[2]: + raise ValueError( + 'spin-orbit run requires identical up/dn band ' + 'windows at every k-point') + _write_oubwin_file(case + '.oubwinup', ifso, win_up) + _write_oubwin_file(case + '.oubwindn', ifso, win_dn) + else: + win = _windows_for_spin( + case + '.almblm', indmftpr, nsort, lmax, e_bot, e_top, proj_mode) + _write_oubwin_file(case + '.oubwin', ifso, win) diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt index c43cbed..05031b5 100644 --- a/test/python/wien2k/CMakeLists.txt +++ b/test/python/wien2k/CMakeLists.txt @@ -42,6 +42,15 @@ add_test(NAME Py_wien2k_symqmc_mixing set_property(TEST Py_wien2k_symqmc_mixing APPEND PROPERTY ENVIRONMENT PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD}) +# Pure-Python case.oubwin band-window generator vs the dmftproj Fortran output. +# Reads the existing compressed almblm inputs from the source fixtures and +# compares both spin outputs against CaOs2.oubwinup in a scratch directory. +add_test(NAME Py_wien2k_oubwin_python + COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_oubwin_python.py + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_property(TEST Py_wien2k_oubwin_python APPEND PROPERTY ENVIRONMENT + PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD}) + add_test(NAME Py_wien2k_kpath_labels COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_kpath_labels.py WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/test/python/wien2k/wien2k_oubwin_python.py b/test/python/wien2k/wien2k_oubwin_python.py new file mode 100644 index 0000000..1b0c17b --- /dev/null +++ b/test/python/wien2k/wien2k_oubwin_python.py @@ -0,0 +1,97 @@ +################################################################################ +# TRIQS: a Toolbox for Research in Interacting Quantum Systems +# (GPL-3.0-or-later) +################################################################################ + +# Verify the pure-Python case.oubwin generator (triqs_dftkit.wien2k.oubwin) +# reproduces the dmftproj Fortran output: regenerate case.oubwinup / case.oubwindn +# from case.almblm{up,dn} + case.indmftpr and compare to the native reference, +# ignoring only trailing whitespace removed from the shipped fixture. +# +# CaOs2 is spin-orbit + spin-polarized, so both outputs match CaOs2.oubwinup. +# Reuse the compressed almblm inputs shipped with the SOC fixture. + +import gzip +import os +import shutil +import tempfile +import numpy as np + +from triqs_dftkit.wien2k import oubwin + +HERE = os.path.dirname(os.path.abspath(__file__)) + + +def _check(case): + tmp = tempfile.mkdtemp() + try: + for ext in ('almblmup', 'almblmdn'): + with gzip.open(os.path.join(HERE, f'{case}.{ext}.gz'), 'rb') as src, \ + open(os.path.join(tmp, f'{case}.{ext}'), 'wb') as dst: + shutil.copyfileobj(src, dst) + shutil.copy(os.path.join(HERE, f'{case}.indmftpr'), tmp) + oubwin.write_oubwin(os.path.join(tmp, case)) + with open(os.path.join(HERE, f'{case}.oubwinup')) as f: + ref = [line.rstrip() for line in f] + for ext in ('oubwinup', 'oubwindn'): + with open(os.path.join(tmp, f'{case}.{ext}')) as f: + got = [line.rstrip() for line in f] + assert got == ref, f'{case}.{ext} differs from dmftproj reference' + finally: + shutil.rmtree(tmp) + + +# Independent expected band ranges for set_projections.f's branch ordering. +for nbmin, energies, expected in ( + (10, [-3., -2., -1., 0.5], (False, 0, 0)), # opens only at nbmax + (5, [0.5], (False, 0, 0)), # nbmin == nbmax + (10, [0., 0.5, 1., 2.], (True, 11, 12)), # strict lower, inclusive upper + (10, [0.5, 1.], (True, 10, 11)), # open window reaches nbmax + (10, [-1., 0., 2.], (False, 0, 0))): # no in-window band + got = oubwin._select_window(nbmin, nbmin + len(energies) - 1, + np.array(energies), 0., 0., 1.) + assert got == expected, (energies, got, expected) + +# Three k-points with distinct weights, a Fermi shift of 10 Ry, and one +# excluded point. The lowest-band weight is used even if that band is outside +# the window; an excluded point emits no band range or weight record. +with tempfile.TemporaryDirectory() as tmp: + almblm = os.path.join(tmp, 'case.almblm') + with open(almblm, 'w') as f: + f.write('''0 +3 +0 +10.0 +1.0 +0 +IK = 1 +0 0 0 +1 10 13 +0.25 9.0 +0.80 10.25 +0.60 10.75 +0.25 12.0 +IK = 2 +0 0 0 +2 20 21 +0.10 8.0 +0.10 9.0 +IK = 3 +0 0 0 +3 30 31 +0.75 10.25 +0.50 10.75 +''') + windows = oubwin._windows_for_spin(almblm, None, 1, 0, 0., 1., 0) + assert windows == [(True, 11, 12, 0.25), (False, 0, 0, 0.10), + (True, 30, 31, 0.75)], windows + output = os.path.join(tmp, 'case.oubwin') + oubwin._write_oubwin_file(output, 0, windows) + with open(output) as f: + assert [line.split() for line in f] == [ + ['3'], ['0'], ['1'], ['11', '12'], ['0.2500000000000000'], + ['0'], ['1'], ['30', '31'], ['0.7500000000000000']] + +_check('CaOs2') + +print('wien2k_oubwin_python: ok') From 20ae97b4c64c283364d9505ad38c61aba20f9151 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 12 Oct 2026 19:36:43 +0200 Subject: [PATCH 2/3] Read each correlated-shell irep record before the SO flag --- python/triqs_dftkit/wien2k/oubwin.py | 14 +++++++++++--- test/python/wien2k/wien2k_oubwin_python.py | 6 ++++++ 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/python/triqs_dftkit/wien2k/oubwin.py b/python/triqs_dftkit/wien2k/oubwin.py index b5452c8..0316b60 100644 --- a/python/triqs_dftkit/wien2k/oubwin.py +++ b/python/triqs_dftkit/wien2k/oubwin.py @@ -112,6 +112,7 @@ def _read_so_flag(indmftpr): nsort = int(raw[0].split()[0]) i = 3 so = 0 + seen_correlated = False for _ in range(nsort): basis = raw[i].split()[0] i += 1 @@ -122,10 +123,17 @@ def _read_so_flag(indmftpr): ireps = [int(x) for x in raw[i].split()] i += 1 correlated = any(v == 2 for v in l_inc) - if any(n > 0 for n in ireps): - i += 1 + for l, included in enumerate(l_inc): + if included == 2 and ireps[l] > 0: + i += 1 if correlated: - so = int(raw[i].split()[0]) + sort_so = int(raw[i].split()[0]) + if sort_so not in (0, 1): + raise ValueError('spin-orbit flags must be 0 or 1') + if seen_correlated and sort_so != so: + raise NotImplementedError('different spin-orbit flags across correlated sorts are not implemented') + so = sort_so + seen_correlated = True i += 1 return so diff --git a/test/python/wien2k/wien2k_oubwin_python.py b/test/python/wien2k/wien2k_oubwin_python.py index 1b0c17b..93a364f 100644 --- a/test/python/wien2k/wien2k_oubwin_python.py +++ b/test/python/wien2k/wien2k_oubwin_python.py @@ -94,4 +94,10 @@ def _check(case): _check('CaOs2') +with tempfile.TemporaryDirectory() as tmp: + path = os.path.join(tmp, 'multi.indmftpr') + with open(path, 'w') as fh: + fh.write('1\n1\n2\ncomplex\n0 2 2\n0 1 1\n1\n1\n0\n-2 3\n') + assert oubwin._read_so_flag(path) == 0 + print('wien2k_oubwin_python: ok') From 87be8ca73c4366191455f5d546ab9b2e31aedc9c Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 12 Oct 2026 19:39:25 +0200 Subject: [PATCH 3/3] Accept inline comments after the projection window --- python/triqs_dftkit/wien2k/oubwin.py | 3 ++- test/python/wien2k/wien2k_oubwin_python.py | 3 +++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/python/triqs_dftkit/wien2k/oubwin.py b/python/triqs_dftkit/wien2k/oubwin.py index 0316b60..a575536 100644 --- a/python/triqs_dftkit/wien2k/oubwin.py +++ b/python/triqs_dftkit/wien2k/oubwin.py @@ -85,7 +85,8 @@ def _read_window(indmftpr): Try list-directed `e_bot e_top proj_mode`; on failure fall back to `e_bot e_top` with proj_mode=0 (dmftproj.f:182-197).""" - lines = [l for l in open(indmftpr).read().splitlines() if l.strip()] + lines = [l.split('!')[0].strip() for l in open(indmftpr)] + lines = [l for l in lines if l] toks = lines[-1].split() e_bot, e_top = _to_float(toks[0]), _to_float(toks[1]) proj_mode = int(toks[2]) if len(toks) >= 3 else 0 diff --git a/test/python/wien2k/wien2k_oubwin_python.py b/test/python/wien2k/wien2k_oubwin_python.py index 93a364f..26b427a 100644 --- a/test/python/wien2k/wien2k_oubwin_python.py +++ b/test/python/wien2k/wien2k_oubwin_python.py @@ -99,5 +99,8 @@ def _check(case): with open(path, 'w') as fh: fh.write('1\n1\n2\ncomplex\n0 2 2\n0 1 1\n1\n1\n0\n-2 3\n') assert oubwin._read_so_flag(path) == 0 + with open(path, 'a') as fh: + fh.write('-2 3 ! default energy-window mode\n') + assert oubwin._read_window(path) == (-2., 3., 0) print('wien2k_oubwin_python: ok')