1#!/bin/bash 2# 3# Test the accuracy of the computed stress tensor by comparing it with a 4# finite difference calculation of dE=tr(u*sigma) 5# 6# use: ./check_stress.sh [energy_term sigma_name] file.r 7# If energy_term and sigma_name are not specified, 8# the terms ekin, eps, 9# example: ./check_stress.sh ekin sigma_ekin gs.r 10# 11# The output file file.r must contain the following sequence 12# 1) load a sample from a ground state calculation performed 13# using stress=ON, ref_cell, ecuts 14# 2) the following commands (with arbitrary strain parameters) 15# set wf_dyn LOCKED 16# strain 0.001 0.001 0.001 0 0 0 17# run 0 18# strain -inverse 0.001 0.001 0.001 0 0 0 19# strain -inverse 0.001 0.001 0.001 0 0 0 20# run 0 21 22 23function check_stress_term 24{ 25 # check the accuracy of a single term of the energy and stress tensor 26 eterm=$1 27 etag='<'$eterm'>' 28 sigma_name=$2 29 file=$3 30 31 echo $eterm $sigma_name 32 33 # read volume 34 volume=$(grep '<unit_cell_volume>' $file | awk '{print $2}') 35 echo volume: $volume 36 37 # compute energy difference E(+u)-E(-u) 38 ep=$(grep '<'$eterm'>' $file | awk 'NR==2 {print $2}') 39 em=$(grep '<'$eterm'>' $file | awk 'NR==3 {print $2}') 40 echo ep: $ep 41 echo em: $em 42 de=$(echo "$ep - $em" | bc -l) 43 echo de: $de 44 45 # read strain 46 u=$(grep '<cmd>' $file | grep strain | head -1 | \ 47 sed "s/<cmd>//" | sed "s/<\/cmd>//" | sed "s/\[qbox\]//" | sed "s/strain//") 48 # echo $u 49 strain=($u) 50 uxx=$(echo ${strain[0]}) 51 uyy=$(echo ${strain[1]}) 52 uzz=$(echo ${strain[2]}) 53 uxy=$(echo ${strain[3]}) 54 uyz=$(echo ${strain[4]}) 55 uxz=$(echo ${strain[5]}) 56 echo strain: $uxx $uyy $uzz $uxy $uyz $uxz 57 58 # read computed stress tensor 59 sxx=$(grep '<'${sigma_name}_xx'>' $file |head -1| awk '{print $2}') 60 syy=$(grep '<'${sigma_name}_yy'>' $file |head -1| awk '{print $2}') 61 szz=$(grep '<'${sigma_name}_zz'>' $file |head -1| awk '{print $2}') 62 sxy=$(grep '<'${sigma_name}_xy'>' $file |head -1| awk '{print $2}') 63 syz=$(grep '<'${sigma_name}_yz'>' $file |head -1| awk '{print $2}') 64 sxz=$(grep '<'${sigma_name}_xz'>' $file |head -1| awk '{print $2}') 65 echo stress: $sxx $syy $szz $sxy $syz $sxz 66 67 # compute energy change from 2 * volume * tr(sigma*u) 68 de_comp=$(echo "-2.0 * $volume * ($uxx * $sxx + $uyy * $syy + $uzz * $szz + \ 69 2.0 * ( $uxy * $sxy + $uyz * $syz + $uxz * $sxz) )" | bc -l) 70 printf "E(+u)-E(-u): %15.6e\n" $de 71 printf "2*v*tr(sigma*u): %15.6e\n" $de_comp 72 73 # compute absolute error 74 abs_err=$(echo "$de - $de_comp" | bc -l) 75 printf "abs error: %15.6e\n" $abs_err 76 rel_err=$(echo "$abs_err / $de" | bc -l) 77 printf "rel error: %15.6e\n" $rel_err 78 79 echo "===============================================================" 80} 81 82if (( $# == 3 )) 83then 84 # check_stress.sh term sigma_name file.r 85 echo "===============================================================" 86 echo " check_stress.sh" ${*} 87 echo "===============================================================" 88 check_stress_term $1 $2 $3 89 exit 90fi 91 92if (( $# == 1 )) 93then 94 # check_stress.sh file.r 95 echo "===============================================================" 96 echo " check_stress.sh" ${*} 97 echo "===============================================================" 98 check_stress_term "ekin" "sigma_ekin" $1 99 check_stress_term "eps" "sigma_eps" $1 100 check_stress_term "econf" "sigma_econf" $1 101 check_stress_term "enl" "sigma_enl" $1 102 check_stress_term "exc" "sigma_exc" $1 103 check_stress_term "esr" "sigma_esr" $1 104 check_stress_term "etotal" "sigma_eks" $1 105 exit 106fi 107 108# incorrect number of arguments. Print usage. 109echo " use: $0 [energy_term sigma_name] file.r" 110