Skip to content

Instantly share code, notes, and snippets.

@t-nissie
Last active November 9, 2022 19:41
Show Gist options
  • Star 4 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save t-nissie/662ce194be356b5456c563195a052427 to your computer and use it in GitHub Desktop.
Save t-nissie/662ce194be356b5456c563195a052427 to your computer and use it in GitHub Desktop.
Calculation of SnO2 with Quantum Espresso

Calculation of SnO2 with Quantum Espresso

SnO2 is a crystal which have rutile structure.

Original files are in https://gist.github.com/t-nissie/662ce194be356b5456c563195a052427/ .

Preparation

  • Clone files in the gist as git clone https://gist.github.com/662ce194be356b5456c563195a052427.git SnO2-LDA-vc-relax-03-gist
  • cd SnO2-LDA-vc-relax-03-gist.
  • Get O.pw-mt_fhi.UPF and Sn.pw-mt_fhi.UPF form http://www.quantum-espresso.org/pseudopotentials/ .
  • Make symlinks to bands.x, dos.x, plotband.x and pw.x in the current directory, i.e. ln -s /SOMEWHERE/espresso-5.4.0/bin/bands.x.

relaxation of cell and atomic positions

Using an input file of rutile.pw.in, we perform optimization of cell parameters and atomic positions. If you have 12 cores in your computer,

$ mpirun -np 12 ./pw.x -in rutile.pw.in | tee rutile.pw.out
$ grep 'unit-cell volume  ' rutile.pw.out | tail -1
     unit-cell volume          =     429.5592 (a.u.)^3
$ grep '!    total energy' rutile.pw.out | tail -1
!    total energy              =    -142.69934279 Ry
$ grep -B13 'End final coordinates' rutile.pw.out

CELL_PARAMETERS (alat=  8.84000000)
   0.974988481   0.000000000   0.000000000
   0.000000000   0.974988481   0.000000000
   0.000000000   0.000000000   0.654134780

ATOMIC_POSITIONS (crystal)
Sn       0.000000000   0.000000000   0.000000000
Sn       0.500000000   0.500000000   0.500000000
O        0.305707826   0.305707826   0.000000000
O        0.694292174   0.694292174   0.000000000
O        0.194292174   0.805707826   0.500000000
O        0.805707826   0.194292174   0.500000000
End final coordinates

Bulk modulus

Convert the unit of bulk modulus by yourself.

$ ./rutile.pw.sh
$ gnuplot eos.gp
       :
Final set of parameters            Asymptotic Standard Error
=======================            ==========================

B0              = 0.01641          +/- 4.396e-06    (0.02679%)
B0p             = 4.99822          +/- 0.02974      (0.595%)
V0              = 430.014          +/- 0.004041     (0.0009397%)
Emin            = -142.699         +/- 1.232e-06    (8.637e-07%)

$ gv -watch eos.eps &

If you do not have gv(1), you can use evince(1) instead for preview of generated eos.eps.

eos.jpg

SCF calculation

Although the last calculation with rutile.pw.in gives SCF results in its end, we perform an SCF calculation with rutile.scf.in. Optimized cell parameters and atom positions were copied from rutile.pw.out into rutile.scf.in.

$ mpirun -np 12 ./pw.x -in rutile.scf.in | tee rutile.scf.out

Confirm that stress and forces are small enough.

You will get a directory of rutile.save/.

Plot band structure

In rutile.bands.in, crystal_b is used. How to write input files of rutile.bands.in and bands.in is described around

Anyway, we do non-SCF calculations at each k-point along A-M-Gamma-X-M-Gamma-Z-R-A-Z with pw.x, and then we execute bands.x with bands.sh as

$ submit EJCF general -jcf rutile.bands.jcf
$ gv rutile.bands.ps &

If you do not have gv(1), you can use evince(1) instead for preview of generated rutile.bands.ps.

Plot DOS

$ submit EJCF general -jcf rutile.dos.jcf
$ gnuplot
gnuplot> set ytics 2
gnuplot> plot [0:10] [-10:14] 'rutile.dos' u 2:($1-7.2) w l lw 2

SnO2.band.png

