Add test of stress tensor

test/stress/Makefile

+include ../Makefile.inc
+all: test.r
+clean:
+	rm -f gs.r gs.xml gs_diag.r gs_offdiag.r test.r
+
+gs.r: params.i
+gs_diag.r: params.i
+gs_offdiag.r: params.i
+
+gs.r: gs.i
+	$(QBOX_CMD) gs.i > gs.r
+
+gs_diag.r: gs_diag.i gs.r
+	$(QBOX_CMD) gs_diag.i > gs_diag.r
+
+gs_offdiag.r: gs_offdiag.i gs.r
+	$(QBOX_CMD) gs_offdiag.i > gs_offdiag.r
+
+test.r: gs_diag.r gs_offdiag.r
+	./check_stress.sh gs_diag.r > test.r
+	./check_stress.sh gs_offdiag.r >> test.r

test/stress/README

+--------------------------------------------------------------------------------
+ Test of the accuracy of the stress tensor
+--------------------------------------------------------------------------------
+Notes:
+- To change the choice of xc, ecut, ecuts, edit the file params.i 
+- In order to test off-diagonal elements of the stress tensor,
+  the starting unit cell in gs.i is distorted arbitrarily.
+--------------------------------------------------------------------------------

test/stress/Si_VBC_LDA-1.0.xml

+../potentials/Si_VBC_LDA-1.0.xml
\ No newline at end of file

test/stress/check_stress.sh

