How to Plot Band Plots & VDOS Plots

General Steps

Band Plot

  1. pw.x
  2. ph.x
  3. q2r.x
  4. matdyn.x

VDOS

  1. pw.x
  2. ph.x
  3. q2r.x
  4. matdyn.x

Scripts

Bash Scipt (plotband.sh)

Run this script on supercomputer

Remember to edit:

  • amass in matdyn input
  • nk in matdyn input
#!/bin/sh

# run from directory where this script is
# cd `echo $0 | sed 's/\(.*\)\/.*/\1/'` # extract pathname
WORKING_DIR=`pwd -LP`

# check whether echo has the -e option
if test "`echo -e`" = "-e" ; then ECHO=echo ; else ECHO="echo -e" ; fi

$ECHO
$ECHO "GET THAT FUCKING BAND PLOT AND VDOS!"
$ECHO "===================================="
date
$ECHO
$ECHO ": ${WORKING_DIR}: starting"
$ECHO 

# how to run executables
PW_COMMAND="pw.x"
PH_COMMAND="ph.x"
MATDYN_COMMAND="matdyn.x"
Q2R_COMMAND="q2r.x"
PLOTBAND_COMMAND="plotband.x"
GP_COMMAND="gnuplot"

$ECHO "  running pw.x as:     $PW_COMMAND"
$ECHO "  running ph.x as:     $PH_COMMAND"
$ECHO "  running q2r.x as:    $Q2R_COMMAND"
$ECHO "  running matdyn.x as: $MATDYN_COMMAND"
$ECHO "  running plotband.x as: $PLOTBAND_COMMAND"
$ECHO "  running gnuplot as: $GP_COMMAND"
$ECHO

CALC_PREFIX="AlOOH"

for PHONON_SCF_DIR in `ls | grep 'scf-*[0-9]'`
#for PHONON_SCF_DIR in `ls | grep 'scf10'`
do

cd $PHONON_SCF_DIR

$ECHO ": entering directory: `pwd -P`"
$ECHO

cat > q2r-bandplot.in <<EOF
 &input
   fildyn='dyn'
   zasr='simple'
   flfrc='${CALC_PREFIX}-bandplot.fc'
 /
 8,8,8
`ls | grep '^dyn[1-9]' | wc -l`
`ls | grep '^dyn[1-9]'`
EOF

$ECHO "  transforming C(q) => C(R)...\c"
$ECHO "  $Q2R_COMMAND < q2r-bandplot.in > q2r-bandplot.out...\c"
$Q2R_COMMAND < q2r-bandplot.in > q2r-bandplot.out
$ECHO " done"

cat > matdyn-bandplot.in <<EOF
 &input
    asr='simple'
    amass(1) = 26.98154
    amass(2) = 15.99900
    amass(3) = 1.00800
    flfrc='${CALC_PREFIX}-bandplot.fc'
    flfrq='${CALC_PREFIX}-bandplot.freq'
    q_in_band_form=.true.
    q_in_cryst_coord = .true.
 /
 10
    0.00000000  0.00000000  0.00000000    100.0 # gG
    0.50000000  0.00000000  0.00000000    100.0 # Z
    0.50000000  0.50000000  0.00000000    100.0 # T
    0.00000000  0.50000000  0.00000000    100.0 # Y
    0.00000000  0.00000000  0.00000000    100.0 # gG
    0.00000000  0.00000000  0.50000000    100.0 # X
    0.50000000  0.00000000  0.50000000    100.0 # U
    0.50000000  0.50000000  0.50000000    100.0 # R
    0.00000000  0.50000000  0.50000000    100.0 # S
    0.00000000  0.00000000  0.50000000    1.0   # X
EOF

$ECHO "  recalculating omega(q) from C(R)...\c"
$ECHO "  $MATDYN_COMMAND < matdyn-bandplot.in > matdyn-bandplot.out...\c"
$MATDYN_COMMAND < matdyn-bandplot.in > matdyn-bandplot.out
$ECHO "  done"

