diff --git a/doc/ChangeLog.md b/doc/ChangeLog.md index f70657d..bc151f3 100644 --- a/doc/ChangeLog.md +++ b/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 diff --git a/doc/h5structure.rst b/doc/h5structure.rst index 0a9349c..b12206b 100644 --- a/doc/h5structure.rst +++ b/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:: diff --git a/python/triqs_dftkit/wien2k/converter.py b/python/triqs_dftkit/wien2k/converter.py index 2a6ac9a..37e25f3 100644 --- a/python/triqs_dftkit/wien2k/converter.py +++ b/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] diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt index d9968e3..05cd5ad 100644 --- a/test/python/wien2k/CMakeLists.txt +++ b/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}) diff --git a/test/python/wien2k/test_kpath_labels.py b/test/python/wien2k/test_kpath_labels.py new file mode 100644 index 0000000..b2f1f45 --- /dev/null +++ b/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 . +# +################################################################################ + +""" +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)