How to run an IR Job

This gist tells you how to generate IR calulation's required input and how to run an IR job.

Make IR Input (make_ir.sh)

This script is based on VC and PHonon's input. Run this make_ir.sh after you have generated phonon's scf input.

The sum rule parameter (asr) need to be double checked.

#!/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 "MAKE THAT FUCKING IR!!!"
$ECHO "======================="
date
$ECHO
$ECHO ": ${WORKING_DIR}: starting"
$ECHO 

for PH_DIRNAME in `ls -d scf*/`
do
    $ECHO ": finding phonon dir ${PH_DIRNAME}"

    IR_DIRNAME=${PH_DIRNAME/scf/ir}
    
    $ECHO "  creating ${IR_DIRNAME} .."
        
    mkdir $IR_DIRNAME
    
    $ECHO "  copying ${PH_DIRNAME:0:-1}.in as scf.in ..\c"
    cp $PH_DIRNAME/${PH_DIRNAME:0:-1}.in $IR_DIRNAME/scf.in
    $ECHO " done"
    $ECHO "  entering ${IR_DIRNAME} .."
    
    cd $IR_DIRNAME
    
    $ECHO "  creating ph.in ..\c"
    cat > ph.in << EOF
Phonon
&INPUTPH
    tr2_ph=1.0d-14
    verbosity = 'high'
    epsil = .true.
    fildyn = 'dyn'
    outdir = './tmp'
    amass(1) = 26.98154
    amass(2) = 15.99900
    amass(3) = 1.00800
/
0.0 0.0 0.0
EOF
    $ECHO " done"

    $ECHO "  creating dynmat.in ..\c"
    cat > dynmat.in << EOF
&INPUT
    fildyn='dyn'
/
EOF
    $ECHO " done"
    cd ..
done

Submit a SCF + PHonon + DYNMAT Job

Remember to CHANGE THE HEADER!!

This part, you may need to check the number of IR you are running. The script will automatically run all IR jobs parallelly at the same time, but the job might not have so many processors. You can change number of processors you are using or set IR_DIRS manually, like

IR_DIR=(ir10 ir15 ir20)

This actual job file is (run-ir.sh), you should submit it (or you need to remove ibrun and run it by shell or bash):

#!/bin/sh
#SBATCH -p skx-dev # Queue (partition) name
#SBATCH -N 4               # Total # of nodes 
#SBATCH -n 24              # Total # of mpi tasks
#SBATCH -t 2:00:00        # Run time (hh:mm:ss)
#SBATCH --mail-user=xxxxxx@xxxxxx.edu
#SBATCH --mail-type=all    # Send email at begin and end of job
#SBATCH -A XXXXXX # Allocation name (req'd if you have more than 1)

# 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 "RUN THAT FUCKING IR!!!"
$ECHO "======================"
date
$ECHO
$ECHO ": ${WORKING_DIR}: starting"
$ECHO

IR_DIRS=`ls -d ir*/`

for IR_DIRNAME in $IR_DIRS
do
    $ECHO ": entering ${IR_DIRNAME}"
    cd $IR_DIRNAME
    $ECHO "  starting self-consistent calculation in ${IR_DIRNAME}..\c"
    ibrun pw.x -in scf.in > scf.out &
    $ECHO " done"
    sleep 2
    cd ..
    $ECHO ""
done
$ECHO "! (`date`) all self-consistent calculation has been submitted."

wait
$ECHO "! (`date`) all self-consistent calculation has been finished."

$ECHO "! (`date`) starting phonon calculation."
for IR_DIRNAME in $IR_DIRS
do
    $ECHO ": entering ${IR_DIRNAME}"
    cd $IR_DIRNAME
    $ECHO "  starting phonon calculation in ${IR_DIRNAME}..\c"
    ibrun ph.x -in ph.in > ph.out &
    $ECHO " done"
    sleep 2
    cd ..
    $ECHO ""
done
$ECHO "! (`date`) all phonon calculation has been submitted."

wait
$ECHO "! (`date`) all phonon calculation has been finished."

$ECHO "! (`date`) starting matdyn calculation."
for IR_DIRNAME in $IR_DIRS
do
    $ECHO ": entering ${IR_DIRNAME}"
    cd $IR_DIRNAME
    $ECHO "  starting dynmat calculation in ${IR_DIRNAME}..\c"
    dynmat.x -in dynmat.in > dynmat.out
    $ECHO " done"
    sleep 2
    cd ..
    $ECHO ""
done
$ECHO "! (`date`) all dynmat calculation has been finished."

Get IR modes from DYNMAT output

Run this script (get_ir_output.py) and you will get an ir_modes.json.

import glob
import re
import json

def get_mode_irs(filename):
    with open(filename) as fp:
        lines = fp.readlines()
    for idx, line in enumerate(lines):
        if re.search(r"# mode\s+\[cm-1\]\s+\[THz\]\s+IR", line):
            break
    for line in lines[idx+1:]:
        data = [float(d) for d in line.strip().split()]
        if len(data) != 0:
            yield data[1:]
        else:
            return

if __name__ == '__main__':
    ir_modes = []
    for dir_name in glob.glob("ir*/"):
        ir_modes.append({
            "title": dir_name,
            "data": list(get_mode_irs(dir_name + "/dynmat.out"))
        })
    with open("ir_modes.json", "w") as fp:
        json.dump(ir_modes, fp)

Plot Seperate Plots

Plots will be generated in ./plots/ subdirectory.

import matplotlib.pyplot as plt
import json

if __name__ == '__main__':
    with open('ir_modes.json') as fp:
        all_ir_modes = json.load(fp)
    for ir_modes in all_ir_modes:
        title = ir_modes['title']
        modes = ir_modes['data']
        plt.figure(figsize=(5, 2))
        plt.plot(*sum([
            [[freq_cm, freq_cm], [0, ir_intensity]]
            for freq_cm, freq_thz, ir_intensity in modes
        ], []), c='k')
        plt.ylim(ymin=0)
        plt.gca().invert_xaxis()
        plt.xlabel('Frequency / cm$^{-1}$')
        plt.tight_layout()
        plt.savefig('plots/' + title[:-1] + '.png', dpi=300)

Plot Everything Together

import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
import json
import re

if __name__ == '__main__':
    with open('ir_modes.json') as fp:
        all_ir_modes = json.load(fp)
    all_ir_modes = sorted(
        all_ir_modes,
        key=lambda ir_modes: float(re.search(r'(-?[0-9\.]+)', ir_modes['title']).group(1))
    )
    pressures = [float(re.search(r'(-?[0-9\.]+)', ir_modes['title']).group(1)) for ir_modes in all_ir_modes]
    norm = Normalize(vmin=min(pressures), vmax=max(pressures))
    plt.figure(figsize=(6, 2.25))
    legends = []
    for p, ir_modes in zip(pressures, all_ir_modes):
        title = ir_modes['title']
        modes = ir_modes['data']
        lines = plt.plot(*sum([
            [[freq_cm, freq_cm], [0, ir_intensity]]
            for freq_cm, freq_thz, ir_intensity in modes
        ], []), c='k', alpha=norm(p))
        legends.append(lines[0])
    plt.ylim(ymin=0)
    plt.gca().invert_xaxis()
    plt.xlabel('Frequency / cm$^{-1}$')
    plt.legend(legends, [
        '$P$ = %.0f GPa' % p for p in pressures
    ], fontsize='xx-small')
    plt.tight_layout()
    plt.savefig('plots/ir_plot.png', dpi=300)