Script for Born effective charges calculation using VASP

gistfile1.sh

#!/usr/bin/env bash

# Script initialization
script=$(basename "$0")
desc='Perform Born effective charges calculation, which are useful for
analyzing spontaneous polarization, using the Vienna Ab initio simulation
package, VASP.'

declare -iA errno
errno[EUNK]=1
errno[EARG]=2
errno[ENOVASP]=4
errno[ECHGFAILED]=8
errno[EBERRYFAILED]=16

# Helpers
usage() { echo "usage: $script [-h] [-x <float>] [-y <float>] [-z <float>]\
 [-n <int>] [-c \"<float> <float> <float>\"] idx [idx ...]" >&2; }
usage_and_exit() { usage; exit "${2:-1}"; }
help() {
    usage;
    echo -e "\n$desc\n"

    echo 'positional arguments:'
    echo ' idx            zero-based indices of one or more atoms to move.'

    echo ''

    echo 'optional arguments:'
    echo ' -x DELTAX      ion displacement delta over x in angstroms.'
    echo ' -y DELTAY      ion displacement delta over y in angstroms.'
    echo ' -z DELTAZ      ion displacement delta over z in angstroms.'
    echo ' -n STEPS       number of displacements. Default is 2. Note: a total'
    echo '                of n+1 polarization calculations are done including'
    echo '                the undistorted system.'
    echo ' -k KPOINTS     number of points on the strings in the IGPAR'
    echo '                direction. See NPPSTR flag in the VASP manual.'
    echo '                Default is 20.'
    echo ' -d DIPOL       center of cell (fractional coordinates). Specifies'
    echo '                the origin with respect to which the ionic'
    echo '                contribution to the dipole moment in the cell is'
    echo '                calculated. See DIPOL VASP keyword for more info.'
    echo '                Default is "0.5 0.5 0.5".'
    echo ' -c CMD         Command to execute vasp. Useful when using a resource'
    echo '                manager like SLURM (e.g. srun -o OUTPUT -e ERR vasp)'
    echo '                or when using custom vasp executable path. Default '
    echo '                is "vasp".'
    echo ' -o OUTPUT      Redirect stdout to the named file.'
    echo ' -e ERR         Redirect stderr to the named file.'
    echo ' -h             show this help message and exit.'

    echo ''

    echo 'return codes'
    echo ' 1  EUNK        Unknown error.'
    echo ' 2  EARG        Bad argument format/value.'
    echo ' 4  ENOVASP     No VASP executable was found.'
    echo ' 8  ECHGFAILED  CHGCAR SP calculation failed.'
    echo '16  EBERRYFAILED  Electronic polarization calculation failed.'
    exit 1
}
warn() { echo "WARNING! $1."; }
err() { echo "ERROR! $1." >&2; }
err_and_exit() { err "$1"; exit "${2:-1}"; }

function check_type {
    [[ $4 =~ $2 ]] && retno=0 || retno=1
    if [[ $retno -ne 0 ]]; then
        err_and_exit "Expected $1 value for -$3 option, got: $4" ${errno[EARG]}
    fi
    return $retno
}
check_int() { check_type 'integer' ^[+-]?[0-9]+$ "$1" "$2"; }
check_float() { check_type 'float' ^[+-]?[0-9]+\.?[0-9]*$ "$1" "$2"; }

