![[PukiWiki]](image/logo_state-senri3_80x80.png "[PukiWiki]"){kind=logo}
シリコン結晶 (セル最適化) †ここでは結晶構造の最適化例としてシリコン結晶の平衡格子定数を求めましょう。 電子状態計算における最適化とは一般的には全エネルギーを最小化する電子の波動関数や構造を求めることを意味しており、エネルギーを最小化する格子定数を平衡格子定数と呼んでいます。 現在のワーキングディレクトリが~/STATE/examplesであるとして、Siディレクトリに移動します。 $ cd Si ここでlsを実行し、nfで始まる入力ファイル群（nf*）や.shで終わるジョブスクリプト（*.sh）が含まれていることを確認します。 以下を実行して入力ファイルnfinp_scfを再確認しましょう。 $ cat nfinp_scf #
# Crystalline silicon in the diamond structure
#
WF_OPT DAV
NTYP 1
NATM 2
TYPE 2
NSPG 227
GMAX 4.00
GMAXP 8.00
KPOINT_MESH 8 8 8
KPOINT_SHIFT OFF OFF OFF
WIDTH 0.0002
EDELTA 0.5000D-09
NEG 8
CELL 10.30 10.30 10.30 90.00 90.00 90.00
&ATOMIC_SPECIES
Si 28.0900 pot.Si_pbe1
&END
&ATOMIC_COORDINATES CRYSTAL
0.000000000000 0.000000000000 0.000000000000 1 1 1
0.250000000000 0.250000000000 0.250000000000 1 1 1
&END
同様にsb100で使用するジョブスクリプトqsub_sb100.shも確認しましょう。 $ cat qsub_sb100.sh #$ -S /bin/sh
#$ -cwd
#$ -q sb.q
#$ -pe x6 6
#$ -N Si
module load intel/2021.2.0
module load intelmpi/2021.2.0
# Disable OPENMP parallelism
export OMP_NUM_THREADS=1
unset I_MPI_TCP_NETMASK
# Set the execuable of the STATE code
ln -fs ${HOME}/STATE/src/state-5.6.9/src/STATE .
# Set the pseudopotential data
ln -fs ../gncpp/pot.Si_pbe1
# Set the input/output file
INPUT_FILE=nfinp_scf
OUTPUT_FILE=nfout_scf
# Run!
mpirun ./STATE < ${INPUT_FILE} > ${OUTPUT_FILE}
結晶構造の最適化は以下のような手順で行います。
エネルギーの最小値を与える体積（格子定数）を平衡体積（平衡格子定数）といいます。 CELL 10.30 10.30 10.30 90.00 90.00 90.00 ここで2-4コラムはBohr半径を単位とした格子定数A, B, Cを表し、5-7コラムは格子ベクトルA2とA3, A3とA1, A1とA2のなす角度（度）を表しています。
入力ファイルが準備できたらジョブスクリプト（qsub_sb100.sh）の以下の行を編集します。 # Set the input/output file INPUT_FILE=nfinp_scf OUTPUT_FILE=nfout_scf
(これを10.50まで繰り返す)~
収束した全エネルギーは出力ファイルの以下の部分の TOTAL ENERGY AND ITS COMPONENTS
TOTAL ENERGY = -7.87286120 A.U.
KINETIC ENERGY = 3.05638778 A.U.
HARTREE ENERGY = 0.53752775 A.U.
XC ENERGY = -2.41643958 A.U.
LOCAL ENERGY = -0.76652974 A.U.
NONLOCAL ENERGY = 0.16607187 A.U.
EWALD ENERGY = -8.44987928 A.U.
PC ENERGY = 0.00000000 A.U.
ENTROPIC ENERGY = 0.00000000 A.U.
TOTAL ENERGYの行にA.U./Hartreeで示されていますので、nfout_scf_10.20-10.50の全ての全エネルギーを取り出すと以下のようになります。 -7.87286120 -7.87332156 -7.87355833 -7.87364725 -7.87354557 -7.87325835 -7.87284975 格子定数の関数として表示すると以下のようになります。 10.20 -7.87286120 10.25 -7.87332156 10.30 -7.87355833 10.35 -7.87364725 10.40 -7.87354557 10.45 -7.87325835 10.50 -7.87284975 これをgnuplotで表示すると以下のようになります。  結晶の体積は出力ファイルの"PRIM. CELL VOLUME"に以下のように表示されています（Bohr半径の三乗）。 PRIM. CELL VOLUME : 0.265302E+03 この情報を使うと全エネルギーを体積の関数として以下のように表示することも可能です。 0.265302E+03 -7.87286120 0.269223E+03 -7.87332156 0.273182E+03 -7.87355833 0.277179E+03 -7.87364725 0.281216E+03 -7.87354557 0.285292E+03 -7.87325835 0.289406E+03 -7.87284975 データを取り出すのが面倒になったらstate2ev.shというユーティリティーを使ってみましょう。 $ state2ev.sh nfout_scf_* 出力データをetot.datというファイルに保存します。 $ state2ev.sh nfout_scf_* > etot.dat これをgnuplotで表示します。  最後に得られたエネルギー・体積のデータをMurnaghan状態方程式にフィットします。ここではeosfitというユーティリティーを使いましょう。 $ eosfit etot.dat これを実行するとフィットされたエネルギー・体積の曲線が以下のように表示されます。  またフィッティングによって得られたパラメーターがeosfit.paramに出力されます。 |
