PAOFLOW.gen.paoflow_driver

Interactive generator for PAOFLOW main.py driver scripts.

Given the output of a Quantum ESPRESSO run (a <prefix>.save directory) and the pseudopotential used for it, this CLI asks a few questions and writes a minimal, clearly-commented main.py that runs PAOFLOW.

Two workflows are supported:

1. ACBN0 / eACBN0 self-consistent Hubbard U (and intersite V) using the standard basis for the PAO projections.

2. Full property run from the QE output, where the user picks the properties to compute from a menu. Standard properties use the standard basis; optical properties (dielectric tensor) need the extended basis and the non-local velocity correction, so they are emitted as a separate PAOFLOW run in the same script.

Two further workflows generate multi-phase driver scripts:

3. Harmonic phonons (main.phonon.py) – finite-displacement phonons via phonopy (generate supercells -> run pw.x -> analyse forces).

4. Electron-phonon (main.elphon.py) – the pseudo-atomic-orbital (Agapito & Bernardi) interpolation of QE’s DFPT coupling into alpha^2F / lambda / Tc (write the ph.x phonon + AHC inputs -> run QE -> analyse).

The generated script is static and heavily commented so it is easy to tweak afterwards.

Attributes

HEADER
BASIS_GEN_BLOCK
ENERGY_RANGE_BLOCK
PROPERTY_MENU
SPIN_PROPERTIES
ENERGY_WINDOW_PROPERTIES
ELPHON_TEMPLATE
ELPHON_PLOT_TEMPLATE
PLOT_HEADER

Functions

ask(prompt[, default])	Ask for a free-text value with an optional default.
ask_yes_no(prompt[, default])	Ask a yes/no question.
ask_int(prompt, default)	Ask for an integer value with a default.
ask_float(prompt, default)	Ask for a floating-point value with a default.
ask_choice(prompt, choices, default)	Ask the user to pick one item from choices (list of strings).
parse_laser_list(spec)	Parse a laser-wavelength specification into a list of floats (nm).
detect_savedir(workdir)
detect_upfs(workdir[, savedir])
detect_prefix(workdir, savedir)
detect_v_cutoff(workdir)	Read the suggested eACBN0 intersite-V cutoff from a QE input header.
select_properties()	Show the property menu and return the ordered set of chosen keys.
build_run_script(cfg)	Assemble a full property-run main.py from the collected config.
build_acbn0_script(cfg)	Assemble an ACBN0 / eACBN0 main.py from the collected config.
build_phonon_script(cfg)	Assemble a main.phonon.py harmonic-phonon workflow from the config.
build_elphon_script(cfg)	Assemble a main.elphon.py electron-phonon (PAO route) workflow from the config.
build_elphon_plot_script(cfg)	Assemble a plot.elphon.py for the Eliashberg alpha^2F / lambda output.
build_plot_script(cfg)	Assemble a plot.py that mimics the property selection in cfg.
build_acbn0_plot_script(cfg)	Assemble a plot.acbn0.py that overlays the ACBN0 / eACBN0 band structures.
build_phonon_plot_script(cfg)	Assemble a plot.phonon.py that plots the dispersion, DOS and thermals.
build_raman_script(cfg)	Assemble a main.raman.py non-resonant (Placzek) Raman workflow.
build_raman_plot_script(cfg)	Assemble a plot.raman.py that overlays the broadened Raman spectra.
collect_common(args, workdir)	Prompt for configuration shared by both workflows.
collect_run(common)	Prompt for the full property-run configuration.
collect_phonon(common)	Prompt for the harmonic-phonon (finite-displacement) configuration.
collect_elphon(common)	Prompt for the electron-phonon (PAO route) configuration.
collect_acbn0(common)	Prompt for the ACBN0 / eACBN0 configuration.
main([argv])

Module Contents

PAOFLOW.gen.paoflow_driver.ask(prompt, default=None)

Ask for a free-text value with an optional default.

PAOFLOW.gen.paoflow_driver.ask_yes_no(prompt, default=False)

Ask a yes/no question.

PAOFLOW.gen.paoflow_driver.ask_int(prompt, default)

Ask for an integer value with a default.

PAOFLOW.gen.paoflow_driver.ask_float(prompt, default)

Ask for a floating-point value with a default.

PAOFLOW.gen.paoflow_driver.ask_choice(prompt, choices, default)

Ask the user to pick one item from choices (list of strings).