cat > plotband.in <<EOF
${CALC_PREFIX}-bandplot.freq
0 2500
${CALC_PREFIX}-bandplot.plot
${CALC_PREFIX}-bandplot.ps
0.0
50.0 0.0
EOF

$ECHO "  writing the phonon dispersions in ${CALC_PREFIX}-bandplot.plot...\c"
$ECHO "  $PLOTBAND_COMMAND < plotband.in > plotband.out...\c"
$PLOTBAND_COMMAND < plotband.in > plotband.out
$ECHO "  done"

cat > matdyn-dos.in <<EOF
 &input
    asr='simple'
    amass(1) = 26.98154
    amass(2) = 15.99900
    amass(3) = 1.00800
    dos=.true.
    flfrc='${CALC_PREFIX}-bandplot.fc'
    fldos='${CALC_PREFIX}-dos.dos'
    flfrq='${CALC_PREFIX}-dos.freq'
    q_in_band_form=.true.
    nk1=8
    nk2=8
    nk3=8
 /
EOF

$ECHO "  recalculating rho(omega) from C(R)... (DOS)\c"
$ECHO "  $MATDYN_COMMAND < matdyn-dos.in > matdyn-dos.out...\c"
$MATDYN_COMMAND < matdyn-dos.in > matdyn-dos.out
$ECHO "  done"
$ECHO "  exit scf directory"
$ECHO ""

cd ..

done

$ECHO
$ECHO ": collecting DOS and band data..."
$ECHO

python3 find_dos.py

$ECHO "  output writting to bands.json"
$ECHO

Python (find_dos.py)

This script will be called by the plotband.sh.

import os
import io
import glob
import json
import numpy

prefix = 'AlOOH'
kpts_pos = [
    0.0000,
    0.5000,
    1.0604,
    1.5604,
    2.1209,
    2.9487,
    3.4487,
    4.0091,
    4.5091,
    5.0696,
]
kpts_labels = ["$\Gamma$", "Z", "T", "Y", "$\Gamma$", "X", "U", "R", "S", "X"]

def read_kpts_pos(fp: io.TextIOBase):
    import re
    for line in fp:
        result = re.search("x coordinate\s+([0-9\.]+)", line.strip())
        if result is not None:
            yield float(result.group(1))

def read_band(fp: io.TextIOBase):
    prev_ptr = 0
    while True:
        line = fp.readline().strip()
        if not line:
            curr_ptr = fp.tell()
            if curr_ptr == prev_ptr: break
            fp.seek(prev_ptr)
            chunk = fp.read(curr_ptr - prev_ptr)
            yield numpy.loadtxt(io.StringIO(chunk)).transpose()
            prev_ptr = fp.tell()

def read_dos(dirname):
    return numpy.loadtxt(os.path.join(dirname, '%s-dos.dos' % prefix)).transpose()

def read_xmgr(dirname):
    with open(os.path.join(dirname, '%s-bandplot.plot' % prefix)) as fp:
        return list(read_band(fp))

def get_opts(work_dirs):
    for work_dir in work_dirs:
        yield {
                'work_dir': work_dir,
                # 'title': '%s GPa' % work_dir.split('_')[1],
                'title': '%s GPa' % work_dir.lstrip('scf'),
                'bands': [band.tolist() for band in read_xmgr(work_dir)],
                'dos': read_dos(work_dir).tolist(),
                'kpts_pos': kpts_pos,
                'kpts_labels': kpts_labels
        } 

