Phonopy band unfolding tools

0-toc.md

Table of contents:

• unfolded-band-string.py
• merge-yaml.py

merge-yaml.py

#!/usr/bin/env python
from __future__ import print_function

try:
    from StringIO import StringIO
except ImportError:
    from io import StringIO

import sys
import yaml

def main():
    from argparse import ArgumentParser
    parser = ArgumentParser(
        description="Transform a supercell band.yaml to have the unit cell's qpoint info")
    parser.add_argument("PRIM", help='primitive band.yaml')
    parser.add_argument("SUPER", help='supercell band.yaml')
    args = parser.parse_args()

    if sys.stdout.isatty():
        info("NOTE: Output will be written to terminal.")
        info("      You probably want to redirect it. (Ctrl-C to cancel)")

    # simple error checking before continuing with bigger work
    info("performing sanity check...")
    p = read_only_nqpoint(args.PRIM)
    s = read_only_nqpoint(args.SUPER)
    if p != s:
        parser.error('incompatible number of qpoints: {} vs {}'.format(p, s))

    info("splitting primitive yaml...")
    (prim_pre, prim_rest) = partial_parse_band_yaml(open(args.PRIM))
    info("splitting supercell yaml...")
    (super_pre, super_rest) = partial_parse_band_yaml(open(args.SUPER))

    info("merging header portion of yaml files...")
    pre = merge_by_schema(prim_pre, super_pre, SCHEMA)
    info("merging phonon portion of yaml files...")
    rest = merge_phonon_section(prim_rest, super_rest)

    def write_to(f):
        yaml.dump(pre, f)
        print(file=f)
        f.writelines(rest)

    info("writing output")
    write_to(sys.stdout)

#-----------------------------------------------------
# Fast sanity checks
    
def read_only_nqpoint(path):
    with open(path) as f:
        # expect nqpoint before phonon:
        for line in f:
            prefix = 'nqpoint:'
            if line.startswith(prefix):
                return int(line[len(prefix):].strip())

            if line.startswith('phonon'):
                info("Cannot find nqpoint in {}".format(path))
                sys.exit(1)

#-----------------------------------------------------
# Splitting the yaml file

# use a yaml parser on everything except the 'phonon' section
def partial_parse_band_yaml(f):
    lines = list(f)

    # NOTE: this is slower than necessary as it scans the whole
    #       file but it helps ensure that we don't miss anything
    #       new that might have gotten added after "phonon:"
    pairs = list(yaml_outer_mapping_line_indices(lines))
    assert pairs[-1][1] == "phonon"

    index = pairs[-1][0]
    initial = lines[:index]
    initial = yaml.load(StringIO(''.join(initial)))
    rest = lines[index:]
    return (initial, rest)

# Identify the lines that contain keys for the outermost mapping
#  in the yaml file.
def yaml_outer_mapping_line_indices(lines):
    import string
    alpha = set(string.ascii_letters)
    for (i, line) in enumerate(lines):

        if line[0] in alpha:
            if ':' not in line:
                info('strange line: ' + repr(line))
                info('not sure what to do, aborting.')
                sys.exit(1)
            key = line[:line.index(':')]
            yield (i, key)

#-----------------------------------------------------
# merging the header portion

USE_PRIM = object()
USE_SUPER = object()
SCHEMA = {
    # NOTE: At the time of writing this (v 1.11.14)
    #       bandplot does not use any of the following.
    'nqpoint': USE_PRIM,
    'npath': USE_PRIM,
    'reciprocal_lattice': USE_PRIM,
    'natom': USE_SUPER,
    'lattice': USE_PRIM,
    'points': USE_SUPER,
    'supercell_matrix': USE_SUPER,

    # NOTE: bandplot DOES use the following
    'segment_nqpoint': USE_PRIM,
}

def merge_by_schema(prim, supercell, schema):
    # base cases
    if schema is USE_PRIM:
        return prim
    elif schema is USE_SUPER:
        return supercell

    # recursive case
    elif isinstance(schema, dict):
        assert isinstance(prim, dict)
        assert isinstance(supercell, dict)

        prim = dict(prim)
        supercell = dict(supercell)
        out = {}
        for k in schema:
            out[k] = merge_by_schema(prim.pop(k), supercell.pop(k), schema[k])

        def expect_empty(d):
            if d:
                print("warning: unexpected keys in yaml file:", file=sys.stderr)
                for key in d:
                    print(" - {!r}".format(key), file=sys.stderr)

        expect_empty(prim)
        expect_empty(supercell)
        return out
    
    else: assert False, "bug: incomplete switch"

#-----------------------------------------------------
# merging the phonon sectino

def merge_phonon_section(primitive, supercell):
    # the less copying necessary, the better
    assert isinstance(primitive, list)
    assert isinstance(supercell, list)

    prefixes = [
        '  distance:',
        '- q-position:',
    ]
    should_use_primitive = lambda s: any(s.startswith(pre) for pre in prefixes)

    prim_ids = list(all_indices(primitive, should_use_primitive))
    super_ids = list(all_indices(supercell, should_use_primitive))

    # sanity check: make sure at least one of each prefix was found
    # (in case the fields were renamed or indentation adjusted)
    for pre in prefixes:
        assert any(primitive[i].startswith(pre) for i in prim_ids[:2])
        assert any(supercell[i].startswith(pre) for i in super_ids[:2])

    assert len(prim_ids) == len(super_ids)
    for (p, s) in zip(prim_ids, super_ids):
        supercell[s] = primitive[p]

    return supercell

#-----------------------------------------------------
# utils

def all_indices(it, pred):
    for (i, x) in enumerate(it):
        if pred(x):
            yield i

# stdout is our output file, so put diagnostics on stderr
info = lambda *args: print(*args, file=sys.stderr)

#-----------------------------------------------------

if __name__ == '__main__':
    main()

unfolded-band-string.py

from __future__ import print_function
import numpy as np

# NOTE: This has been updated to fix the final k-point

#------------------------------
# edit this to match the params given to phonopy for the unit cell

# sample density
BAND_POINTS = 200
# high symmetry path
BAND = [
    [0.0, 0.0, 0.0],
    [0.5, 0.0, 0.0],
    [1./3., 1./3., 0.0],
    [0.0, 0.0, 0.0],
]

# supercell dimensions
DIMS = [6, 6, 1]
#------------------------------

def unit_qs(path, density):
    out = []
    for q1, q2 in zip(path, path[1:]):
        out.append(np.array([
            np.linspace(c1, c2, density)
            for (c1, c2) in zip(q1, q2)
        ]).T)
        
    # add one more arbitrary qpoint to the very end, which will be ignored
    #  if we are using BAND_POINTS = 1.
    out.append([[0, 0, 0]])
    return np.vstack(out)
unit_qs = unit_qs(BAND, BAND_POINTS)

def super_qs(unit_qs, sc_matrix):
    # convert from unitcell-reciprocal-fractional coordinates
    # into supercell-reciprocal-fractional
    qs = np.dot(unit_qs, sc_matrix.T)
    # reduce into the first brillouin zone
    # (actually, not quite into the FBZ, but close enough)
    qs %= 1
    return qs

super_qs = super_qs(unit_qs, np.diag(DIMS))

# leave some blank space so things can easily be added to the top
print()
print()
print("BAND_POINTS = 1")
print("BAND = ", " ".join(str(x) for x in super_qs.flat))