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.
トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS