#!/bin/csh
if ( x$1 == "x-h" ) then
  echo Calculates bulk modulus.
  echo Syntax: calcbulk mins maxs ds runcommand
  echo "where mins and maxs are the minimum and maximum strains (e.g. 0.00 0.01)"
  echo "ds is the step size (e.g. 0.005)"
  echo "runcommand is the command to run the ab initio code (e.g. runstruct_vasp)"
  exit
endif
set mins=$1
set maxs=$2
set ds=$3
set strains=`echo $mins $maxs $ds | awk '{for (s=$1; s<=$2+$3/2; s+=$3) {print s}}'`
set runcommand="$argv[4-]"
cp str.out str_relax.out
echo -n > x.tmp
echo -n > y.tmp
foreach strain ( $strains )
  if ( ! -e s_$strain/energy ) then
    mkdir s_$strain
    cellcvrt -sc=`calc 1+$strain` < str_relax.out > s_$strain/str.out
    pushd s_$strain >& /dev/null
    $runcommand
    if ( -e error ) then
	popd >& /dev/null
	touch error
	break
    endif
    popd >& /dev/null
    cellcvrt -sc=`calc "1/(1+$strain)"` < s_$strain/str_relax.out > str_relax.out
  endif
  echo $strain >> x.tmp
  cat s_$strain/energy >> y.tmp
end
if ( ! -e error ) then
  lsfit -x=x.tmp -y=y.tmp -1 -pw=2 > eos.tmp
  echo `cat eos.tmp` `cellcvrt -pv < str.out` |\
    awk '{s=-$2/(2*$3); print s > "strain.out"; print (2./9.)*$3/($4*(1+s)^3) > "bulkmod"; print $1-($2^2)/(4*$3) > "energy"}'
  set ropts=`paste y.tmp x.tmp | awk '{print $1+0,$2}' | sort -n | head -1 | awk '{print $2}'`
  set opts=`cat strain.out`
  cat s_$ropts/str_relax.out | cellcvrt -sc=`calc "(1+$opts)/(1+$ropts)"` > str_relax.out
  cp s_$ropts/dos.out .
endif
