#!/bin/csh -f

## *** フォルダ構造 ***
## hcpCo/─┬─analysis/
##         ├─in/
##         ├─out/
##         ├─data/
##         ├─template/─hcpCo_Template.in
##         └─qnewt.sh

## *** プロジェクト名 ***
set PROJECT="hcpCo"

## *** c/aの設定 ***
## c/aの下限
set ETA0="1.5"
## c/aの上限
set ETA1="1.7"



## *** 標準入力から値を読み取る ***
## 格子体積
set OMEGA=$1
## c/aの初期値
set ETA_INI=$2

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

## *** 計算するc/aのリストの作成 ***
set ETA_INI=`echo "scale=8; ${ETA_INI}*1/1" | bc -l`
set ETA_INI99=`echo "scale=8; ${ETA_INI}*0.99/1" | bc -l`
set ETA_INI101=`echo "scale=8; ${ETA_INI}*1.01/1" | bc -l`
set ETA_LIST=( ${ETA_INI} ${ETA_INI99} ${ETA_INI101} )

## *** 初期値の計算 ***
set ETA=${ETA_INI}
## ポテンシャルが無ければコピー
if ( ! -e data/${PROJECT}_${OMEGA}_${ETA} ) then
  if ( -e data/${PROJECT}_${ETA} ) then
    cp data/${PROJECT}_${OMEGA} data/${PROJECT}_${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'`
## テンプレートから入力ファイルを作成
if ( ! -e in/${PROJECT}_${OMEGA}_${ETA}.in ) then
  sed 's/'ABOHR'/'${A}'/g' template/${PROJECT}_Template.in | sed 's/'ETA'/'${ETA}'/g' | sed 's/'RMT'/'${RMTA}'/g' | sed 's/'OMEGA'/'${OMEGA}'/g' > in/${PROJECT}_${OMEGA}_${ETA}.in
endif
## 第一原理計算実行
specx < in/${PROJECT}_${OMEGA}_${ETA}.in > out/${PROJECT}_${OMEGA}_${ETA}.out
## ポテンシャルをコピー
cp data/${PROJECT}_${OMEGA}_${ETA} data/${PROJECT}_${OMEGA}
## 結果の読み出し
set ENE_INI=`tail -n 1 data/${PROJECT}_${OMEGA}_${ETA}.info | awk '{print $2}'`
set MAG_INI=`tail -n 1 data/${PROJECT}_${OMEGA}_${ETA}.info | awk '{print $3}'`
echo ${ETA_INI} ${ENE_INI} ${MAG_INI}

## *** 初期値の99%の計算 ***
set ETA=${ETA_INI99}
## ポテンシャルが無ければコピー
if ( ! -e data/${PROJECT}_${OMEGA}_${ETA} ) then
  if ( -e data/${PROJECT}_${OMEGA} ) then
    cp data/${PROJECT}_${OMEGA} data/${PROJECT}_${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'`
## テンプレートから入力ファイルを作成
if ( ! -e in/${PROJECT}_${OMEGA}_${ETA}.in ) then
  sed 's/'ABOHR'/'${A}'/g' template/${PROJECT}_Template.in | sed 's/'ETA'/'${ETA}'/g' | sed 's/'RMT'/'${RMTA}'/g' | sed 's/'OMEGA'/'${OMEGA}'/g' > in/${PROJECT}_${OMEGA}_${ETA}.in
endif
## 第一原理計算実行
specx < in/${PROJECT}_${OMEGA}_${ETA}.in > out/${PROJECT}_${OMEGA}_${ETA}.out
## 結果の読み出し
set ENE_INI99=`tail -n 1 data/${PROJECT}_${OMEGA}_${ETA}.info | awk '{print $2}'`
set MAG_INI99=`tail -n 1 data/${PROJECT}_${OMEGA}_${ETA}.info | awk '{print $3}'`
echo ${ETA_INI99} ${ENE_INI99}  ${MAG_INI99}

## *** 初期値の101%の計算 ***
set ETA=${ETA_INI101}
## ポテンシャルが無ければコピー
if ( ! -e data/${PROJECT}_${OMEGA}_${ETA} ) then
  if ( -e data/${PROJECT}_${OMEGA} ) then
    cp data/${PROJECT}_${OMEGA} data/${PROJECT}_${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'`
## テンプレートから入力ファイルを作成
if ( ! -e in/${PROJECT}_${OMEGA}_${ETA}.in ) then
  sed 's/'ABOHR'/'${A}'/g' template/${PROJECT}_Template.in | sed 's/'ETA'/'${ETA}'/g' | sed 's/'RMT'/'${RMTA}'/g' | sed 's/'OMEGA'/'${OMEGA}'/g' > in/${PROJECT}_${OMEGA}_${ETA}.in
endif
## 第一原理計算実行
specx < in/${PROJECT}_${OMEGA}_${ETA}.in > out/${PROJECT}_${OMEGA}_${ETA}.out
## 結果の読み出し
set ENE_INI101=`tail -n 1 data/${PROJECT}_${OMEGA}_${ETA}.info | awk '{print $2}'`
set MAG_INI101=`tail -n 1 data/${PROJECT}_${OMEGA}_${ETA}.info | awk '{print $3}'`
echo ${ETA_INI101} ${ENE_INI101} ${MAG_INI101}

## *** 推定値の計算 ***
echo "Guess:"
set ENERGY_LIST=( ${ENE_INI} ${ENE_INI99} ${ENE_INI101} )
set ETA_GUESS=`echo "scale=8; (($ETA_LIST[2]*$ETA_LIST[2]-$ETA_LIST[1]*$ETA_LIST[1])*$ENERGY_LIST[3]+($ETA_LIST[1]*$ETA_LIST[1]-$ETA_LIST[3]*$ETA_LIST[3])*$ENERGY_LIST[2]+($ETA_LIST[3]*$ETA_LIST[3]-$ETA_LIST[2]*$ETA_LIST[2])*$ENERGY_LIST[1])/((2*$ETA_LIST[2]-2*$ETA_LIST[1])*$ENERGY_LIST[3]+(2*$ETA_LIST[1]-2*$ETA_LIST[3])*$ENERGY_LIST[2]+(2*$ETA_LIST[3]-2*$ETA_LIST[2])*$ENERGY_LIST[1])" | bc -l`
set ETA=${ETA_GUESS}
## ポテンシャルが無ければコピー
if ( ! -e data/${PROJECT}_${OMEGA}_${ETA} ) then
  if ( -e data/${PROJECT}_${OMEGA} ) then
    cp data/${PROJECT}_${OMEGA} data/${PROJECT}_${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'`
## テンプレートから入力ファイルを作成
if ( ! -e in/${PROJECT}_${OMEGA}_${ETA}.in ) then
  sed 's/'ABOHR'/'${A}'/g' template/${PROJECT}_Template.in | sed 's/'ETA'/'${ETA}'/g' | sed 's/'RMT'/'${RMTA}'/g' | sed 's/'OMEGA'/'${OMEGA}'/g' > in/${PROJECT}_${OMEGA}_${ETA}.in
endif
## 第一原理計算実行
specx < in/${PROJECT}_${OMEGA}_${ETA}.in > out/${PROJECT}_${OMEGA}_${ETA}.out
## ポテンシャルをコピー
cp data/${PROJECT}_${OMEGA}_${ETA} data/${PROJECT}_${OMEGA}
## 結果の読み出し
set ENE_GUESS=`tail -n 1 data/${PROJECT}_${OMEGA}_${ETA}.info | awk '{print $2}'`
set MAG_GUESS=`tail -n 1 data/${PROJECT}_${OMEGA}_${ETA}.info | awk '{print $3}'`
echo ${ETA_GUESS} ${ENE_GUESS} ${MAG_GUESS}
echo ${ETA_GUESS} ${ENE_GUESS} ${MAG_GUESS} >> analysis/${PROJECT}_${OMEGA}.txt
echo "Next:"
echo "./qnewt.sh" ${OMEGA} ${ETA_GUESS}