sum() { awk 'BEGIN {t=0; for (i in ARGV) t+=ARGV[i]; print t}' "$@"; }
sub() { awk 'BEGIN {t=ARGV[1]*2; for (i in ARGV) t-=ARGV[i]; print t}' "$@"; }
mean() { sum=$(sum "$@"); awk -v sum="$sum" "BEGIN {print sum/(ARGC-1)}" "$@"; }
transpose() {
    local args=("$@")
    local -a vx; local -a vy; local -a vz
    for i in $(seq 0 3 $((${#args[@]}-1))); do
        vx+=(${args[$i]})
        vy+=(${args[$((i+1))]})
        vz+=(${args[$((i+2))]})
    done
    echo ${vx[*]} ${vy[*]} ${vz[*]}
}
apply_v() {
    local action=($1); shift
    local arr_t=($(transpose "$@"))
    local sum_arr=()
    local len=${#arr_t[@]}
    local step=$(bc <<< "$len/3")
    for i in $(seq 0 "$step" $((len-1))); do
        cmd=("$action" "${arr_t[@]:$i:$step}")
        sum_arr+=($("${cmd[@]}"))
    done
    echo "${sum_arr[@]}"
}
sum_v() { apply_v 'sum' "$@"; }
sub_v() { apply_v 'sub' "$@"; }
mean_v() { apply_v 'mean' "$@"; }

exec_vasp() { cmd=(${opts[cmd]}); "${cmd[@]}"; }

# Parse arguments
declare -A opts
opts[dx]=0
opts[dy]=0
opts[dz]=0
opts[steps]=2
opts[kpts]=20
opts[dcenter]='.5 .5 .5'
opts[cmd]='vasp'

[[ $# -eq 0 ]] && usage_and_exit
while getopts ":x:y:z:n:k:d:c:o:e:h" flag; do
    case $flag in
        x) check_float "$flag" "$OPTARG" && opts[dx]=$OPTARG;;
        y) check_float "$flag" "$OPTARG" && opts[dy]=$OPTARG;;
        z) check_float "$flag" "$OPTARG" && opts[dz]=$OPTARG;;
        n) check_int "$flag" "$OPTARG" && opts[steps]=$OPTARG;;
        k) check_int "$flag" "$OPTARG" && opts[kpts]=$OPTARG;;
        d)
            for i in $OPTARG; do
                check_float "$flag" "$i"
                if [[ $(bc <<< "$i < 0 || $i > 1") -eq 1 ]]; then
                    err "Center coordinate must be fractional, got: $i"
                    exit ${errno[EARG]}
                fi
            done
            opts[dcenter]=$OPTARG
            ;;
        c) opts[cmd]=$OPTARG;;
        o) opts[outfile]=$OPTARG;;
        e) opts[errfile]=$OPTARG;;
        h) help;;
        \?)
            err "Invalid option: -$OPTARG"
            usage_and_exit
            ;;
        :)
            err_and_exit "Option -$OPTARG requires an argument";;
    esac
done
shift $((OPTIND-1))
idxs=("$@")

# Check number of remaining (non-option) arguments
if [[ ${#idxs[@]} -eq 0 ]]; then
    err 'No atom indices were given' ${errno[EARG]}
    usage_and_exit
fi

# Check whether one or more delta are non-zero
if [[ $(bc <<< "${opts[dx]} == 0 && ${opts[dy]} == 0 \
        && ${opts[dz]} == 0") -eq 1 ]]; then
    err_and_exit 'One or more delta must be non-zero' ${errno[EARG]}
fi

# Set output and err files
[[ -n "${opts[outfile]}" ]] && exec 1> "${opts[outfile]}"
[[ -n "${opts[errfile]}" ]] && exec 2> "${opts[errfile]}"

# Check whether vasp executable is present
if [[ "${opts[cmd]}" == "vasp" ]]; then
    which vasp &> /dev/null
    [[ $? -ne 0 ]] && err_and_exit 'VASP executable not found' ${errno[ENOVASP]}
fi

# Check that INCAR file represents a SP calculation
if grep -Eq '^ *[^#]? *IBRION *= *[0-9]' INCAR; then
    err 'Calculation appears to be an optimization. Please check the INCAR' \
        'file for IBRION keyword'
    exit ${errno[EARG]}
fi

# Workspace setup
rootdir=$(pwd)
workdir=$rootdir/${script}_workdir
mkdir -p "$workdir"
# outdir=$rootdir/Output
# mkdir -p "$outdir"

start_time=$(date +%s)

# Check whether CHGCAR exists; run an SP otherwise
if [[ -f CHGCAR ]]; then
    echo 'Using existing CHGCAR file.'
else
    warn 'CHGCAR file was not found'
    printf 'Running SP to calculate CHGCAR...'

    chgdir=$workdir/chgdir
    mkdir -p "$chgdir"
    cp KPOINTS POSCAR POTCAR "$chgdir"
    # ensure to print out the CHGCAR file
    sed -r 's/(LCHARG *= *)[\.A-Z]+(.+)/\1.TRUE.\2/' INCAR > "$chgdir/INCAR"

    cd "$chgdir"
    exec_vasp

    if [[ $? -ne 0 || -s ERR || ! -f CHGCAR ]]; then
        echo -e \\n
        err 'CHGCAR calculation did fail'
        echo 'Here are the last 20 lines from the stderr buffer:' >&2
        tail -20 ERR >&2
        rm -r "$chgdir"
        exit ${errno[ECHGFAILED]}
    fi

    cp CHGCAR "$rootdir"
    rm -r "$chgdir"
    echo ' done.'
fi

# Print out parameters
echo -e \\n'Parameters'
printf ' delta=( %.3f %.3f %.3f ) Angst\n' "${opts[dx]}" "${opts[dy]}" "${opts[dz]}"
printf ' k-points=%i\n' "${opts[kpts]}"
printf ' dipole center=( %.2f %.2f %.2f )\n' ${opts[dcenter]}
printf ' cmd=%s\n\n' "${opts[cmd]}"

# Declare global array holders
declare -a er_arr # holds the total polarization for each system

for step in $(seq -w 0 "${opts[steps]}"); do
    [[ 10#$step -eq 0 ]] && job=undistorted || job="distorted step $step"
    echo "> Evaluating system $job..."

    [[ 10#$step -eq 0 ]] && tmpdir=undistorted || tmpdir=distorted_$step
    tmpdir=$workdir/$tmpdir

    unset -v er_ev_arr # holds electronic polarization term in each axis
    unset -v er_bp_arr # holds electronic polarization term in each axis
    unset -v ion_p # holds ionic dipole moment in each axis

    for axis in $(seq 3); do
        cd "$rootdir"

        tmpdir2=${tmpdir}_$axis
        mkdir -p "$tmpdir2"

        jobname=$(basename "$tmpdir2")
        echo "Setting up $jobname..."
        cp KPOINTS POSCAR POTCAR CHGCAR "$tmpdir2"

        # Add displacements if needed
        if [[ 10#$step -ne 0 ]]; then
            dx=$(bc -l <<< "${opts[dx]}*$step")
            dy=$(bc -l <<< "${opts[dy]}*$step")
            dz=$(bc -l <<< "${opts[dz]}*$step")

            start=$(awk '/(Direct|Cartesian)/{ print NR; exit }' POSCAR)
            for idx in $idxs; do
                lnum=$((start+idx+1))
                line=$(sed "${lnum}q;d" POSCAR)
                repl=$(echo "$line" | awk -v dx="$dx" -v dy="$dy" -v dz="$dz" \
                       '{printf "%.4f %.4f %.4f", $1+dx, $2+dy, $3+dz}')
                sed -i "${lnum}s/.*/$repl/" "$tmpdir2/POSCAR"
            done
        fi

        # Modify INCAR file
        regex='(ISTART|ICHARG|INIWAV|LCHARG|LWAVE|^ *#.+)'
        sed -r "/$regex/d" INCAR | sed '/^$/d' > "$tmpdir2/INCAR"
        echo "LWAVE = .FALSE." >> "$tmpdir2/INCAR"
        echo "LCHARG = .FALSE." >> "$tmpdir2/INCAR"
        echo 'LBERRY = .TRUE.' >> "$tmpdir2/INCAR"
        echo "NPPSTR = ${opts[kpts]}" >> "$tmpdir2/INCAR"
        echo "DIPOL = ${opts[dcenter]}" >> "$tmpdir2/INCAR"
        echo "IGPAR = ${axis}" >> "$tmpdir2/INCAR"

        cd "$tmpdir2"
        printf 'Executing electronic polarization job...'
        exec_vasp

        # Verify that everything went fine
        retno=$?
        if [[ $retno -ne 0 ]]; then
            err "Something went wrong with job $jobname, exit code: $retno"
            echo 'Here are the last 20 lines from the stderr buffer:' >&2
            tail -20 ERR >&2
            rm -r "$tmpdir2"
            exit ${errno[EBERRYFAILED]}
        fi

        # Print some status info
        echo ' done.'
        grep 'Total CPU' OUTCAR | \
         awk '{printf "Job took about %.2f minutes of CPU time.\n\n", $6/60}'

        # Get values
        er_ev_arr+=($(awk '/e<r>_ev/{print $2, $3, $4}' OUTCAR))
        er_bp_arr+=($(awk '/e<r>_bp/{print $2, $3, $4}' OUTCAR))
        ion_p=($(awk '/ionic dipole/{print $5, $6, $7}' OUTCAR))

        #rm -r "$tmpdir2"
    done

    # Calculate total polarization for current system
    er_ev=($(mean_v "${er_ev_arr[@]}"))
    er_bp=($(sum_v "${er_bp_arr[@]}"))
    er_el=($(sum_v "${er_ev[@]}" "${er_bp[@]}"))
    er=($(sum_v "${er_el[@]}" "${ion_p[@]}"))
    er_arr+=("${er[@]}")

    echo "Total polarization for $job system:"
    printf " e<r>_ev=( %10.5f %10.5f %10.5f ) e*Angst\n" "${er_ev[@]}"
    printf " e<r>_bp=( %10.5f %10.5f %10.5f ) e*Angst\n" "${er_bp[@]}"
    printf " e<r>_el=( %10.5f %10.5f %10.5f ) e*Angst\n" "${er_el[@]}"
    printf "  p[ion]=( %10.5f %10.5f %10.5f ) e*Angst\n" "${ion_p[@]}"
    printf "    e<r>=( %10.5f %10.5f %10.5f ) e*Angst\n" "${er[@]}"
    echo -e ''\\n
done

echo -e "--------------------------------------------------"\\n
echo -e "All electronic polarization calculations are done."\\n
echo "Total change in polarization (d_er = er_dist - er_undist):"

er_undist=("${er_arr[@]:0:3}")
declare -a d_er_arr
for i in $(seq 3 3 $((${#er_arr[@]}-1))); do
    er_dist=("${er_arr[@]:$i:3}")
    d_er=($(sub_v "${er_dist[@]}" "${er_undist[@]}"))
    d_er_arr+=("${d_er[@]}")
    step=$((i/3))
    dx=$(bc -l <<< "${opts[dx]}*$step")
    dy=$(bc -l <<< "${opts[dy]}*$step")
    dz=$(bc -l <<< "${opts[dz]}*$step")

    printf ' d=( %6.3f %6.3f %6.3f ) d_e<r>=( %10.5f %10.5f %10.5f ) e*Angst\n' \
     "$dx" "$dy" "$dz" "${d_er[@]}"
done
echo ''

echo "Born effective charges (Z* = d_er / delta):"
for i in $(seq 0 3 $((${#d_er_arr[@]}-1))); do
    d_er=("${d_er_arr[@]:$i:3}")
    step=$((i/3+1))
    dx=$(bc -l <<< "${opts[dx]}*$step")
    dy=$(bc -l <<< "${opts[dy]}*$step")
    dz=$(bc -l <<< "${opts[dz]}*$step")

    [[ $(bc -l <<< "$dx != 0") -eq 1 ]] && zx=$(bc -l <<< "${d_er[0]}/$dx") || zx=0
    [[ $(bc -l <<< "$dy != 0") -eq 1 ]] && zy=$(bc -l <<< "${d_er[1]}/$dy") || zy=0
    [[ $(bc -l <<< "$dz != 0") -eq 1 ]] && zz=$(bc -l <<< "${d_er[2]}/$dz") || zz=0

    printf ' d=( %6.3f %6.3f %6.3f ) Z*=( %10.5f %10.5f %10.5f ) e*Angst\n' \
     "$dx" "$dy" "$dz" "$zx" "$zy" "$zz"
done
echo ''

end_time=$(date +%s)
printf "Born effective charges calculation took a total of %.2f minutes.\n" \
 "$(bc -l <<< "($end_time-$start_time)/60")"

cd "$rootdir"
[[ ! -s ERR ]] && rm ERR

#rm -r "$workdir"
exit 0

Hi,
Is it possible to add a detail step-by-step description on how to use the code? Thank you.