#!/bin/csh -f

## *** フォルダ構造 ***
## bctFe/─┬─analysis/
##         ├─in/
##         ├─out/
##         ├─data/
##         ├─template/─bctFe_Template.in
##         ├─bctFe.sh
##         └─bctFe-Result.sh

setenv OMP_STACKSIZE    200M
limit stacksize unlimited
setenv OMP_NUM_THREADS 4

## *** プロジェクト名 ***
set PROJECT="bctFe"
## ポテンシャルファイル名
set POTENTIAL=${PROJECT}

## *** 格子体積のリスト (Bohr^3) ***
set OMEGA_LIST=( 190 185 180 175 170 165 160 155 150 145 140 135 130 125 120 115 110 105 100 95 90 )
## *** c/aのリスト ***
##set ETA_LIST=( 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 )
set ETA_LIST=( 1.00 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.10 1.10 1.11 1.12 1.13 1.14 1.15 1.16 1.17 1.18 1.19 1.20 1.21 1.22 1.23 1.24 1.25 1.26 1.27 1.28 1.29 1.30 1.31 1.32 1.33 1.34 1.35 1.36 1.37 1.38 1.39 1.40 1.41 1.42 1.43 1.44 1.45 1.46 1.47 1.48 1.49 1.50 )

## *** 第一原理計算 ****
foreach OMEGA ( ${OMEGA_LIST} )
  ## *** MT半径の決定 ***
  ## 左端のc/aでMT半径を計算
  set ETA0=`echo $ETA_LIST[1]`
  set A0=`echo "e((1/3)*l(${OMEGA}/${ETA0}))" | bc -l`
  if (`echo "${ETA0} > sqrt(2)" | bc -l` == 1) then
    set RMTB0=`echo "${A0}/2" | bc -l`
  else
    set RMTB0=`echo "${A0}*sqrt(2+${ETA0}^2)/4" | bc -l`
  endif
  ## 右端のc/aでMT半径を計算
  set ETA1=`echo $ETA_LIST[$#ETA_LIST]`
  set A1=`echo "e((1/3)*l(${OMEGA}/${ETA1}))" | bc -l`
  if (`echo "${ETA1} > sqrt(2)" | bc -l` == 1) then
    set RMTB1=`echo "${A1}/2" | bc -l`
  else
    set RMTB1=`echo "${A1}*sqrt(2+${ETA1}^2)/4" | bc -l`
  endif
  ## 小さい方のMT半径を採用
  if (`echo "${RMTB0} < ${RMTB1}" | bc -l` == 1) then
    set RMTB=`echo $RMTB0`
  else
    set RMTB=`echo $RMTB1`
  endif

  ## *** 第一原理計算 ***
  foreach ETA ( ${ETA_LIST} )
     if ( ! -e data/${POTENTIAL}_${OMEGA}_${ETA} ) then  ## ポテンシャルが無ければコピー
       if ( -e data/${POTENTIAL}_${ETA} ) then
         cp data/${POTENTIAL}_${ETA} data/${POTENTIAL}_${OMEGA}_${ETA}
       endif
     endif
    ## 格子体積とc/aから格子定数aを計算
    set A=`echo "scale=7; e((1/3)*l(2*${OMEGA}/(sqrt(3)*${ETA})))" | bc -l | sed -e 's/^\./0./g'`
    ## MT半径の単位換算
    set RMTA=`echo "scale=5; ${RMTB} / e((1/3)*l(2*${OMEGA}/(sqrt(3)*${ETA})))" | bc -l | sed -e 's/^\./0./g'`
    ## テンプレートから入力ファイルを作成
    sed 's/'ABOHR'/'${A}'/g' template/${PROJECT}_Template.in | sed 's/'ETA'/'${ETA}'/g' | sed 's/'RMTA'/'${RMTA}'/g' | sed 's/'OMEGA'/'${OMEGA}'/g' > in/${PROJECT}_${OMEGA}_${ETA}.in
    ## 第一原理計算実行
    specx < in/${PROJECT}_${OMEGA}_${ETA}.in > out/${PROJECT}_${OMEGA}_${ETA}.out
     ## 次の初期ポテンシャルへコピー
     cp data/${POTENTIAL}_${OMEGA}_${ETA} data/${POTENTIAL}_${ETA}
  end
end