if __name__ == '__main__':
    dos_dirs =  [os.path.dirname(filename) for filename in glob.glob(r'*/%s-dos.dos'  % prefix)]
    xmgr_dirs = [os.path.dirname(filename) for filename in glob.glob(r'*/%s-bandplot.plot' % prefix)]
    
    print("  * finding bandplot data..")
    print("  * DOS directory: ", dos_dirs)
    print("  * band plot directory: ", xmgr_dirs)
    
    common_dirs = [dirname for dirname in dos_dirs if dirname in xmgr_dirs]
    
    with open(glob.glob(r'*/plotband.out')[0]) as fp:
        kpts_pos = list(read_kpts_pos(fp))
    
    print("  * k-point x coords: ", kpts_pos)
    
    if len(kpts_pos) != len(kpts_labels):
        print ("  ! k-point labels might be wrong, please double check!")
    
    with open('bands.json', 'w') as fp:
        json.dump(list(get_opts(common_dirs)), fp)

Python

Get bands.json and run this script on local computer.

import json
import os
import numpy
import io
import json

import matplotlib
import matplotlib.pyplot as plt

prefix = 'AlOOH'
kpts_pos = [0.0000, 0.0558, 0.1184, 0.1742, 0.2368, 0.3293 ,0.3851, 0.4476, 0.5035]
kpts_labels = ["$\Gamma$", "Z", "T", "Y", "$\Gamma$", "X", "U", "S", "X"]

def read_band(fp: io.TextIOBase):
    prev_ptr = 0
    while True:
        line = fp.readline().strip()
        if not line:
            curr_ptr = fp.tell()
            if curr_ptr == prev_ptr: break
            fp.seek(prev_ptr)
            chunk = fp.read(curr_ptr - prev_ptr)
            yield numpy.loadtxt(io.StringIO(chunk)).transpose()
            prev_ptr = fp.tell()

def read_dos(dirname):
    return numpy.loadtxt(os.path.join(dirname, '%s.dos' % prefix)).transpose()

def read_xmgr(dirname):
    with open(os.path.join(dirname, '%s.xmgr' % prefix)) as fp:
        return list(read_band(fp))

def plot_band(fig: matplotlib.figure.Figure,
              ax_bands: matplotlib.axes.Axes,
              ax_dos: matplotlib.axes.Axes,
              bands: list, kpts_pos, kpts_labels, dos, title=None):
    for band_q, band_freq in bands:
        ax_bands.plot(band_q, band_freq, c='k', alpha=.5)
    for kpt_pos in kpts_pos:
        ax_bands.axvline(kpt_pos, c='k', lw=.5)
    ax_bands.set_ylabel(r'Frequency')
    ax_bands.set_xticks(kpts_pos)
    ax_bands.set_xticklabels(kpts_labels)
    ax_bands.set_xlim((numpy.min(kpts_pos), numpy.max(kpts_pos)))

    dos_freq, dos_density = dos
    ax_dos.plot(dos_density, dos_freq, c='k', lw=1)
    ax_dos.set_xticks([])
    ax_dos.set_yticks([])

    fig.tight_layout(rect=[0, 0, 1, .95])
    if title: fig.suptitle(title, va='top')

def plot_band_dir(dirname):
    fig = plt.figure()
    (ax_bands, ax_dos) = fig.subplots(1, 2, gridspec_kw = {
        'width_ratios': (5, 1)
    })
    plot_band(fig, ax_bands, ax_dos, read_xmgr(dirname), kpts_pos, kpts_labels, read_dos(dirname))
    fig.savefig('%s.pdf' % dirname)

def plot_band_option(option, out_dir='.'):
    fig = plt.figure()
    (ax_bands, ax_dos) = fig.subplots(1, 2, gridspec_kw = {
        'width_ratios': (5, 1)
    })
    bands = [numpy.array(band) for band in option['bands']]
    dos = numpy.array(option['dos'])
    kpts_labels = option['kpts_labels']
    kpts_pos = option['kpts_pos']
    title = option['title']
    plot_band(fig, ax_bands, ax_dos, bands, kpts_pos, kpts_labels, dos, title=title)
    fig.savefig(os.path.join(out_dir, '%s.png' % option['work_dir']), dpi=300)

if __name__ == '__main__':
    with open('bands.json') as fp:
        options = json.load(fp)
    for option in options:
        plot_band_option(option, 'plots')