* シリコン結晶 [#v77e93e6]
この例ではダイアモンド構造のシリコンを取り上げ、セル（格子定数）の最適化および全エネルギーと格子定数の平面波カットオフ依存性を調べる方法を記述する。擬ポテンシャルにはSi_pbe (ノルム保存擬型) を用いる。
** セルの最適化 [#e566418c]
現在のSTATEにはストレスの計算は実装されておらず、単位胞の形状および体積を自動で最適化することはできない。そのため全エネルギーを格子定数の関数として計算し、フィッティングにより平衡格子定数および体積を決定する必要がある。以下にその手順を記述する。~
ここでは以下のスクリプトを使用し、自動的に入力ファイルを生成し、STATEを実行する。
#!/bin/sh
#
# Set the number of CPUs
#
NCPU=4
#
# Set the source directory of STATE
#
ROOTDIR=${HOME}/STATE/src/state/
SRCDIR=${ROOTDIR}/src
PPDIR=${HOME}/STATE/gncpp
MPIRUN="mpirun -np ${NCPU}"
STATE=STATE
#
# Prepare the STATE executable
#
ln -fs ${SRCDIR}/${STATE}
#
# Prepare the pseudopotential data
#
ln -fs ${PPDIR}/pot.Si_pbe1 fort.37
#
# Set the minimum and maximum lattice parameter and mesh
#
AMIN=10.25
AMAX=10.45
MA=20
#
for (( i=0; i <= ${MA}; i++ ))
do
#
# Set the lattice parameters
#
A=`echo "scale=2; ${AMIN} + ( ${AMAX} - ${AMIN} ) * ${i} / ${MA}" | bc`
B=${A}
C=${A}
#
# Input/output files
#
INPUT_FILE=nfinp_a${A}
OUTPUT_FILE=nfout_a${A}
cat > ${INPUT_FILE} << EOF
0 0 0 0 0 0 : INPUT_CTRL(1:6) (DUMMY)
4.00 8.00 1 2 2 : GMAX GMAXP NTYP NATM NATM2
227 2 : num_space_group, type
${A} ${B} ${C} 90.00 90.00 90.00 : A B C ALPHA BETA GAMMA
4 4 4 2 2 2 : N1 N2 N3 M1 M2 M3
0 0 : NCORD, NINV
0.00 0.00 0.00 1 1 1 : CPS(1,1:3) IWEI IMDTYP ITYP
0.25 0.25 0.25 1 1 1 : CPS(2,1:3) IWEI IMDTYP ITYP
14 0.5000 28.09 6 1 0.0 : ATOMN ALFA AMION ILOC IVAN ZETA1
0 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACC
0 1 : IPRE IPRI
200 200 0 84200.00 0 : NMD1 NMD2 LAST_ITER CPUMAX IFSTOP
6 1 : WAY_MIX MIX_WHAT
0 20 0.60 : ITER_START KBXMIX MIX_ALPHA
0.20 0.30 0.20 0.20 0.20 : DTIM1 DTIM2 DTIM3 DTIM4 DTIM
300.00 4 1 0.50D-09 : DTIO IMDALG IEXPL EDELTA
-0.0020 0.50D+03 0 : WIDTH FORCCR ISTRESS
ggapbe 1 : xctype,kspin
2.00 : DETSTM
102 : NBZTYP
0 0 0 : NKX NKY NKZ (DUMMY)
0 0 0 : NKX2 NKY2 NKZ2 (DUMMY)
8 : KEG
1 : NEXTST
0 : (DUMMY)
2 : IMSD
0 : EVALUATE_EKO_DIFF
0 : NPDOSAO
0 0.00 : SM_N DOPING
EOF
#
# Run!
#
${MPIRUN} ${STATE} < $INPUT_FILE > $OUTPUT_FILE
#
done
このスクリプトを環境に合わせて修正し実行すると、nfinp_a*、nfout_a*が得られる。~
次に
for f in nfout_a*
do a=${f##nfout_a}
etot=`grep 'TOTAL ENERGY' $f | grep A.U. | awk '{print $4}'`
echo $a $etot
done
を実行すると格子定数の関数として全エネルギーが表示される。~
ここでは添付のstate2ev.shを用いて全エネルギーを体積の関数としてプロットしよう。以下を実行する。
state2ev.sh nfout_a* > etot.dat
次に添付のeosfitを使用し、全エネルギーをMurnaghanの状態方程式でフィットする。
eosfit etot.dat
この例での格子定数の間隔はとても小さいので全エネルギーは滑らかな曲線になっていない（滑らかな曲線を得るためには10x10x10程度のメッシュが必要）。~
eosfitを実行するとeosfit.paramというファイルが生成される。
そこで示されたb0, b0', v0, e0はそれぞれ体積弾性率、体積弾性率の微分、平衡体積、平衡位置での全エネルギーを示す。今の例では平衡体積として41.0191148309001 Angstrom^3が得られ、平衡格子定数は5.4746 Angstromとなる（実験値は5.416 Angstrom）。
** 全エネルギーの平面波カットオフエネルギー依存生 [#x74524ae]
以下のようなスクリプトを実行し格子定数一定で全エネルギーの平面波カットオフ依存生を調べる。
#!/bin/sh
#
# Set the number of CPUs
#
NCPU=4
#
# Set the source directory of STATE
#
ROOTDIR=${HOME}/STATE/src/state/
SRCDIR=${ROOTDIR}/src
PPDIR=${HOME}/STATE/gncpp
MPIRUN="mpirun -np ${NCPU}"
STATE=STATE
#
# Prepare the STATE executable
#
ln -fs ${SRCDIR}/${STATE}
#
# Prepare the pseudopotential data
#
ln -fs ${PPDIR}/pot.Si_pbe1 fort.37
#
# Set the minimum and maximum wave vector cutoffs
#
GMAX_MIN=2.0
GMAX_MAX=10.0
MG=8
#
for (( i=0; i <= ${MG}; i++ ))
do
#
# Set the cutoff energies
GMAX=`echo "scale=2; ${GMAX_MIN} + ( ${GMAX_MAX} - ${GMAX_MIN} ) * ${i} / ${MG}" | bc`
GMAXP=`echo "scale=2; ${GMAX} * 2.0" | bc`
#
# Input/output files
#
INPUT_FILE=nfinp_gmax${GMAX}
OUTPUT_FILE=nfout_gmax${GMAX}
cat > ${INPUT_FILE} << EOF
0 0 0 0 0 0 : INPUT_CTRL(1:6) (DUMMY)
${GMAX} ${GMAXP} 1 2 2 : GMAX GMAXP NTYP NATM NATM2
227 2 : num_space_group, type
10.2347 10.2347 10.2347 90.0 90.0 90.0 : A B C ALPHA BETA GAMMA
4 4 4 2 2 2 : N1 N2 N3 M1 M2 M3
0 0 : NCORD, NINV
0.00 0.00 0.00 1 1 1 : CPS(1,1:3) IWEI IMDTYP ITYP
0.25 0.25 0.25 1 1 1 : CPS(2,1:3) IWEI IMDTYP ITYP
14 0.5000 28.09 6 1 0.0 : ATOMN ALFA AMION ILOC IVAN ZETA1
0 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACC
0 1 : IPRE IPRI
200 200 0 84200.00 0 : NMD1 NMD2 LAST_ITER CPUMAX IFSTOP
6 1 : WAY_MIX MIX_WHAT
0 20 0.60 : ITER_START KBXMIX MIX_ALPHA
0.20 0.30 0.20 0.20 0.20 : DTIM1 DTIM2 DTIM3 DTIM4 DTIM
300.00 4 1 0.50D-09 : DTIO IMDALG IEXPL EDELTA
-0.0020 0.50D+03 0 : WIDTH FORCCR ISTRESS
ggapbe 1 : xctype,kspin
2.00 : DETSTM
101 : NBZTYP
0 0 0 : NKX NKY NKZ (DUMMY)
0 0 0 : NKX2 NKY2 NKZ2 (DUMMY)
8 : KEG
1 : NEXTST
0 : (DUMMY)
2 : IMSD
0 : EVALUATE_EKO_DIFF
0 : NPDOSAO
0 0.00 : SM_N DOPING
EOF
#
# Run!
#
${MPIRUN} ${STATE} < $INPUT_FILE > $OUTPUT_FILE
#
done
#
###############
# END OF FILE #
###############
(gmax=8,9は恐らく実空間のグリッドの問題で計算が行えなかった)~
得られた全エネルギーとカットオフの関係は以下のようになる。
#GMAX ETOT (Hartree)
2.00 -7.68389161
3.00 -7.83436760
4.00 -7.87325631
5.00 -7.87585304
6.00 -7.87874066
7.00 -7.87922178
8.00 -7.87928813
9.00 -7.87935331
10.00 -7.87935680
カットオフエネルギーを増加していくと変分原理によりエネルギーが減少していくことが分かる。
![[PukiWiki]](image/logo_state-senri3_80x80.png "[PukiWiki]"){kind=logo}
