STATE_tutorial/Si_jp - 第一原理分子動力学プログラム STATE Senri Wiki

シリコン結晶 †

このチュートリアルでは固体の計算におけるカットオフエネルギーとk点についての収束性のテスト、およびセルパラメーターの最適化の方法について説明します。

収束性のテスト †

以下では全エネルギーのカットオフエネルギーについての収束性を調べます。ここで格子定数は10.2347 Bohrとしブリルアンゾーン積分は 4 x 4 x 4 Monkhorst-Pack k点メッシュを用います。ノルム保存型擬ポテンシャル (pot.Si_pbe1) を使用し、電荷密度の平面波展開のカットオフエネルギーは波動関数のものの4倍（カットオフ波数ベクトルは2倍）になっています。

入力ファイル

 0  0  0  0  0  0                      : INPUT_CTRL(1:6) (DUMMY)
 5.00  10.00 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

上記の入力ファイルをテンプレートとして、全エネルギーをGMAXとGMAXP (GMAXP=GMAX*2) の関数として計算すると以下が得られます。

#GMAX   ETOT (Ha)
 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

↑

ユニットセルの最適化 †

現在のSTATEではストレステンソルの計算がまだ実装されていないので、セルの最適化は手動で行う必要があります。全エネルギーをセルパラメーターの関数として計算し、適当な関数にフィッティングすることで平衡格子定数（セルパラメーター）が得られます。この例ではカットオフ波数ベクトル（カットオフエネルギー）は6 (36 Ry)、k点は4 x 4 x 4 Monkhorst-Packメッシュを採用し、以下のようなスクリプトを使って計算を自動化しています。

#!/bin/sh
#
# Path to STATE executable and pseudopotential(s)
#
ROOTDIR=${HOME}/STATE/src/state-5.6.6/
SRCDIR=${ROOTDIR}/src
PPDIR=${HOME}/STATE/gncpp
MPIRUN="mpirun -np 4"
STATE=STATE
#
# Prepare the STATE executable
#
ln -fs ${SRCDIR}/${STATE}
#
# Prepare the pseudopotential data
#
ln -fs ${PPDIR}/pot.Si_pbe1 fort.37
#
# Minimum and maximum lattice parameter and mesh
#
AMIN=10.25
AMAX=10.45
MA=20
#
for (( i=0; i <= ${MA}; i++ ))
do
#
# Lattice parameter
#
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

このスクリプトは使用する環境（キューイングシステム）によって変更する必要があります。このスクリプトを用いることで以下のような体積の関数としての全エネルギーが得られます（格子定数ではないことに注意）。

#volume(Bohr^3) ftot(Ha)
0.269223E+03 -7.87884234
0.270011E+03 -7.87890143
0.270802E+03 -7.87895926
0.271593E+03 -7.87900043
0.272387E+03 -7.87904967
0.273182E+03 -7.87907649
0.273978E+03 -7.87910231
0.274776E+03 -7.87911557
0.275576E+03 -7.87912945
0.276377E+03 -7.87912776
0.277179E+03 -7.87912853
0.277984E+03 -7.87911256
0.278789E+03 -7.87910050
0.279597E+03 -7.87907182
0.280406E+03 -7.87905614
0.281216E+03 -7.87901302
0.282028E+03 -7.87897140
0.282842E+03 -7.87891793
0.283657E+03 -7.87886603
0.284473E+03 -7.87879889
0.285292E+03 -7.87873581

このデータは添付のgnuplotのスクリプトeosfitを使用してMurnaghan状態方程式にフィットすることで平衡体積の276.381913141324 (Bohr^3) (40.9555250048756 Angstrom^3)と平衡格子定数 10.34 Bohr (5.47 Angstrom) が得られます。この値はk点およびカットオフエネルギーについても収束性を調べる必要があります。

↑

演習問題 †

平衡格子定数のk点およびカットオフエネルギーについての収束性を調べてみましょう。

第一原理分子動力学プログラム STATE Senri Wiki

Menu

シリコン結晶 †

収束性のテスト †

ユニットセルの最適化 †

演習問題 †