PAOFLOW.gen.paoflow_driver.parse_laser_list(spec)

Parse a laser-wavelength specification into a list of floats (nm).

Accepts a comma- or space-separated list ('488, 514.5, 532') or a restricted Python expression that evaluates to an iterable of numbers, e.g. '[n for n in range(450, 650, 5)]'. The expression is evaluated with no builtins and only range exposed, so it cannot import modules or call arbitrary functions.

PAOFLOW.gen.paoflow_driver.detect_savedir(workdir)

PAOFLOW.gen.paoflow_driver.detect_upfs(workdir, savedir=None)

PAOFLOW.gen.paoflow_driver.detect_prefix(workdir, savedir)

PAOFLOW.gen.paoflow_driver.detect_v_cutoff(workdir)

Read the suggested eACBN0 intersite-V cutoff from a QE input header.

paoflow-gen-qe writes ! Suggested eACBN0 intersite V cutoff: <x> Angstrom into the header of the <compound>.scf.in it generates. Returns that cutoff (in Angstrom) from the first matching input file, or None if no input file or comment is found.

PAOFLOW.gen.paoflow_driver.HEADER = Multiline-String

Show Value

"""#!/usr/bin/env python3
"""PAOFLOW driver script (generated by paoflow_gen.py).

Edit the constants below and the property calls in the run function(s) to
customize the calculation.  Run with:

    python main.py
    # or, in parallel:
    mpirun -np <N> python main.py
"""

import os
import sys

from PAOFLOW import PAOFLOW
from PAOFLOW.basis_gen import generate_basis_for_pseudo
from PAOFLOW.basis_gen.driver import _default_shells
from PAOFLOW.inputs.read_upf import UPF as _UPFParser

try:
    from mpi4py import MPI

    RANK = MPI.COMM_WORLD.Get_rank()
except ImportError:
    RANK = 0
"""

PAOFLOW.gen.paoflow_driver.BASIS_GEN_BLOCK = Multiline-String

Show Value

"""
def ensure_basis(preset="extended"):
    """Generate the pseudo-atom basis under BASISPATH for every species.

    The 'extended' preset is a superset of 'standard' and 'minimal', so
    generating it once is enough for every configuration used below.
    """
    if RANK != 0:
        return
    for upf_path in UPFS:
        upf = _UPFParser(upf_path)
        element = upf.element.strip()
        elem_dir = os.path.join(BASISPATH, element)
        expected = _default_shells(upf, preset=preset)
        missing = [
            s for s in expected
            if not os.path.exists(os.path.join(elem_dir, "{}.dat".format(s)))
        ]
        if missing:
            print("Generating pseudo-atom basis for {} under {} ...".format(
                element, BASISPATH))
            generate_basis_for_pseudo(
                upf_path, BASISPATH.rstrip(os.sep), preset=preset, verbose=True
            )
        else:
            print("Using existing pseudo-atom basis for {} under {}".format(
                element, BASISPATH))
"""

PAOFLOW.gen.paoflow_driver.ENERGY_RANGE_BLOCK = Multiline-String

Show Value

"""
def suggest_energy_window(p):
    """Print the full PAO band range (eV, rel. to E_F) as a suggested window.

    Must be called after p.pao_eigh().  The eigenvalues ('E_k') are distributed
    across MPI ranks, so the local extrema are reduced to give every process the
    same global range.  This only prints a hint; the calculation uses the
    user-defined EMIN / EMAX / NE constants above.
    """
    import numpy as np

    E_k = p.data_controller.data_arrays['E_k']
    emin = float(np.amin(E_k))
    emax = float(np.amax(E_k))
    if "MPI" in globals():
        comm = MPI.COMM_WORLD
        emin = comm.allreduce(emin, op=MPI.MIN)
        emax = comm.allreduce(emax, op=MPI.MAX)
    if RANK == 0:
        print('Suggested energy window from PAO bands: '
              'EMIN={:.4f}, EMAX={:.4f} eV'.format(emin, emax))
        print('  (currently using EMIN={}, EMAX={}, NE={})'.format(EMIN, EMAX, NE))
    return emin, emax
"""

PAOFLOW.gen.paoflow_driver.PROPERTY_MENU = [('bands', 'Band structure'), ('dos', 'DOS / projected DOS'), ('transport', 'Boltzmann transport...

PAOFLOW.gen.paoflow_driver.SPIN_PROPERTIES

PAOFLOW.gen.paoflow_driver.ENERGY_WINDOW_PROPERTIES

PAOFLOW.gen.paoflow_driver.select_properties()

Show the property menu and return the ordered set of chosen keys.

PAOFLOW.gen.paoflow_driver.build_run_script(cfg)

Assemble a full property-run main.py from the collected config.

PAOFLOW.gen.paoflow_driver.build_acbn0_script(cfg)

Assemble an ACBN0 / eACBN0 main.py from the collected config.

PAOFLOW.gen.paoflow_driver.build_phonon_script(cfg)

Assemble a main.phonon.py harmonic-phonon workflow from the config.

The generated script drives the finite-displacement phonon workflow in three phases (generate displaced supercells -> run pw.x -> analyse forces), which can be run together or one at a time (handy on HPC schedulers).

PAOFLOW.gen.paoflow_driver.ELPHON_TEMPLATE = Multiline-String

Show Value

"""#!/usr/bin/env python3
"""PAOFLOW electron-phonon (Eliashberg) workflow -- PAO route (generated by paoflow_gen.py).

Pseudo-atomic-orbital (Agapito & Bernardi, Phys. Rev. B 97, 235146 (2018)) interpolation
of Quantum ESPRESSO's DFPT electron-phonon coupling.  PAOFLOW reads QE's *full*
coarse-grid coupling (no potential reconstruction), rotates it into the PAO gauge
and Wigner-Seitz interpolates electrons + vertex to a dense grid for alpha^2F,
lambda, omega_log and Tc.

Two phases:

    python main.elphon.py inputs    # write the ph.x phonon + AHC input templates
    python main.elphon.py analyse   # PAO interpolation -> alpha^2F, lambda, Tc

The ``analyse`` phase parallelises the per-q interpolation over MPI ranks, so
on a cluster launch it with mpirun for an (up to nq-fold) speedup::

    mpirun -np N python main.elphon.py analyse

The dense electron diagonalisation is done redundantly on every rank, so it is
usually best to keep N at or below the number of q-points and let BLAS threads
(OMP_NUM_THREADS) fill the remaining cores.

Between them, run the two QE ph.x steps in the same outdir (typically on HPC):

    1. phonon:  ph.x < <prefix>.ph.in    # full DFPT dvscf on the q-grid
    2. ahc:     ph.x < <prefix>.ahc.in   # electron_phonon='ahc' -> ahc_dir/

``fildvscf`` and ``fildyn`` MUST match between the two ph.x steps.  The AHC path
(SOURCE='ahc') is for norm-conserving pseudopotentials; for ultrasoft / PAW use
the patched-QE ``el_ph_mat`` dump (SOURCE='elphmat').
"""

import argparse
import os
import sys

import numpy as np
from mpi4py import MPI

from PAOFLOW import PAOFLOW
from PAOFLOW.elphon.do_pao_eph import eliashberg_from_qe_coupling
from PAOFLOW.elphon.elph_bloch import read_nscf


# ----------------------------------------------------------------------- #
# Configuration  (edit freely -- masses / NELEC / NBND are system-specific) #
# ----------------------------------------------------------------------- #
HERE = os.path.dirname(os.path.abspath(__file__))
SAVEDIR = os.path.join(HERE, __SAVEDIR__)
OUTPUTDIR = __OUTPUTDIR__
PREFIX = __PREFIX__
BASISDIR = os.path.join(HERE, __BASISDIR__)

# Coupling source:
#   'ahc'     -> unpatched QE AHC dumps (ahc_dir/ahc_gkk_iq<iq>.bin); NC pseudos.
#   'elphmat' -> patched-QE el_ph_mat dumps (elph_dir/elphmat.<iq>.dat); any pseudo.
SOURCE = __SOURCE__
COUPLING_DIR = os.path.join(HERE, __COUPLING_DIR__)

KGRID = __KGRID__          # SCF / coupling k-grid (== pw.x K_POINTS)
QGRID = __QGRID__          # phonon q-grid (nq1, nq2, nq3)
NBND = __NBND__            # bands in the nscf / AHC run (nbnd = ahc_nbnd)
MASSES_AMU = __MASSES__    # atomic masses (amu), one per atom in the cell
NELEC = __NELEC__          # valence electrons (dense E_F recompute)
NK_DENSE = __NK_DENSE__    # dense interpolation grid
SIGMA_RY = __SIGMA__       # Fermi-surface smearing (Ry)
MU_STAR = __MU_STAR__      # Coulomb pseudopotential for Tc
PTHR = __PTHR__            # projectability threshold
DYNPREFIX = PREFIX         # <DYNPREFIX>.dyn<iq> phonon files

# For SOURCE='elphmat' (patched, symmetry-reduced) set the irreducible-q star
# weights here; leave empty (or use SOURCE='ahc') to sum the full q-grid with
# unit weights.
Q_WEIGHTS = __QWEIGHTS__


def _phonon_input():
    return '\n'.join([
        '&inputph',
        "  prefix='%s'," % PREFIX,
        "  outdir='./',",
        "  fildyn='%s.dyn'," % PREFIX,
        "  fildvscf='dvscf',",
        '  tr2_ph=1.0d-12,',
        '  ldisp=.true., nq1=%d, nq2=%d, nq3=%d,' % (QGRID[0], QGRID[1], QGRID[2]),
        '/',
        '',
    ])


def _ahc_input():
    return '\n'.join([
        '&inputph',
        "  prefix='%s'," % PREFIX,
        "  outdir='./',",
        "  fildyn='%s.dyn'," % PREFIX,
        "  fildvscf='dvscf',",
        "  electron_phonon='ahc',",
        '  ahc_nbnd=%d, ahc_nbndskip=0, skip_upperfan=.true.,' % NBND,
        "  ahc_dir='ahc_dir/',",
        '  ldisp=.true., nq1=%d, nq2=%d, nq3=%d,' % (QGRID[0], QGRID[1], QGRID[2]),
        '/',
        '',
    ])


def inputs():
    """Phase 1: write the ph.x phonon and AHC input templates."""
    ph = os.path.join(HERE, PREFIX + '.ph.in')
    ahc = os.path.join(HERE, PREFIX + '.ahc.in')
    with open(ph, 'w') as fh:
        fh.write(_phonon_input())
    with open(ahc, 'w') as fh:
        fh.write(_ahc_input())
    print('Wrote %s' % ph)
    print('Wrote %s' % ahc)
    print('Run (same outdir), then `analyse`:')
    print('  1) ph.x < %s' % os.path.basename(ph))
    print('  2) ph.x < %s' % os.path.basename(ahc))


def analyse():
    """Phase 2: PAO interpolation of the QE coupling -> alpha^2F, lambda, Tc."""
    if not os.path.isdir(SAVEDIR):
        sys.exit('%s not found. Run pw.x (nscf) first.' % SAVEDIR)
    if not os.path.isdir(COUPLING_DIR):
        sys.exit('%s not found. Run the QE phonon + AHC steps first (phase: inputs).' % COUPLING_DIR)

    pf = PAOFLOW.PAOFLOW(workpath=HERE, outputdir=OUTPUTDIR, savedir=SAVEDIR, verbose=False)
    pf.projections(configuration='standard', basispath=BASISDIR)
    pf.projectability(pthr=PTHR)
    # Grab the projection matrices A_k BEFORE pao_hamiltonian (which deletes them).
    A = pf.data_controller.data_arrays['U'][:, :, :, 0].copy()
    pf.pao_hamiltonian()
    HRs = pf.data_controller.data_arrays['HRs']
    info = read_nscf(SAVEDIR)

    if SOURCE == 'ahc' or not Q_WEIGHTS:
        nq = QGRID[0] * QGRID[1] * QGRID[2]
        q_weights = [1.0] * nq
    else:
        q_weights = list(Q_WEIGHTS)
        nq = len(q_weights)
    dyn_paths = [os.path.join(HERE, '%s.dyn%d' % (DYNPREFIX, i + 1)) for i in range(nq)]

    out = eliashberg_from_qe_coupling(
        A, HRs, info['kpts_cryst'], info['bg'], info['at'],
        COUPLING_DIR, q_weights, tuple(KGRID), dyn_paths,
        source=SOURCE, masses_amu=MASSES_AMU, nk_dense=NK_DENSE,
        sigmas_ry=[SIGMA_RY], nelec=NELEC, mu_star=MU_STAR,
    )

    # The Eliashberg result is identical on every rank (Allreduce inside the
    # driver); only rank 0 reports and writes the output files.
    if MPI.COMM_WORLD.Get_rank() != 0:
        return

    kB = 8.617333262e-5  # eV/K
    print('PAO-route Eliashberg (%s, source=%s, Nk=%d):' % (PREFIX, SOURCE, NK_DENSE))
    print('  N(E_F)   = %.3f states/spin/Ry' % out['dos_ef'].mean())
    print('  lambda   = %.4f' % out['lambda'])
    print('  <w_log>  = %.1f K' % (out['omega_log'] / kB))
    print('  Tc (McM) = %.2f K  (mu* = %.3f)' % (out['Tc_mcmillan'], MU_STAR))
    print('  Tc (AD)  = %.2f K  (mu* = %.3f)' % (out['Tc_allen_dynes'], MU_STAR))

    outdir = os.path.join(HERE, OUTPUTDIR)
    os.makedirs(outdir, exist_ok=True)
    a2f = os.path.join(outdir, 'alpha2F.dat')
    np.savetxt(a2f, np.column_stack([out['omega'] * 1.0e3, out['a2F']]),
               header='omega(meV)   alpha^2F(omega)   '
                      '(lambda=%.4f, omega_log=%.2fK, Tc_McM=%.3fK, Tc_AD=%.3fK, mu*=%.3f)'
                      % (out['lambda'], out['omega_log'] / kB, out['Tc_mcmillan'],
                         out['Tc_allen_dynes'], MU_STAR))
    npz = os.path.join(outdir, 'eliashberg.npz')
    np.savez(npz, **out)
    print('  wrote %s' % a2f)
    print('  wrote %s' % npz)


def main():
    parser = argparse.ArgumentParser(description='PAOFLOW electron-phonon (PAO) workflow.')
    parser.add_argument('phase', nargs='?', choices=['inputs', 'analyse', 'all'],
                        default='all', help='Workflow phase to run (default: all).')
    args = parser.parse_args()
    if args.phase in ('inputs', 'all'):
        inputs()
    if args.phase in ('analyse', 'all'):
        analyse()


if __name__ == '__main__':
    main()
"""

PAOFLOW.gen.paoflow_driver.build_elphon_script(cfg)

Assemble a main.elphon.py electron-phonon (PAO route) workflow from the config.

PAOFLOW.gen.paoflow_driver.ELPHON_PLOT_TEMPLATE = Multiline-String

Show Value

"""#!/usr/bin/env python3
"""Plot the Eliashberg alpha^2F(omega) and cumulative lambda(omega).

    python plot.elphon.py

Reads OUTPUTDIR/eliashberg.npz written by main.elphon.py (analyse).
"""

import os

import numpy as np
import matplotlib.pyplot as plt

HERE = os.path.dirname(os.path.abspath(__file__))
OUTPUTDIR = __OUTPUTDIR__
NPZ = os.path.join(HERE, OUTPUTDIR, 'eliashberg.npz')


def main():
    if not os.path.isfile(NPZ):
        raise SystemExit('%s not found; run main.elphon.py analyse first.' % NPZ)
    d = np.load(NPZ)
    omega = d['omega'] * 1e3   # eV -> meV
    a2F = d['a2F']
    lam = float(d['lambda'])
    tc_ad = float(d['Tc_allen_dynes']) if 'Tc_allen_dynes' in d else None
    tc_mcm = float(d['Tc_mcmillan']) if 'Tc_mcmillan' in d else None
    mu = float(d['mu_star']) if 'mu_star' in d else None
    # Cumulative lambda(omega) = 2 * integral_0^omega a2F(w)/w dw.
    w = d['omega']
    with np.errstate(divide='ignore', invalid='ignore'):
        integrand = np.where(w > 0, 2.0 * a2F / w, 0.0)
    lam_cum = np.concatenate([[0.0], np.cumsum(0.5 * (integrand[1:] + integrand[:-1]) * np.diff(w))])

    fig, ax1 = plt.subplots(figsize=(6, 4))
    ax1.plot(omega, a2F, color='C0', label=r'$\alpha^2F(\omega)$')
    ax1.set_xlabel(r'$\omega$ (meV)')
    ax1.set_ylabel(r'$\alpha^2F(\omega)$', color='C0')
    ax1.set_xlim(left=0.0)
    ax1.set_ylim(bottom=0.0)
    ax2 = ax1.twinx()
    ax2.plot(omega, lam_cum, color='C3', label=r'$\lambda(\omega)$')
    ax2.set_ylabel(r'$\lambda(\omega)$', color='C3')
    ax2.set_ylim(bottom=0.0)
    title = r'Eliashberg spectral function ($\lambda = %.3f$)' % lam
    if tc_mcm is not None and tc_ad is not None:
        title += '\n' + r'$T_c^{McM} = %.2f$ K,  $T_c^{AD} = %.2f$ K ($\mu^* = %.2f$)' % (tc_mcm, tc_ad, mu)
    ax1.set_title(title)
    fig.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()
"""

PAOFLOW.gen.paoflow_driver.build_elphon_plot_script(cfg)

Assemble a plot.elphon.py for the Eliashberg alpha^2F / lambda output.

PAOFLOW.gen.paoflow_driver.PLOT_HEADER = Multiline-String

Show Value

"""#!/usr/bin/env python3
"""Plotting helper for the PAOFLOW run (generated by paoflow_gen.py).

This script mirrors the property selection made when the PAOFLOW driver was
generated.  Run it after the PAOFLOW calculation has produced its output files
and pick what to plot from the menu:

    python plot.py

The energy window and property y-axis limits can be set on the command line:

    python plot.py --emin -5 --emax 5 --ymin 0 --ymax 100

Pass ``--all`` to plot every available quantity without any prompting (the
menu and axis-limit prompts are skipped and the generated defaults are used):

    python plot.py --all

and are also asked for interactively when the script runs (press Enter to keep
the default / command-line value).  EMIN/EMAX set the energy axis (vertical for
band plots, horizontal for the energy-resolved property plots) and default to
the window identified when the driver was generated.  YMIN/YMAX set the property
y-axis and default to automatic scaling.  Optical spectra always start at 0 on
the energy axis.

The plots use the helpers in PAOFLOW.GPAO.
"""

import argparse
import glob
import os
import sys

from PAOFLOW import GPAO


HERE = os.path.dirname(os.path.abspath(__file__))
OUTPUTDIR = os.path.join(HERE, __OUTPUTDIR__)

pplt = GPAO.GPAO()

# Default energy window (eV, relative to E_F) identified when the driver was
# generated.  Overridable on the command line via --emin/--emax.  YMIN/YMAX
# default to None (automatic axis limits) and can be set via --ymin/--ymax.
EMIN = __EMIN__
EMAX = __EMAX__
YMIN = None
YMAX = None


def _ewin():
    """Energy-axis window (EMIN, EMAX) used for band / energy-resolved plots."""
    return (EMIN, EMAX)


def _ewin_optical():
    """Energy-axis window for optical spectra (always starts at 0)."""
    return (0.0, EMAX)


def _ylim():
    """Property-axis limits (YMIN, YMAX), or None for automatic scaling."""
    if YMIN is None and YMAX is None:
        return None
    return (YMIN, YMAX)


def _dos_mag_lim(fnames):
    """Magnitude-axis limits for DOS/PDOS plots.

    When the user sets YMIN/YMAX those win.  Otherwise the magnitude axis is
    rescaled to the data that falls inside the chosen energy window [EMIN, EMAX]
    so that limiting the energy range actually zooms the curve (important for
    metals, whose largest DOS peaks may sit far outside the window).
    """
    if YMIN is not None or YMAX is not None:
        return _ylim()
    if isinstance(fnames, str):
        fnames = [fnames]
    import numpy as np

    peak = 0.0
    for fn in fnames:
        if not fn:
            continue
        try:
            data = np.loadtxt(fn)
        except (OSError, ValueError):
            continue
        if data.ndim != 2 or data.shape[1] < 2:
            continue
        es = data[:, 0]
        mag = data[:, 1:]
        mask = (es >= EMIN) & (es <= EMAX)
        if not mask.any():
            continue
        peak = max(peak, float(np.nanmax(np.abs(mag[mask]))))
    if peak <= 0.0:
        return None
    return (0.0, 1.1 * peak)



def _ask_float(prompt, default):
    """Prompt for a float, accepting Enter to keep *default* (may be None)."""
    shown = 'auto' if default is None else default
    try:
        raw = input('{} [{}]: '.format(prompt, shown)).strip()
    except EOFError:
        return default
    if raw == '':
        return default
    if raw.lower() in ('auto', 'none'):
        return None
    try:
        return float(raw)
    except ValueError:
        print('  (invalid number, keeping {})'.format(shown))
        return default



def _one(pattern):
    """Return the first OUTPUTDIR file matching *pattern* (or None).

    Patterns are written with an optional '<prefix>.' in front (e.g.
    '*.bands_0.dat').  PAOFLOW may write the files either with that prefix
    ('Si.bands_0.dat') or without it ('bands_0.dat'), so a '*.' pattern
    falls back to the prefix-less form.
    """
    hits = sorted(glob.glob(os.path.join(OUTPUTDIR, pattern)))
    if not hits and pattern.startswith('*.'):
        hits = sorted(glob.glob(os.path.join(OUTPUTDIR, pattern[2:])))
    return hits[0] if hits else None


def _many(pattern):
    """Return all OUTPUTDIR files matching *pattern* (sorted).

    As with _one, a '*.' pattern also matches prefix-less files.
    """
    hits = sorted(glob.glob(os.path.join(OUTPUTDIR, pattern)))
    if not hits and pattern.startswith('*.'):
        hits = sorted(glob.glob(os.path.join(OUTPUTDIR, pattern[2:])))
    return hits


def _missing(*names):
    """Warn about missing output files; return True if anything is missing."""
    gone = [n for n in names if n is None]
    if gone:
        print("  (output file(s) not found in {}; run PAOFLOW first)".format(OUTPUTDIR))
        return True
    return False


def _file_energy_range(fname):
    """Return (min, max) of the energy column of *fname*, or None."""
    if not fname:
        return None
    import numpy as np

    try:
        data = np.loadtxt(fname)
    except (OSError, ValueError):
        return None
    if data.ndim != 2 or data.shape[0] == 0:
        return None
    es = data[:, 0]
    return (float(np.nanmin(es)), float(np.nanmax(es)))


def _default_energy_window():
    """Default energy window for the prompts.

    When a DOS file is present its energy grid is used so that the DOS fills
    the side-by-side panel (the DOS is usually computed on a narrower grid than
    the bands).  Otherwise the window identified when the driver was generated
    is kept.
    """
    return _file_energy_range(_one('*.dosdk_0.dat'))


_SPIN_COLORS = ['tab:blue', 'tab:red', 'tab:green', 'tab:orange']


def _spin_channels(pattern):
    """Return ([files], [labels]) for the spin channels of *pattern*.

    The generated patterns target the first spin channel through a '_0' tag.
    A spin-polarized (nspin=2) run also writes the '_1' channel; when present
    both files are returned so the two channels can be overlaid in one figure.
    """
    f0 = _one(pattern)
    if f0 is None:
        return [], []
    files, labels = [f0], ['spin up']
    f1 = _one(pattern.replace('_0.', '_1.'))
    if f1 is not None:
        files.append(f1)
        labels.append('spin down')
    return files, labels


def _overlay_transport(files, labels, title, y_label, scale=1.0, min_zero=False):
    """Overlay the diagonal-averaged transport tensor of several files.

    Used for spin-polarized conductivity / Seebeck output, where each spin
    channel is written to its own file: both channels are drawn on a single
    figure (one line per spin and temperature).
    """
    import numpy as np
    import matplotlib.pyplot as plt
    from PAOFLOW.inputs.read_pao_output import read_transport_PAO

    fig = plt.figure()
    fig.suptitle(title)
    ax = fig.add_subplot(111)
    ic = 0
    ymax = 0.0
    for fn, base in zip(files, labels):
        enes, temps, tensors = read_transport_PAO(fn)
        for it, temp in enumerate(temps):
            trace = scale * np.einsum('eii->e', tensors[it]) / 3.0
            lbl = base if len(temps) == 1 else '{}, T={:g}'.format(base, temp)
            ax.plot(enes, trace, color=_SPIN_COLORS[ic % len(_SPIN_COLORS)], label=lbl)
            ymax = max(ymax, float(np.nanmax(np.abs(trace))))
            ic += 1
    ax.set_xlim(*_ewin())
    yl = _ylim()
    if yl is not None:
        ax.set_ylim(*yl)
    elif min_zero and ymax > 0.0:
        ax.set_ylim(0.0, 1.1 * ymax)
    ax.set_xlabel('Energy (eV)')
    ax.set_ylabel(y_label)
    ax.legend()
    plt.show()


def _overlay_bands_dos(band_files, dos_files, labels, sym_file, title):
    """Bands beside DOS with both spin channels overlaid on one figure.

    Reproduces plot_dos_beside_bands but loops over the spin channels so
    that the two channels share a single bands panel and a single DOS panel.
    """
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib import gridspec
    from PAOFLOW.inputs.read_pao_output import (
        read_bands_PAO,
        read_dos_PAO,
        read_band_path_PAO,
    )

    sym_points = read_band_path_PAO(sym_file) if sym_file else None
    y_lim = _ewin()

    fig = plt.figure()
    spec = gridspec.GridSpec(ncols=2, nrows=1, width_ratios=[5, 1])
    fig.suptitle(title)
    ax_b = fig.add_subplot(spec[0])
    ax_d = fig.add_subplot(spec[1])

    nk = 1
    dos_peak = 0.0
    for i, (fb, fd, lbl) in enumerate(zip(band_files, dos_files, labels)):
        col = _SPIN_COLORS[i % len(_SPIN_COLORS)]
        bands = read_bands_PAO(fb)
        nk = bands.shape[1]
        for j, b in enumerate(bands):
            ax_b.plot(b, color=col, label=lbl if j == 0 else None)
        es, dos = read_dos_PAO(fd)
        ax_d.plot(dos, es, color=col, label=lbl)
        mask = (es >= y_lim[0]) & (es <= y_lim[1])
        if mask.any():
            dos_peak = max(dos_peak, float(np.nanmax(np.abs(dos[mask]))))

    ax_b.set_xlim(0, nk - 1)
    ax_b.set_ylim(*y_lim)
    if sym_points is None:
        ax_b.xaxis.set_visible(False)
    else:
        ax_b.set_xticks(sym_points[0])
        ax_b.set_xticklabels(sym_points[1])
        ax_b.vlines(sym_points[0], y_lim[0], y_lim[1], color='gray')
    ax_b.set_ylabel(r'$\epsilon$($\mathbf{k}$) (eV)', fontsize=12)
    ax_b.legend()

    xl = _ylim()
    if xl is not None:
        ax_d.set_xlim(*xl)
    elif dos_peak > 0.0:
        ax_d.set_xlim(0.0, 1.1 * dos_peak)
    ax_d.set_ylim(*y_lim)
    ax_d.yaxis.set_visible(False)
    ax_d.set_xlabel('DOS')

    plt.tight_layout()
    plt.show()
"""

PAOFLOW.gen.paoflow_driver.build_plot_script(cfg)

Assemble a plot.py that mimics the property selection in cfg.

PAOFLOW.gen.paoflow_driver.build_acbn0_plot_script(cfg)

Assemble a plot.acbn0.py that overlays the ACBN0 / eACBN0 band structures.

The ACBN0 driver writes one band-structure file per converged case (bands_U for the on-site-U solution and, for eACBN0, bands_UV for the joint U+V solution). This script overlays whichever of those cases are present so they can be compared directly on a single figure.

PAOFLOW.gen.paoflow_driver.build_phonon_plot_script(cfg)

Assemble a plot.phonon.py that plots the dispersion, DOS and thermals.

Reads the phonon_band.dat / phonon_dos.dat / phonon_band.labels / phonon_thermal.dat files written by main.phonon.py and renders them through PAOFLOW.GPAO.GPAO.

PAOFLOW.gen.paoflow_driver.build_raman_script(cfg)

Assemble a main.raman.py non-resonant (Placzek) Raman workflow.

The generated script displaces the primitive cell by +/-delta along every optical zone-centre eigenvector (read from the harmonic-phonon FORCE_SETS), writes an SCF-only QE input per displacement, and – in the analyse phase – runs the PAOFLOW internal-projection optical pipeline on each displaced cell to finite-difference the Raman tensor.

PAOFLOW.gen.paoflow_driver.build_raman_plot_script(cfg)

Assemble a plot.raman.py that overlays the broadened Raman spectra.

PAOFLOW.gen.paoflow_driver.collect_common(args, workdir)

Prompt for configuration shared by both workflows.

PAOFLOW.gen.paoflow_driver.collect_run(common)

Prompt for the full property-run configuration.

PAOFLOW.gen.paoflow_driver.collect_phonon(common)

Prompt for the harmonic-phonon (finite-displacement) configuration.

PAOFLOW.gen.paoflow_driver.collect_elphon(common)

Prompt for the electron-phonon (PAO route) configuration.

PAOFLOW.gen.paoflow_driver.collect_acbn0(common)

Prompt for the ACBN0 / eACBN0 configuration.

PAOFLOW.gen.paoflow_driver.main(argv=None)