[wien2k] Store high-symmetry k-path labels in dft_bands_input

doc/ChangeLog.md

@@ -27,10 +27,13 @@ Find below an itemized list of changes in this release.
 * Add a VASP driver for charge self-consistent DFT+DMFT calculations
 * Add `KPOINTS_OPT` band conversion from `vaspout.h5`: when `LOCPROJ_OPT` data are available, the converter writes `dft_bands_input` for band/spectral workflows, applying the same PLO config settings (`EWINDOW`, `TRANSFORM`, `NORMALIZE`, and optional `EFERMI`) as the regular VASP conversion, and stores the high-symmetry k-path labels
 * Warn on misplaced or unknown tags in the PLOVASP configuration
+* Read VASP `ICHARG=5` miscellaneous input from `vaspout.h5`
 
 ### Wannier90
 * Add ABINIT support to `Wannier90Converter` for charge self-consistent calculations
-* Read VASP `ICHARG=5` miscellaneous input from `vaspout.h5`
+
+### Wien2k
+* Read the high-symmetry k-path labels from the end of `case.outband` and store them as `kpts_labels` / `kpts_labels_idx` in `dft_bands_input`, matching the VASP band conversion
 
 ### Fix
 * Fix a bug in the `deltaN` write for the Quantum Espresso and Abinit interfaces

doc/h5structure.rst

@@ -603,6 +603,17 @@ counting the k-points along the path.
        ``max(n_parproj)``, ``max(shells['dim'])``, ``max(n_orbitals)``]
      - As in ``dft_parproj_input``, along the path. Elk writes a dummy
        ``array([0])``.
+   * - ``kpts_labels``
+     - list of string
+     - Names of the high-symmetry points of the path (e.g. ``'GAMMA'``, ``'X'``),
+       for labelling the ticks of a band plot. Only written by Wien2k and VASP,
+       and only when the underlying DFT output provides the labels.
+   * - ``kpts_labels_idx``
+     - numpy.array.int, dim [``len(kpts_labels)``]
+     - Position of each entry of ``kpts_labels`` along the path, as a 0-based
+       index into the ``n_k`` k-points. Each converter validates these indices
+       as far as its own DFT output allows, so consumers should not rely on
+       them being sorted or in range without checking.
 
 .. note::
 

python/triqs_dftkit/wien2k/converter.py

@@ -30,6 +30,91 @@
 import os.path
 
 