+#!/bin/bash
+#
+# Test the accuracy of the computed stress tensor by comparing it with a 
+# finite difference calculation of dE=tr(u*sigma)
+#
+# use: ./check_stress.sh [energy_term sigma_name] file.r
+#      If energy_term and sigma_name are not specified, 
+#      the terms ekin, eps, 
+# example: ./check_stress.sh ekin sigma_ekin gs.r
+#
+# The output file file.r must contain the following sequence
+# 1) load a sample from a ground state calculation performed
+#    using stress=ON, ref_cell, ecuts
+# 2) the following commands (with arbitrary strain parameters)
+#    set wf_dyn LOCKED
+#    strain 0.001 0.001 0.001 0 0 0
+#    run 0
+#    strain -inverse 0.001 0.001 0.001 0 0 0
+#    strain -inverse 0.001 0.001 0.001 0 0 0
+#    run 0
+
+
+function check_stress_term
+{
+  # check the accuracy of a single term of the energy and stress tensor
+  eterm=$1
+  etag='<'$eterm'>'
+  sigma_name=$2
+  file=$3
+
+  echo $eterm $sigma_name
+  
+  # read volume
+  volume=$(grep '<unit_cell_volume>' $file | awk '{print $2}')
+  echo volume: $volume
+  
+  # compute energy difference E(+u)-E(-u)
+  ep=$(grep '<'$eterm'>' $file | awk 'NR==2 {print $2}')
+  em=$(grep '<'$eterm'>' $file | awk 'NR==3 {print $2}')
+  echo ep: $ep
+  echo em: $em
+  de=$(echo "$ep - $em" | bc -l)
+  echo de: $de
+  
+  # read strain 
+  u=$(grep '<cmd>' $file | grep strain | head -1 | \
+    sed "s/<cmd>//" | sed "s/<\/cmd>//" | sed "s/\[qbox\]//" | sed "s/strain//")
+  # echo $u
+  strain=($u)
+  uxx=$(echo ${strain[0]})
+  uyy=$(echo ${strain[1]})
+  uzz=$(echo ${strain[2]})
+  uxy=$(echo ${strain[3]})
+  uyz=$(echo ${strain[4]})
+  uxz=$(echo ${strain[5]})
+  echo strain: $uxx $uyy $uzz $uxy $uyz $uxz
+  
+  # read computed stress tensor
+  sxx=$(grep '<'${sigma_name}_xx'>' $file |head -1| awk '{print $2}')
+  syy=$(grep '<'${sigma_name}_yy'>' $file |head -1| awk '{print $2}')
+  szz=$(grep '<'${sigma_name}_zz'>' $file |head -1| awk '{print $2}')
+  sxy=$(grep '<'${sigma_name}_xy'>' $file |head -1| awk '{print $2}')
+  syz=$(grep '<'${sigma_name}_yz'>' $file |head -1| awk '{print $2}')
+  sxz=$(grep '<'${sigma_name}_xz'>' $file |head -1| awk '{print $2}')
+  echo stress: $sxx $syy $szz $sxy $syz $sxz
+  
+  # compute energy change from 2 * volume * tr(sigma*u)
+  de_comp=$(echo "-2.0 * $volume * ($uxx * $sxx + $uyy * $syy + $uzz * $szz + \
+   2.0 * ( $uxy * $sxy + $uyz * $syz + $uxz * $sxz) )" | bc -l)
+  printf "E(+u)-E(-u):     %15.6e\n" $de
+  printf "2*v*tr(sigma*u): %15.6e\n" $de_comp
+  
+  # compute absolute error
+  abs_err=$(echo "$de - $de_comp" | bc -l)
+  printf "abs error:       %15.6e\n" $abs_err
+  rel_err=$(echo "$abs_err / $de" | bc -l)
+  printf "rel error:       %15.6e\n" $rel_err
+  
+  echo "==============================================================="
+}
+
+if (( $# == 3 ))
+then
+  # check_stress.sh term sigma_name file.r
+  echo "==============================================================="
+  echo " check_stress.sh" ${*}
+  echo "==============================================================="
+  check_stress_term $1 $2 $3
+  exit
+fi
+
+if (( $# == 1 ))
+then
+  # check_stress.sh file.r
+  echo "==============================================================="
+  echo " check_stress.sh" ${*}
+  echo "==============================================================="
+  check_stress_term "ekin"    "sigma_ekin" $1
+  check_stress_term "eps"     "sigma_eps"  $1
+  check_stress_term "econf"   "sigma_econf"  $1
+  check_stress_term "enl"     "sigma_enl"  $1
+  check_stress_term "exc"     "sigma_exc"  $1
+  check_stress_term "esr"     "sigma_esr"  $1
+  check_stress_term "etotal"  "sigma_eks"  $1
+  exit
+fi
+
+# incorrect number of arguments. Print usage.
+echo " use: $0 [energy_term sigma_name] file.r"

test/stress/gs.i

+# Si4 ground state
+set cell 16 0 0  0 16 0  0 0 16
+set ref_cell 18 0 0  0 18 0  0 0 18
+species silicon Si_VBC_LDA-1.0.xml
+atom Si1 silicon  3.700 -0.100  0.300
+atom Si2 silicon -0.100  2.800 -0.200
+atom Si3 silicon -3.700  0.100  0.300
+atom Si4 silicon -0.100 -2.800 -0.200
+strain 0.02 0.04 0.06 0.03 0.05 0.07
+params.i
+set wf_dyn PSDA
+set stress ON
+set debug STRESS
+randomize_wf
+set scf_tol 1.e-9
+run 0 300
+save gs.xml

test/stress/gs_diag.i

+load gs.xml
+params.i
+set wf_dyn PSDA
+set stress ON
+set debug STRESS
+run 0
+set wf_dyn LOCKED
+strain 0.0001 0.0002 0.0003  0 0 0
+run 0
+strain -inverse 0.0001 0.0002 0.0003  0 0 0
+strain -inverse 0.0001 0.0002 0.0003  0 0 0
+run 0

test/stress/gs_offdiag.i

+load gs.xml
+params.i
+set wf_dyn PSDA
+set stress ON
+set debug STRESS
+run 0
+set wf_dyn LOCKED
+strain 0 0 0  0.001 0.002 -0.001
+run 0
+strain -inverse 0 0 0  0.001 0.002 -0.001
+strain -inverse 0 0 0  0.001 0.002 -0.001
+run 0

test/stress/params.i

+set ecut 12
+set ecuts 10
+set xc LDA

test/stress/sample.xsd

+../../xml/sample.xsd
\ No newline at end of file

test/stress/species.xsd

+../../xml/species.xsd
\ No newline at end of file