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

• バックアップ一覧
• 差分 を表示
• 現在との差分 を表示
• ソース を表示
• 計算コードの使い方/STATE/Si_opt へ行く。
  • 1 (2019-01-02 (水) 22:25:21)
  • 2 (2019-01-03 (木) 00:04:14)
  • 3 (2019-01-03 (木) 18:24:03)

シリコン結晶 †

この例ではダイアモンド構造のシリコンを取り上げ、セル（格子定数）の最適化および全エネルギーと格子定数の平面波カットオフ依存性を調べる方法を記述する。擬ポテンシャルにはSi_pbe (ノルム保存擬型) を用いる。

セルの最適化 †

現在のSTATEにはストレスの計算は実装されておらず、単位胞の形状および体積を自動で最適化することはできない。そのため全エネルギーを格子定数の関数として計算し、フィッティングにより平衡格子定数および体積を決定する必要がある。以下にその手順を記述する。

ここでは以下のスクリプトを使用し、自動的に入力ファイルを生成し、STATEを実行する。

#!/bin/sh
# 
# Number of CPUs
# 
NCPU=4
# 
# 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
#
# Minimum and maximum lattice parameter and mesh
AMIN=10.25
AMAX=10.45
MA=20
#
for (( i=0; i <= ${MA}; i++ ))
do
#
# 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

を実行する。