+def _parse_kpath_label_line(line):
+    """
+    Parse one line of the high-symmetry point block that dmftproj appends to
+    case.outband. The format is the fixed Fortran (2i6,a): columns 1-6 hold a
+    running counter, columns 7-12 the 1-based position of the point along the
+    band path, and columns 13 onwards the label.
+
+    Returns (counter, pos, label), or None if the line does not have that form.
+    """
+    line = line.rstrip('\n')
+    if len(line) < 13:
+        return None
+    try:
+        counter = int(line[0:6])
+        pos = int(line[6:12])
+    except ValueError:
+        return None
+    label = line[12:].strip()
+    if not label:
+        return None
+    return counter, pos, label
+
+
+def _read_kpath_labels(band_file, n_k):
+    """
+    Read the high-symmetry k-path labels appended to the end of case.outband
+    and map them onto the flattened band k-point index.
+
+    dmftproj writes one line per high-symmetry point after the projector data,
+    in fixed Fortran format (2i6,a), e.g.
+
+            1     1GAMMA
+            2   122X
+
+    Returns (labels, idx) where labels is a list of label strings and idx is a
+    0-based numpy int array giving, for each label, the position of that
+    high-symmetry point in the n_k band path. Returns (None, None) if no label
+    block is present or if the block is incomplete.
+    """
+    # The block is a handful of lines at the very end of a file that holds all
+    # the projectors, so read the tail rather than the whole file.
+    n_bytes = 8192
+    with open(band_file, 'rb') as R:
+        R.seek(0, os.SEEK_END)
+        size = R.tell()
+        R.seek(max(0, size - n_bytes))
+        lines = R.read().decode('utf-8', 'replace').splitlines()
+    if size > n_bytes:
+        # the first line of the chunk is in general cut in the middle
+        lines = lines[1:]
+    while lines and not lines[-1].strip():
+        lines.pop()
+
+    # Walk backwards from the end of the file: the label block is the trailing
+    # run of (2i6,a) lines, and the counter of its first line is 1.
+    labels = []
+    idx = []
+    first_counter = None
+    for line in reversed(lines):
+        parsed = _parse_kpath_label_line(line)
+        if parsed is None:
+            break
+        first_counter, pos, label = parsed
+        labels.append(label)
+        idx.append(pos - 1)
+        if first_counter == 1:
+            break
+    labels.reverse()
+    idx.reverse()
+
+    if not labels:
+        return None, None
+
+    if first_counter != 1:
+        mpi.report("convert_bands_input : WARNING : the high-symmetry point block in %s does not start at counter 1, so it is truncated or corrupted; skipping k-path labels." % band_file)
+        return None, None
+
+    idx = numpy.array(idx, dtype=int)
+    if idx[0] < 0 or idx[-1] >= n_k or numpy.any(numpy.diff(idx) < 0):
+        mpi.report("convert_bands_input : WARNING : inconsistent high-symmetry point indices in %s; skipping k-path labels." % band_file)
+        return None, None
+
+    return labels, idx
+
+
 class Converter(ConverterTools):
     """
     Conversion from Wien2k output to an hdf5 file that can be used as input for the SumkDFT class.
@@ -483,6 +568,8 @@ def convert_bands_input(self):
 
         # Reading done!
 
+        kpts_labels, kpts_labels_idx = _read_kpath_labels(self.band_file, n_k)
+
         # Save it to the HDF:
         with HDFArchive(self.hdf_file, 'a') as ar:
             if not (self.bands_subgrp in ar):
@@ -491,6 +578,9 @@ def convert_bands_input(self):
             # created. If it exists, the data is overwritten!
             things_to_save = ['n_k', 'n_orbitals', 'proj_mat',
                           'hopping', 'n_parproj', 'proj_mat_all']
+            if kpts_labels is not None:
+                things_to_save += ['kpts_labels', 'kpts_labels_idx']
+                mpi.report("  Stored %i high-symmetry k-path labels: %s" % (len(kpts_labels), ', '.join(kpts_labels)))
             for it in things_to_save:
                 ar[self.bands_subgrp][it] = locals()[it]
 

test/python/wien2k/CMakeLists.txt

@@ -11,3 +11,9 @@ add_test(NAME Py_wien2k_convert
          WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
 set_property(TEST Py_wien2k_convert 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})
+set_property(TEST Py_wien2k_kpath_labels APPEND PROPERTY ENVIRONMENT
+             PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})

test/python/wien2k/test_kpath_labels.py

@@ -0,0 +1,74 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
+#
+# 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 <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+"""
+Test the high-symmetry k-path label block that dmftproj appends to case.outband.
+A real case.outband is far too large to ship as a fixture, so this writes a
+synthetic tail in the exact Fortran (2i6,a) format of outband.f:277, including
+the trailing blank the `a` descriptor emits for the declared length of kname.
+"""
+
+import os
+import tempfile
+
+import numpy as np
+
+from triqs_dftkit.wien2k.converter import _read_kpath_labels
+
+N_K = 501
+
+
+def line(counter, pos, label):
+    return '%6d%6d%s' % (counter, pos, label.ljust(10))
+
+
+def read(lines, n_k=N_K):
+    fd, path = tempfile.mkstemp(suffix='.outband')
+    try:
+        with os.fdopen(fd, 'w') as f:
+            f.write('\n'.join(lines) + '\n')
+        return _read_kpath_labels(path, n_k)
+    finally:
+        os.remove(path)
+
+
+# One block exercising the awkward cases at once: '\xG' does not start with a
+# letter (XCrySDen writes Gamma that way), X|Y puts two labels on the same
+# k-point, the projector line in front of the block itself parses as (2i6,a),
+# the padding pushes the block past the tail that is read back, and the file
+# ends on blank lines.
+padding = ['%20.14f%20.14f' % (0.5, 0.25)] * 2000
+block = [line(1, 1, '\\xG'), line(2, 122, 'X'), line(3, 122, 'Y'),
+         line(4, 268, 'L'), line(5, 501, 'K')]
+labels, idx = read(padding + [line(3, 7, '0.52')] + block + ['', '   '])
+
+assert labels == ['\\xG', 'X', 'Y', 'L', 'K'], labels
+np.testing.assert_array_equal(idx, [0, 121, 121, 267, 500])
+
+# A malformed line in the middle of the block truncates the backward walk, so
+# the block no longer starts at counter 1 and is rejected rather than silently
+# returning only its tail.
+broken = list(block)
+broken[2] = 'this line is not a label line'
+assert read(padding + broken) == (None, None)
+
+# A position past the end of the band path is rejected as well.
+assert read(padding + block, n_k=400) == (None, None)