Plotted electronic dispersion and DOS of SnO2. Spikes in the band plot may be due to a bug in pw.x.

ToDo

  • To check convergence w.r.t. k-points
  • To check convergence w.r.t. kinetic energy cutoff

Refereces

  • J. Robertson: J. Phys. C: Solid State Phys. 12, 4767 (1979).
&bands
prefix = 'rutile'
outdir='./'
filband = 'bands.dat'
/
&dos
outdir = './',
prefix='rutile',
fildos='rutile.dos',
DeltaE=0.1
ngauss=0,
degauss=0.01
/
#!/usr/bin/env gnuplot
##
!./eos.sh > eos.dat
set terminal postscript eps enhanced "Times-Roman" color 14
set encoding iso_8859_1
set output "eos.eps"
set size 0.5, 0.5
set xlabel 'volume [Bohr^3]'
set ylabel '{/Times-Italic E}-{/Times-Italic E}_0 [Ry]'
set yrange [-0.001:0.010]
set ytics 0.001
set format y "%.3f"
Etot(x) = B0*x/(B0p*(B0p-1.0)) *(B0p*(1.0-V0/x)+(V0/x)**B0p-1.0)+Emin
# initial values
B0 = 0.02
B0p = 5.6
V0 = 432.0
Emin = -143.014
fit Etot(x) 'eos.dat' using 4:10 via B0, B0p, V0, Emin
plot 'eos.dat' using 4:($10-Emin) t 'calculated data' w p pt 1 ps 2 lw 2,\
Etot(x)-Emin t 'equation of state' w l lt 1 lw 2
set output
!epstopdf eos.eps
#Local variables:
# compile-command: "gnuplot eos.gp"
#End:
#!/bin/sh
##
for d in *[0-9].pw.out; do
energy=`grep '! total' $d | tail -n 1`
volume=`grep ' unit-cell volume ' $d | tail -n 1`
echo $volume $energy
done
&control
calculation = 'bands'
prefix = 'rutile'
pseudo_dir = './'
outdir = './'
/
&system
ibrav = 6
celldm(1) = 8.61889817204
celldm(3) = 0.670915393
space_group = 136
nat = 2
ntyp = 2
ecutwfc = 120.0
ecutrho = 1200.0
nbnd=40
/
&electrons
/
ATOMIC_SPECIES
Sn 118.71 Sn.pw-mt_fhi.UPF
O 16.00 O.pw-mt_fhi.UPF
ATOMIC_POSITIONS (crystal_sg)
Sn 0.000000000 0.000000000 0.000000000
O 0.305707826 0.305707826 0.000000000
K_POINTS (crystal_b)
10
0.500000 0.500000 0.500000 50 !A
0.500000 0.500000 0.000000 50 !M
0.000000 0.000000 0.000000 50 !Gamma
0.500000 0.000000 0.000000 50 !X
0.500000 0.500000 0.000000 50 !M
0.000000 0.000000 0.000000 50 !Gamma
0.000000 0.000000 0.500000 50 !Z
0.500000 0.000000 0.500000 50 !R
0.500000 0.500000 0.500000 50 !A
0.000000 0.000000 0.500000 50 !Z
#!/usr/bin/csh -f
#@ job_name = rutile.bands
#@ output = $(job_name).$(jobid).$(stepid).out
#@ error = $(job_name).$(jobid).$(stepid).err
#@ notification = never
#@ job_type = parallel
#@ total_tasks = 32
#@ stack_limit = 960mb
#@ task_affinity = cpu(1)
#@ cpus_per_core = 1
#@ queue
# Usage: submit EJCF general -jcf rutile.bands.jcf
##
poe ./pw.x < rutile.bands.in > rutile.bands.out
poe ./bands.x < bands.in > bands.out
./plotband.x <<EOF
bands.dat
-2.81, 21.2
rutile.bands.xmgr
rutile.bands.ps
7.2
2., 7.2
EOF
rm -f rutile.wfc*
&control
calculation = 'nscf'
prefix = 'rutile'
pseudo_dir = './'
outdir = './'
wf_collect=.false.
/
&system
ibrav = 6
celldm(1) = 8.61889817204
celldm(3) = 0.670915393
space_group = 136
nat = 2
ntyp = 2
ecutwfc = 120.0
ecutrho = 1200.0
nbnd=40
/
&electrons
conv_thr = 1.0d-8
/
&IONS
/
&cell
cell_factor = 3.0d0
/
ATOMIC_SPECIES
Sn 118.71 Sn.pw-mt_fhi.UPF
O 16.00 O.pw-mt_fhi.UPF
K_POINTS AUTOMATIC
8 8 16 1 1 1
ATOMIC_POSITIONS (crystal_sg)
Sn 0.000000000 0.000000000 0.000000000
O 0.305707826 0.305707826 0.000000000
#!/usr/bin/csh -f
#@ job_name = rutile.dos
#@ output = $(job_name).$(jobid).$(stepid).out
#@ error = $(job_name).$(jobid).$(stepid).err
#@ notification = never
#@ job_type = parallel
#@ total_tasks = 32
#@ stack_limit = 960mb
#@ task_affinity = cpu(1)
#@ cpus_per_core = 1
#@ queue
# Usage: submit EJCF general -jcf rutile.dos.jcf
##
poe ./pw.x < rutile.dos.in > rutile.dos.out
poe ./dos.x < dos.in > dos.out
&control
calculation = 'vc-relax'
restart_mode = 'from_scratch'
prefix = 'rutile-8.84.pw'
pseudo_dir = './'
outdir = './'
wf_collect=.true.
forc_conv_thr = 1d-4
disk_io = 'none'
/
&system
ibrav = 0
celldm(1) = 8.84
nat = 6
ntyp = 2
ecutwfc = 120.0
ecutrho = 1200.0
/
&electrons
conv_thr = 1.0d-8
/
&IONS
/
&cell
cell_factor = 3.0d0
cell_dofree = 'all'
/
ATOMIC_SPECIES
Sn 118.71 Sn.pw-mt_fhi.UPF
O 16.00 O.pw-mt_fhi.UPF
K_POINTS AUTOMATIC
8 8 16 1 1 1
ATOMIC_POSITIONS (crystal)
Sn 0.000000000 0.000000000 0.000000000
Sn 0.500000000 0.500000000 0.500000000
O 0.303875023 0.303875023 0.000000000
O 0.696124977 0.696124977 0.000000000
O 0.196124977 0.803875023 0.500000000
O 0.803875023 0.196124977 0.500000000
CELL_PARAMETERS
0.989189687 0.000000000 0.000000000
0.000000000 0.989189687 0.000000000
0.000000000 0.000000000 0.634961507
#!/usr/bin/csh -f
#@ job_name = rutile.pw
#@ output = $(job_name).$(jobid).$(stepid).out
#@ error = $(job_name).$(jobid).$(stepid).err
#@ notification = never
#@ job_type = parallel
#@ total_tasks = 32
#@ stack_limit = 960mb
#@ task_affinity = cpu(1)
#@ cpus_per_core = 1
#@ queue
# Usage: submit EJCF general -jcf rutile.pw.jcf
##
setenv MPI_COMMAND poe
poe ./pw.x < rutile.pw.in > rutile.pw.out
./rutile.pw.sh
#!/bin/sh
##
if [ -z "$MPI_COMMAND" ]; then
MPI_COMMAND='mpiexec -np 12'
fi
for a in `seq -w 8.70 0.02 9.00`; do
cat <<EOF > $a.pw.in
&control
calculation = 'vc-relax'
restart_mode = 'from_scratch'
prefix = 'rutile-$a.pw'
pseudo_dir = './'
outdir = './'
wf_collect=.true.
forc_conv_thr = 1d-4
disk_io = 'none'
/
&system
ibrav = 0
celldm(1) = $a
!celldm(3) = .64410104
nat = 6
ntyp = 2
!occupations = 'fixed'
ecutwfc = 120.0
ecutrho = 1200.0
!input_dft = 'vdW-DF-C09'
!nosym = .TRUE.
/
&electrons
conv_thr = 1.0d-8
/
&IONS
/
&cell
cell_factor = 3.0d0
cell_dofree = 'shape'
/
ATOMIC_SPECIES
Sn 118.71 Sn.pw-mt_fhi.UPF
O 16.00 O.pw-mt_fhi.UPF
K_POINTS AUTOMATIC
8 8 16 1 1 1
ATOMIC_POSITIONS (crystal)
Sn 0.000000000 0.000000000 0.000000000
Sn 0.500000000 0.500000000 0.500000000
O 0.303875023 0.303875023 0.000000000
O 0.696124977 0.696124977 0.000000000
O 0.196124977 0.803875023 0.500000000
O 0.803875023 0.196124977 0.500000000
CELL_PARAMETERS
0.989189687 0.000000000 0.000000000
0.000000000 0.989189687 0.000000000
0.000000000 0.000000000 0.634961507
EOF
$MPI_COMMAND ./pw.x < $a.pw.in > $a.pw.out
done
rm -rf *.save *wfc*
&control
calculation = 'scf'
restart_mode = 'from_scratch'
tstress = .true.
tprnfor = .true.
prefix = 'rutile'
pseudo_dir = './'
outdir = './'
wf_collect=.true.
/
&system
ibrav = 6
celldm(1) = 8.61889817204
celldm(3) = 0.670915393
space_group = 136
nat = 2
ntyp = 2
ecutwfc = 120.0
ecutrho = 1200.0
/
&electrons
conv_thr = 1.0d-8
/
&IONS
/
&cell
cell_factor = 3.0d0
/
ATOMIC_SPECIES
Sn 118.71 Sn.pw-mt_fhi.UPF
O 16.00 O.pw-mt_fhi.UPF
K_POINTS AUTOMATIC
8 8 16 1 1 1
ATOMIC_POSITIONS (crystal_sg)
Sn 0.000000000 0.000000000 0.000000000
O 0.305707826 0.305707826 0.000000000
&control
calculation = 'scf'
restart_mode = 'from_scratch'
tstress = .true.
tprnfor = .true.
prefix = 'rutile.scf'
pseudo_dir = './'
outdir = './'
wf_collect=.true.
forc_conv_thr = 1d-4
disk_io = 'none'
/
&system
ibrav = 0
celldm(1) = 8.61889817204
nat = 6
ntyp = 2
ecutwfc = 120.0
ecutrho = 1200.0
/
&electrons
conv_thr = 1.0d-8
/
&IONS
/
&cell
cell_factor = 3.0d0
/
ATOMIC_SPECIES
Sn 118.71 Sn.pw-mt_fhi.UPF
O 16.00 O.pw-mt_fhi.UPF
K_POINTS AUTOMATIC
8 8 16 1 1 1
ATOMIC_POSITIONS (crystal)
Sn 0.000000000 0.000000000 0.000000000
Sn 0.500000000 0.500000000 0.500000000
O 0.305707826 0.305707826 0.000000000
O 0.694292174 0.694292174 0.000000000
O 0.194292174 0.805707826 0.500000000
O 0.805707826 0.194292174 0.500000000
CELL_PARAMETERS
1.000000000 0.000000000 0.000000000
0.000000000 1.000000000 0.000000000
0.000000000 0.000000000 0.670915393
#!/usr/bin/csh -f
#@ job_name = rutile.scf
#@ output = $(job_name).$(jobid).$(stepid).out
#@ error = $(job_name).$(jobid).$(stepid).err
#@ notification = never
#@ job_type = parallel
#@ total_tasks = 32
#@ stack_limit = 960mb
#@ task_affinity = cpu(1)
#@ cpus_per_core = 1
#@ queue
# Usage: submit EJCF general -jcf rutile.scf.jcf
##
poe ./pw.x < rutile.scf.in > rutile.scf.out
@hamzaal01407
Copy link

Thank you for your time and thank you especially for the tutorial (very detailed).

I have a question about the level of Fermi:
You had drawn the two graphs "DOS + Bands", with a displacement (with the shifting ) in the graphs with a value equal to the value of Fermi.

I calculated using your INPUT files, I found fremi = ### 9.1874 eV
but I believe that your Fermi value is a little lower. What is the Fermi value of your result ??

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment