Created
October 10, 2014 18:28
-
-
Save franciscoadasme/e83f2b3e9959dd6b0f45 to your computer and use it in GitHub Desktop.
Script for Born effective charges calculation using VASP
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Hi,
Is it possible to add a detail step-by-step description on how to use the code? Thank you.