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