![[PukiWiki]](image/logo_state-senri3_80x80.png "[PukiWiki]"){kind=logo}
第一原理分子動力学プログラム STATE Senri Wiki
- バックアップ一覧
- 差分 を表示
- 現在との差分 を表示
- ソース を表示
- STATE_tutorial/Si へ行く。
- 1 (2019-02-12 (火) 22:01:40)
- 2 (2019-02-12 (火) 22:04:46)
- 3 (2019-02-16 (土) 15:27:08)
- 4 (2019-09-05 (木) 10:26:05)
- 5 (2019-09-07 (土) 09:43:38)
Crystalline silicon †
This tutorial explains how to perform convergence study with respect to the cutoff energy and the number of k-points, and how to obtain the equilibrium cell parameters.
Convergence study †
Let us investigate the convergence of the total energy with respect to the cutoff energy. Here lattice constant of 10.2347 Bohr and 4 x 4 x 4 Monkhorst-Pack k-point mesh are used. Normconserving pseudopotential (pot.Si_pbe1) is used and the cutoff energy for the augmentation charge is 4 times larger than that for the wave functions.
- Input file
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
Using the above input file as a template, we obtain the following total energy as a function of the cutoff energy#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
Cell optimization †
In the current version of STATE, the stress tensor is not yet implemented (there was an implementation in the ancestor of STATE) and thus, the cell parameters should be optimized manually. The equilibrium lattice parameter is obtained by calculating the total energy as a function of cell parameter and by fitting to a certain function. In this example, we use the cutoff wave vector (cutoff energy) of 6 (36 Ry), 4x4x4 Monkhorst-Pack k-point mesh. To automate the calculations, I usually write a script like this:
#!/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
The script differs depending on the environment (queueing system), but the basic is the same.
Using the script, one may obtain the following total energy as a function of _volume_.
#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
This total energy is fittted to the Murnaghan equation of state using a gnuplot script "eosfit" as attached, to obtain the equilbrium volume of 276.381913141324 (Bohr^3) (40.9555250048756 Angstrom^3) , corresponding to the lattice constant of 10.34 Bohr (5.47 Angstrom). For more critical assessment, the convergence of the lattice constant with respect to the cutoff energy and the number of k-points are studied.
Exercise †
- Perform a convergence study of the equilibrium lattice constant on the cutoff energy and the numer of k-points.
