From 2a570e6da6ba076144253b99eadc04f41aa7f1df Mon Sep 17 00:00:00 2001 From: Harrison LaBollita Date: Wed, 2 Sep 2026 08:45:14 -0400 Subject: [PATCH] [wien2k] Add init_dmftpr, the case.indmftpr generator dmftproj reads case.indmftpr, which init_lapw does not produce. Writing it by hand means matching a positional format with one lsort/lnreps entry per l, so this ports the interactive generator from triqs_dft_tools (bin/init_dmftpr): it reads case.struct for the species and multiplicities, asks what to project and how, and writes the file in the layout dmftproj expects. Installed with install(PROGRAMS) so the execute bit is set, unlike the install(FILES) dft_tools used. Fixes carried over the port: * the energy window is always written. Previously a third token that was not exactly 'ev'/'eV'/'Ev' -- e.g. '-0.6 0.14 Ry', or 'EV' -- fell through the unit test without writing anything, silently producing a case.indmftpr with no window line at all. * the irrep default was the string "0 0 0 0\n" passed to list.extend, which iterates it character by character and wrote eight separate lines. This hit every correlated d shell with no t2g/eg split, and every correlated f shell. * the struct file is taken as .struct rather than the first *.struct glob match, which could pick up a stray tmp.struct. * orbital letters are validated instead of raising KeyError. * the file is written only once every answer is in, so aborting part way through cannot leave a half-finished case.indmftpr on disk. * python3 shebang. --- CMakeLists.txt | 3 + bin/CMakeLists.txt | 3 + bin/init_dmftpr | 241 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 247 insertions(+) create mode 100644 bin/CMakeLists.txt create mode 100755 bin/init_dmftpr diff --git a/CMakeLists.txt b/CMakeLists.txt index 81d0f59..ba7e626 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -109,6 +109,9 @@ add_subdirectory(deps) # Build and install the dmftproj executable add_subdirectory(fortran/dmftproj) +# Command line scripts (init_dmftpr generates dmftproj's input file) +add_subdirectory(bin) + # Tests if(Build_Tests) add_subdirectory(test) diff --git a/bin/CMakeLists.txt b/bin/CMakeLists.txt new file mode 100644 index 0000000..94a0fa5 --- /dev/null +++ b/bin/CMakeLists.txt @@ -0,0 +1,3 @@ +# Command line scripts. install(PROGRAMS) sets the execute bit, unlike the +# install(FILES) that dft_tools used for init_dmftpr. +install(PROGRAMS init_dmftpr DESTINATION bin) diff --git a/bin/init_dmftpr b/bin/init_dmftpr new file mode 100755 index 0000000..66c64be --- /dev/null +++ b/bin/init_dmftpr @@ -0,0 +1,241 @@ +#!/usr/bin/env python3 +""" +Interactive generator for case.indmftpr, the input file of dmftproj. + +Ported from triqs_dft_tools (bin/init_dmftpr). Reads case.struct for the atom +species and multiplicities, asks what to project and how, and writes the file in +the layout dmftproj expects. WIEN2k ships a worked example at +$WIENROOT/SRC_templates/case.indmftpr. + +Run it in the WIEN2k case directory, after init_lapw: + + init_dmftpr + +Differences from the dft_tools original, all bug fixes: + * the energy window is always written. Previously a third token that was not + exactly 'ev'/'eV'/'Ev' -- e.g. '-0.6 0.14 Ry', or 'EV' -- fell through the + unit test without writing anything, silently producing a case.indmftpr with + no window line at all. + * the struct file is taken as .struct rather than the first *.struct + glob match, which could pick up a stray tmp.struct. + * orbital letters are validated instead of raising KeyError. + * python3 shebang, and installed with the execute bit set. +""" +import os +import sys + +from numpy import array + +# Occupancy vectors, indexed s, p, d, f. 2 marks a correlated shell, 1 an +# uncorrelated projector. +ORBITALS = {"s": [1, 0, 0, 0], + "p": [0, 1, 0, 0], + "d": [0, 0, 1, 0], + "f": [0, 0, 0, 1]} +CORR_ORBITALS = {"s": [2, 0, 0, 0], + "p": [0, 2, 0, 0], + "d": [0, 0, 2, 0], + "f": [0, 0, 0, 2]} + +EV_PER_RY = 13.605698 + +# dmftproj reads one lsort/lnreps entry per l from 0 to lmax. +LMAX = 3 + + +def parse_orbitals(text, table): + """Sum the occupancy vectors for a string like 'spd'; None if any letter is bad.""" + bad = [c for c in text if c not in table] + if bad: + print("Not an orbital: {}. Use s, p, d or f.".format(", ".join(sorted(set(bad))))) + return None + total = array([0, 0, 0, 0], dtype=int) + for c in text: + total += array(table[c], dtype=int) + return total + + +def parse_window(text): + """ + Turn a projection-window answer into the ' ' line, in Ry. + + Accepts 'emin emax' (Ry) or 'emin emax ' with unit ev/eV/Ry, case + insensitive. Returns None if the window does not straddle the Fermi energy + or cannot be read, so the caller can re-prompt. + """ + fields = text.split() + if len(fields) not in (2, 3): + print("Expected 'emin emax' or 'emin emax UNIT'.") + return None + try: + emin, emax = float(fields[0]), float(fields[1]) + except ValueError: + print("Could not read '{}' as two numbers.".format(text)) + return None + if not (emin < 0 < emax): + print("The energy window ({}) does not contain the Fermi energy!".format(text)) + return None + unit = fields[2].lower() if len(fields) == 3 else "ry" + if unit == "ry": + return "{0:0.5f} {1:0.5f}".format(emin, emax) + if unit == "ev": + return "{0:0.5f} {1:0.5f}".format(emin / EV_PER_RY, emax / EV_PER_RY) + print("Unknown unit '{}'; use Ry or eV.".format(fields[2])) + return None + + +def ask(prompt, allowed=None): + """Prompt until the answer is in `allowed` (or any non-empty answer if None).""" + while True: + answer = input(prompt) + if allowed is None or answer in allowed: + return answer + print("Did not recognize that input. Try again.") + + +def collect_indmftpr(struct_file): + """ + Ask the user what to project and return the case.indmftpr lines. + + Returns a list of lines (no trailing newlines). Nothing is written here, so + an abort part way through cannot leave a half-finished file on disk, and the + generated layout can be checked without driving stdin. + """ + lines = [] + struct = open(struct_file).readlines() + species = [line.split()[0] for line in struct if "NPT" in line] + mult = [line.split("=")[1].split()[0] for line in struct if "MULT" in line] + num_atoms = len(species) + if num_atoms == 0: + print("No atoms found in {} (no NPT lines).".format(struct_file)) + sys.exit(1) + + print("number of atoms = {} ({})\n".format(num_atoms, " ".join(species))) + + lines.append(str(num_atoms)) + lines.append(" ".join(mult)) + lines.append(str(LMAX)) + + for atom in range(num_atoms): + label = "ATOM {} ({})".format(atom + 1, species[atom]) + + sph_harm = ask("What flavor of spherical harmonics do you want to use " + "for {}? (cubic/complex/fromfile)\n".format(label), + ["cubic", "complex", "fromfile"]) + lines.append(sph_harm) + if sph_harm == "fromfile": + while True: + filename = input("name of file defining the basis?\n") + if not os.path.isfile(filename): + print("{} could not be found in the current directory!".format(filename)) + continue + if len(filename) >= 25: # dmftproj's field width + rename = input("{} is too long! Rename the file to: \n".format(filename)) + os.rename(filename, rename) + filename = rename + lines.append(filename) + break + + if ask("Do you want to treat {} as correlated (y/n)?\n".format(label), + ["y", "n"]) == "y": + proj = ask("Specify the correlated orbital? (d,f)\n", ["d", "f"]) + while True: + non_corr = input("projectors for non-correlated orbitals? " + "(type h for help, blank for none)\n") + if non_corr == "h": + print("indicate orbital projectors using (s, p, d, or f). " + "For multiple, combine them (sp, pd, spd, etc.)") + continue + if proj in non_corr: + print("Error: User can not choose orbital {} as both " + "correlated and uncorrelated!".format(proj)) + continue + projectors = array(CORR_ORBITALS[proj], dtype=int) + if non_corr: + extra = parse_orbitals(non_corr, ORBITALS) + if extra is None: + continue + projectors = projectors + extra + lines.append(" ".join(map(str, projectors))) + break + + irrep, to_write = None, ["0 0 0 0"] + if proj == "d": + irrep = ask("Split this orbital into it's irreps? (t2g/eg/n)\n", + ["t2g", "eg", "n"]) + if irrep == "t2g": + to_write = ["0 0 2 0", "01"] + elif irrep == "eg": + to_write = ["0 0 2 0", "10"] + + if ask("Do you want to include soc? (y/n)\n", ["y", "n"]) == "y": + if proj == "d" and irrep in ("t2g", "eg"): + print("Warning: For SOC, dmftproj will use the entire " + "d-shell. Using entire d-shell!") + lines.append("0 0 0 0") + else: + lines.extend(to_write) + lines.append("1") + else: + lines.extend(to_write) + lines.append("0") + + else: # uncorrelated: projectors only + while True: + proj = input("Specify the projectors that you would like to " + "include? (type h for help)\n") + if proj == "h": + print("indicate orbital projectors using (s, p, d, or f). " + "For multiple, combine them (sp, pd, spd, etc.)") + continue + projectors = parse_orbitals(proj, ORBITALS) + if projectors is None: + continue + lines.append(" ".join(map(str, projectors))) + lines.append("0 0 0 0") + break + + while True: + line = parse_window(input("Specify the projection window around eF " + "(default unit is Ry, specify eV with " + "-X.XX X.XX eV)\n")) + if line is not None: + lines.append(line) + break + + return lines + + +def main(): + case = os.path.basename(os.getcwd()) + struct_file = case + ".struct" + out_file = case + ".indmftpr" + + if "-h" in sys.argv or "--help" in sys.argv: + print(__doc__) + return 0 + if not os.path.isfile(struct_file): + print("Could not identify a {} file!".format(struct_file)) + return 1 + if os.path.isfile(out_file): + if ask("Previous {} detected! Continue? (y/n)\n".format(out_file), + ["y", "n"]) == "n": + return 0 + + print("Preparing dmftproj input file : {}\n".format(out_file)) + lines = collect_indmftpr(struct_file) + + # Written only once every answer is in, so an abort leaves the previous file + # (or no file) rather than a truncated one. + with open(out_file, "w") as out: + out.write('\n'.join(lines) + '\n') + print("initialize {} file ok!".format(out_file)) + return 0 + + +if __name__ == "__main__": + try: + sys.exit(main()) + except (KeyboardInterrupt, EOFError): + print("\naborted; no file written") + sys.exit(130)