How to do a lattice parameter optimization of Silicon using script

We need 2 files

algerien1970@linux-wipm:~/abinitio/QE-tutorials/Si_eos> ls
 run_si_eos  Si.pz-vbc.UPF

Script run_si_eos

#!/bin/sh
# reminder: from now on, what follows the character # is a comment
####################################################################
#
# define the following variables according to your needs
#
outdir=OUT_EOS
pseudo_dir=./
# the following is not actually used:
# espresso_dir=top_directory_of_espresso_package
####################################################################

rm -f si.eos.out si.etot_vs_alat
touch si.etot_vs_alat

for alat in 9.8 9.9 10.0 10.1 10.2 10.3 10.4 10.5 10.6 10.7 ; do

# self-consistent calculation
cat > si.eos.in << EOF
 &control
    prefix='silicon',
    pseudo_dir = '$pseudo_dir/',
    outdir='$outdir/'
 /
 &system    
    ibrav=  2, celldm(1) =$alat, nat=  2, ntyp= 1,
    ecutwfc = 20.0, 
 /
 &electrons
 /
ATOMIC_SPECIES
 Si  28.086  Si.pz-vbc.UPF
ATOMIC_POSITIONS
 Si 0.00 0.00 0.00 
 Si 0.25 0.25 0.25 
K_POINTS automatic
   4 4 4 1 1 1
EOF

# If pw.x is not found, specify the correct value for $espresso_dir,
# use $espresso_dir/bin/pw.x instead of pw.x

pw.x < si.eos.in |tee si.eos.out

grep -e 'lattice parameter' -e ! si.eos.out | \
      awk '/lattice/{alat=$(NF-1)}/!/{print alat, $(NF-1)}' >> si.etot_vs_alat

done

Pseudopotential Si.pz-vbc.UPF

You can download the file using the following command

$ wget https://www.quantum-espresso.org/upf_files/Si.pz-vbc.UPF --no-check-certificate

Optimization

algerien1970@linux-wipm:~/abinitio/QE-tutorials/Si_eos> chmod +x run_si_eos

algerien1970@linux-wipm:~/abinitio/QE-tutorials/Si_eos> ./run_si_eos

We get the following file si.etot_vs_alat

algerien1970@linux-wipm:~/abinitio/QE-tutorials/Si_eos> ls
OUT_EOS  run_si_eos  si.eos.in  si.eos.out  si.etot_vs_alat  Si.pz-vbc.UPF

si.etot_vs_alat

9.8000 -15.83345937
9.9000 -15.83971334
10.0000 -15.84397984
10.1000 -15.84661822
10.2000 -15.84754590
10.3000 -15.84700488
10.4000 -15.84512699
10.5000 -15.84209101
10.6000 -15.83789204
10.7000 -15.83273079

You can plot file si.etot_vs_alat using your preferred plotting program, for instance:

$ gnuplot
gnuplot > plot "si.etot_vs_alat" using 1:2 with lines

The experimental lattice parameter for Si is 5.47 Å, or 10.26 a.u.: this is a case where plain simple LDA yields remarkable results You may experiment changing cutoff, k-points, pseudopotential, ... You should find that

• The energy vs lattice parameter E(a) curves are shifted down rather uniformly with increasing cutoff and are not strongly dependent on k-points.
• Structural properties and energy differences converge faster that total energies.

Use code ev.x to fit your results to a phenomenological EOS (e.g. Murnaghan) and to get accurate values for the lattice parameter and for the bulk modulus. The code prompts for some data, reads a file like si.etot vs alat: for cubic systems, rows ai E(ai )

ev.x :

Run interatively. It reads a file containing 2 columns: lattice parameter and
corresponding total energy, it gives in output optimized lattice parameter,
bulk modulus and more.
Contributions by Eyvaz Isaev

Running code ev.x

algerien1970@linux-wipm:~/abinitio/QE-tutorials/Si_eos> ev.x
     Lattice parameter or Volume are in (au, Ang) > au
     Enter type of bravais lattice (fcc, bcc, sc, noncubic) > fcc
     Enter type of equation of state :
     1=birch1, 2=birch2, 3=keane, 4=murnaghan > 4
     Input file > si.etot_vs_alat
 Minimization succeeded
     Output file > si.lattice_out

The content of the file si.lattice_out

# equation of state: murnaghan.        chisq =   0.7026D-09
# a0 = 10.2129 a.u., k0 =  921 kbar, dk0 =  4.20 d2k0 =  0.000 emin =  -15.84753
# a0 =  5.40445 Ang, k0 =  92.1 GPa,  V0 =   266.31 (a.u.)^3,  V0 =   39.46 A^3 

#########################################################################
# Lat.Par       E_calc        E_fit       E_diff    Pressure      Enthalpy
# a.u.            Ry           Ry            Ry        GPa           Ry
#########################################################################
  9.80000     -15.83346     -15.83345    -0.00001      14.96      -15.59421
  9.90000     -15.83971     -15.83972     0.00001      10.53      -15.66610
 10.00000     -15.84398     -15.84403     0.00005       6.67      -15.73066
 10.10000     -15.84662     -15.84658    -0.00004       3.30      -15.78887
 10.20000     -15.84755     -15.84752    -0.00003       0.35      -15.84118
 10.30000     -15.84700     -15.84700    -0.00001      -2.22      -15.88832
 10.40000     -15.84513     -15.84515     0.00002      -4.48      -15.93085
 10.50000     -15.84209     -15.84209    -0.00000      -6.47      -15.96931
 10.60000     -15.83789     -15.83791     0.00002      -8.21      -16.00403
 10.70000     -15.83273     -15.83271    -0.00002      -9.74      -16.03550