Part 2: MgO¶
By the end of this tutorial, you will be able to:
- perform a relaxation of the cell volume
- monitor the convergence of the lattice parameter with the k points
4.1 Task¶
Monitor the convergence of the lattice parameter as a function of the k-points mesh density.
To obtain an accurate lattice parameter, it is important to have a sufficiently large basis set that is controlled by the ENCUT tag in the INCAR file. Additionally, sufficient sampling of the first Brillouin zone is required. It is recommended to check how the quantity of interest (in the current example the relaxed lattice parameter) changes with changing k-point density so that one can choose the smallest density that yields an accurate lattice parameter thus saving computational resources. This is important so that the subsequent more expensive calculations (force constants from finite differences in a supercell) are as computationally efficient as possible.
4.2 Input¶
The input files to run this example are prepared at $TUTORIALS/phonon/e04_*. Check them out!
MgO
4.211
0.0000000000 0.5000000000 0.5000000000
0.5000000000 0.0000000000 0.5000000000
0.5000000000 0.5000000000 0.0000000000
Mg O
1 1
Direct
0.0 0.0 0.0
0.5 0.5 0.5
PREC = Accurate
EDIFF = 1.E-08
ENCUT = 500
ISMEAR = 0
SIGMA = 0.05
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.
# Relax
IBRION = 2
ISIF = 3
NSW = 100
PAW_PBE Mg 13Apr2007
PAW_PBE O 08Apr2002
4.3 Calculation¶
We will run a few lattice relaxation calculations using different k-point meshes to check the convergence. This can be done using the following commands
cd $TUTORIALS/phonon/e04_*
for n in 3 4 6 8 10 12
do
echo "KPOINTS" > KPOINTS
echo "0" >> KPOINTS
echo "Gamma" >> KPOINTS
echo "$n $n $n" >> KPOINTS
mpirun -np 4 vasp_std
cp vaspout.h5 vaspout_$n.h5 #save vaspout.h5
done
These calculations take a while to run. While they are running have a look at the script above and the Python script below and try to understand what they will do.
Click to see the answer!
The script above creates a KPOINTS file specifying a N×N×N for N = {4, 6, 8, 10, 12}, launches a VASP calculation and then copies the vaspout.h5 file to vaspout_$N.h5
The Python script below will read the lattice parameter from different vaspout_$N.h5 files and plot them as a function of the k-point mesh.
Once the calculations are finished, run the following Python script
import py4vasp
kmeshes = [3, 4, 6, 8, 10, 12]
y = []
for k in kmeshes:
mycalc = py4vasp.Calculation.from_file( f"./e04_mgo-relax/vaspout_{k}.h5")
# The conventional cell has 4 primitive cells
y.append((4.0 * mycalc.structure.volume()) ** (1/3))
from py4vasp import plot
plot(kmeshes, y, "lattice constant", xlabel='k point mesh')
Looking at the graph above: what k-point mesh do you consider sufficient to converge the calculation? You can create a POSCAR file directly from your previous calculation with the following code:
#open the POSCAR file where we want to write the relaxed structure
with open("./e05_mgo-fin-diff-force-constants/POSCAR-relaxed",'w') as f:
mycalc = py4vasp.Calculation.from_file(f"./e04_mgo-relax/vaspout_6.h5")
f.write(mycalc.structure.to_POSCAR())
4.4 Questions¶
- Why should you converge the lattice parameter as a function of the k-point density?
By the end of this tutorial, you will be able to:
- create supercells for phonon calculations
- compute the force-constants using the finite differences approach
5.1 Task¶
We will compute the force-constants and phonon modes using the finite differences approach for the MgO-bulk crystal.
5.2 Input¶
First, we create POSCAR files with the desired supercell size. Since we have a bulk material, we study increasing the primitive cell along all three cartesian directions. The cubic symmetry allows us to use a common scaling factor along all directions. In the following Python code, we create supercells for 2×2×2, 3×3×3, and 4×4×4. We also create the KPOINTS files corresponding to a k-point sampling of 6x6x6 in the primitive cell. Execute the following code and verify that the POSCAR and KPOINTS files are what you would expect.
from pymatgen.core import Structure
from pymatgen.io.vasp import Poscar
my_struc = Structure.from_file("./e05_mgo-fin-diff-force-constants/POSCAR")
# make a 2x2x2 supercell
my_struc.make_supercell([2,2,2])
# write supercell to POSCAR format with specified filename
folder = "./e05_mgo-fin-diff-force-constants/supercell_2x2x2"
Poscar(my_struc).write_file(filename=f"{folder}/POSCAR",significant_figures=16)
# A POSCAR file will appear
# at $TUTORIALS/phonon/e05_mgo-fin-diff-force-constants/supercell_2x2x2
# You may need to refresh the file browser to see it
with open(f"{folder}/KPOINTS","w") as f:
f.write("KPOINTS\n0\nGamma\n3 3 3")
my_struc = Structure.from_file("./e05_mgo-fin-diff-force-constants/POSCAR")
my_struc.make_supercell([3,3,3])
folder = "./e05_mgo-fin-diff-force-constants/supercell_3x3x3"
Poscar(my_struc).write_file(filename=f"{folder}/POSCAR",significant_figures=16)
with open(f"{folder}/KPOINTS","w") as f:
f.write("KPOINTS\n0\nGamma\n2 2 2")
my_struc = Structure.from_file("./e05_mgo-fin-diff-force-constants/POSCAR")
my_struc.make_supercell([4,4,4])
folder = "./e05_mgo-fin-diff-force-constants/supercell_4x4x4"
Poscar(my_struc).write_file(filename=f"{folder}/POSCAR",significant_figures=16)
with open(f"{folder}/KPOINTS","w") as f:
f.write("KPOINTS\n0\nGamma\n2 2 2")
Click to see supercell_2x2x2/POSCAR!
Mg8 O8
1.0
0.0000000000000000 4.2458414324158662 4.2458414324158662
4.2458414324158662 0.0000000000000000 4.2458414324158662
4.2458414324158662 4.2458414324158662 0.0000000000000000
Mg O
8 8
direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 Mg
0.0000000000000000 0.0000000000000000 0.5000000000000000 Mg
0.0000000000000000 0.5000000000000000 0.0000000000000000 Mg
0.0000000000000000 0.5000000000000000 0.5000000000000000 Mg
0.5000000000000000 0.0000000000000000 0.0000000000000000 Mg
0.5000000000000000 0.0000000000000000 0.5000000000000000 Mg
0.5000000000000000 0.5000000000000000 0.0000000000000000 Mg
0.5000000000000000 0.5000000000000000 0.5000000000000000 Mg
0.2500000000000000 0.2500000000000000 0.2500000000000000 O
0.2500000000000000 0.2500000000000000 0.7500000000000000 O
0.2500000000000000 0.7500000000000000 0.2500000000000000 O
0.2500000000000000 0.7500000000000000 0.7500000000000000 O
0.7500000000000000 0.2500000000000000 0.2500000000000000 O
0.7500000000000000 0.2499999999999999 0.7500000000000000 O
0.7500000000000000 0.7500000000000000 0.2499999999999999 O
0.7500000000000001 0.7500000000000000 0.7500000000000000 O
Click to see supercell_3x3x3/POSCAR!
Mg27 O27
1.0
0.0000000000000000 6.3687621486237997 6.3687621486237997
6.3687621486237997 0.0000000000000000 6.3687621486237997
6.3687621486237997 6.3687621486237997 0.0000000000000000
Mg O
27 27
direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 Mg
0.0000000000000000 0.0000000000000000 0.3333333333333333 Mg
0.0000000000000000 0.0000000000000000 0.6666666666666666 Mg
0.0000000000000000 0.3333333333333333 0.0000000000000000 Mg
0.0000000000000000 0.3333333333333333 0.3333333333333333 Mg
0.0000000000000000 0.3333333333333333 0.6666666666666666 Mg
0.0000000000000000 0.6666666666666666 0.0000000000000000 Mg
0.9999999999999999 0.6666666666666666 0.3333333333333333 Mg
0.0000000000000000 0.6666666666666666 0.6666666666666666 Mg
0.3333333333333333 0.0000000000000000 0.0000000000000000 Mg
0.3333333333333333 0.0000000000000000 0.3333333333333333 Mg
0.3333333333333333 0.0000000000000000 0.6666666666666667 Mg
0.3333333333333333 0.3333333333333333 0.0000000000000000 Mg
0.3333333333333333 0.3333333333333333 0.3333333333333333 Mg
0.3333333333333333 0.3333333333333333 0.6666666666666667 Mg
0.3333333333333334 0.6666666666666666 0.9999999999999999 Mg
0.3333333333333333 0.6666666666666667 0.3333333333333333 Mg
0.3333333333333334 0.6666666666666666 0.6666666666666667 Mg
0.6666666666666666 0.0000000000000000 0.0000000000000000 Mg
0.6666666666666666 0.9999999999999999 0.3333333333333333 Mg
0.6666666666666666 0.0000000000000000 0.6666666666666666 Mg
0.6666666666666666 0.3333333333333334 0.9999999999999999 Mg
0.6666666666666667 0.3333333333333333 0.3333333333333334 Mg
0.6666666666666666 0.3333333333333334 0.6666666666666665 Mg
0.6666666666666666 0.6666666666666666 0.0000000000000000 Mg
0.6666666666666666 0.6666666666666666 0.3333333333333334 Mg
0.6666666666666666 0.6666666666666666 0.6666666666666666 Mg
0.1666666666666667 0.1666666666666667 0.1666666666666667 O
0.1666666666666667 0.1666666666666667 0.5000000000000000 O
0.1666666666666667 0.1666666666666667 0.8333333333333334 O
0.1666666666666667 0.4999999999999999 0.1666666666666666 O
0.1666666666666666 0.5000000000000000 0.4999999999999999 O
0.1666666666666667 0.4999999999999999 0.8333333333333335 O
0.1666666666666666 0.8333333333333334 0.1666666666666666 O
0.1666666666666667 0.8333333333333334 0.4999999999999999 O
0.1666666666666667 0.8333333333333333 0.8333333333333333 O
0.4999999999999999 0.1666666666666667 0.1666666666666666 O
0.5000000000000000 0.1666666666666666 0.5000000000000001 O
0.4999999999999999 0.1666666666666667 0.8333333333333333 O
0.5000000000000000 0.5000000000000000 0.1666666666666666 O
0.5000000000000000 0.5000000000000000 0.5000000000000000 O
0.5000000000000000 0.5000000000000000 0.8333333333333333 O
0.4999999999999999 0.8333333333333333 0.1666666666666666 O
0.5000000000000000 0.8333333333333333 0.5000000000000001 O
0.5000000000000000 0.8333333333333333 0.8333333333333331 O
0.8333333333333334 0.1666666666666666 0.1666666666666666 O
0.8333333333333334 0.1666666666666667 0.4999999999999999 O
0.8333333333333333 0.1666666666666667 0.8333333333333333 O
0.8333333333333333 0.4999999999999999 0.1666666666666667 O
0.8333333333333333 0.5000000000000000 0.4999999999999999 O
0.8333333333333333 0.5000000000000000 0.8333333333333331 O
0.8333333333333333 0.8333333333333333 0.1666666666666668 O
0.8333333333333333 0.8333333333333333 0.5000000000000000 O
0.8333333333333333 0.8333333333333333 0.8333333333333333 O
Click to see supercell_4x4x4/POSCAR!
Mg64 O64
1.0
0.0000000000000000 8.4916828648317324 8.4916828648317324
8.4916828648317324 0.0000000000000000 8.4916828648317324
8.4916828648317324 8.4916828648317324 0.0000000000000000
Mg O
64 64
direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 Mg
0.0000000000000000 0.0000000000000000 0.2500000000000000 Mg
0.0000000000000000 0.0000000000000000 0.5000000000000000 Mg
0.0000000000000000 0.0000000000000000 0.7500000000000000 Mg
0.0000000000000000 0.2500000000000000 0.0000000000000000 Mg
0.0000000000000000 0.2500000000000000 0.2500000000000000 Mg
0.0000000000000000 0.2500000000000000 0.5000000000000000 Mg
0.0000000000000000 0.2500000000000000 0.7500000000000000 Mg
0.0000000000000000 0.5000000000000000 0.0000000000000000 Mg
0.0000000000000000 0.5000000000000000 0.2500000000000000 Mg
0.0000000000000000 0.5000000000000000 0.5000000000000000 Mg
0.0000000000000000 0.4999999999999999 0.7500000000000000 Mg
0.0000000000000000 0.7500000000000000 0.0000000000000000 Mg
0.0000000000000000 0.7500000000000000 0.2500000000000000 Mg
0.9999999999999999 0.7500000000000001 0.5000000000000001 Mg
0.0000000000000000 0.7500000000000000 0.7500000000000000 Mg
0.2500000000000000 0.0000000000000000 0.0000000000000000 Mg
0.2500000000000000 0.0000000000000000 0.2500000000000000 Mg
0.2500000000000000 0.0000000000000000 0.5000000000000000 Mg
0.2500000000000000 0.0000000000000000 0.7500000000000000 Mg
0.2500000000000000 0.2500000000000000 0.0000000000000000 Mg
0.2500000000000000 0.2500000000000000 0.2500000000000000 Mg
0.2500000000000000 0.2500000000000000 0.5000000000000000 Mg
0.2500000000000000 0.2500000000000000 0.7500000000000000 Mg
0.2500000000000000 0.5000000000000000 0.0000000000000000 Mg
0.2500000000000000 0.5000000000000000 0.2500000000000000 Mg
0.2500000000000001 0.5000000000000000 0.5000000000000000 Mg
0.2500000000000000 0.5000000000000000 0.7500000000000000 Mg
0.2500000000000000 0.7500000000000000 0.0000000000000000 Mg
0.2500000000000000 0.7500000000000000 0.2500000000000000 Mg
0.2499999999999999 0.7500000000000001 0.5000000000000002 Mg
0.2500000000000000 0.7500000000000000 0.7500000000000000 Mg
0.5000000000000000 0.0000000000000000 0.0000000000000000 Mg
0.5000000000000000 0.0000000000000000 0.2500000000000000 Mg
0.5000000000000000 0.0000000000000000 0.5000000000000000 Mg
0.4999999999999999 0.0000000000000000 0.7500000000000000 Mg
0.5000000000000000 0.2500000000000000 0.0000000000000000 Mg
0.5000000000000000 0.2500000000000000 0.2500000000000000 Mg
0.5000000000000000 0.2500000000000000 0.5000000000000000 Mg
0.5000000000000000 0.2500000000000001 0.7500000000000000 Mg
0.5000000000000000 0.5000000000000000 0.0000000000000000 Mg
0.5000000000000000 0.5000000000000000 0.2500000000000000 Mg
0.5000000000000000 0.5000000000000000 0.5000000000000000 Mg
0.4999999999999999 0.5000000000000000 0.7500000000000000 Mg
0.4999999999999999 0.7500000000000000 0.0000000000000000 Mg
0.5000000000000000 0.7499999999999999 0.2500000000000001 Mg
0.4999999999999998 0.7500000000000001 0.5000000000000000 Mg
0.4999999999999999 0.7500000000000000 0.7500000000000000 Mg
0.7500000000000000 0.0000000000000000 0.0000000000000000 Mg
0.7500000000000000 0.0000000000000000 0.2500000000000000 Mg
0.7500000000000001 0.0000000000000000 0.5000000000000001 Mg
0.7500000000000001 0.0000000000000000 0.7500000000000000 Mg
0.7500000000000000 0.2500000000000000 0.0000000000000000 Mg
0.7500000000000000 0.2500000000000000 0.2500000000000000 Mg
0.7500000000000000 0.2499999999999999 0.5000000000000000 Mg
0.7500000000000000 0.2499999999999999 0.7500000000000000 Mg
0.7500000000000001 0.5000000000000000 0.9999999999999999 Mg
0.7500000000000001 0.5000000000000001 0.2499999999999999 Mg
0.7500000000000001 0.5000000000000000 0.4999999999999999 Mg
0.7500000000000001 0.5000000000000000 0.7499999999999999 Mg
0.7500000000000001 0.7500000000000000 0.0000000000000000 Mg
0.7500000000000000 0.7500000000000000 0.2499999999999999 Mg
0.7500000000000000 0.7500000000000001 0.5000000000000002 Mg
0.7500000000000001 0.7500000000000000 0.7500000000000000 Mg
0.1250000000000000 0.1250000000000000 0.1250000000000000 O
0.1250000000000000 0.1250000000000000 0.3750000000000000 O
0.1250000000000000 0.1250000000000000 0.6250000000000000 O
0.1250000000000000 0.1250000000000000 0.8750000000000000 O
0.1250000000000000 0.3750000000000000 0.1250000000000000 O
0.1250000000000000 0.3750000000000000 0.3750000000000000 O
0.1250000000000000 0.3749999999999999 0.6250000000000000 O
0.1250000000000000 0.3750000000000000 0.8750000000000000 O
0.1250000000000000 0.6250000000000000 0.1250000000000000 O
0.1250000000000000 0.6250000000000000 0.3750000000000000 O
0.1250000000000001 0.6250000000000000 0.6250000000000000 O
0.1250000000000001 0.6249999999999999 0.8749999999999998 O
0.1250000000000000 0.8750000000000000 0.1250000000000000 O
0.1250000000000000 0.8750000000000000 0.3750000000000000 O
0.1250000000000000 0.8750000000000000 0.6250000000000000 O
0.1250000000000000 0.8750000000000000 0.8750000000000000 O
0.3750000000000000 0.1250000000000000 0.1250000000000000 O
0.3750000000000000 0.1250000000000000 0.3750000000000000 O
0.3750000000000000 0.1250000000000000 0.6250000000000000 O
0.3749999999999999 0.1250000000000000 0.8750000000000000 O
0.3750000000000000 0.3750000000000000 0.1250000000000000 O
0.3750000000000001 0.3750000000000000 0.3750000000000000 O
0.3750000000000000 0.3750000000000000 0.6250000000000000 O
0.3750000000000000 0.3750000000000001 0.8750000000000000 O
0.3750000000000000 0.6250000000000000 0.1250000000000000 O
0.3750000000000000 0.6250000000000000 0.3750000000000000 O
0.3750000000000000 0.6250000000000000 0.6250000000000000 O
0.3750000000000001 0.6249999999999999 0.8749999999999998 O
0.3749999999999999 0.8750000000000000 0.1250000000000000 O
0.3750000000000000 0.8749999999999999 0.3750000000000001 O
0.3749999999999999 0.8750000000000000 0.6250000000000000 O
0.3749999999999999 0.8750000000000000 0.8750000000000000 O
0.6250000000000000 0.1250000000000000 0.1250000000000000 O
0.6250000000000000 0.1250000000000000 0.3750000000000000 O
0.6250000000000000 0.1250000000000001 0.6250000000000000 O
0.6250000000000000 0.1250000000000001 0.8750000000000000 O
0.6250000000000000 0.3749999999999999 0.1250000000000000 O
0.6250000000000000 0.3750000000000000 0.3750000000000000 O
0.6250000000000000 0.3750000000000000 0.6250000000000000 O
0.6249999999999999 0.3750000000000001 0.8750000000000000 O
0.6250000000000000 0.6250000000000000 0.1250000000000000 O
0.6250000000000000 0.6250000000000000 0.3750000000000001 O
0.6249999999999999 0.6250000000000000 0.6250000000000000 O
0.6250000000000000 0.6249999999999999 0.8749999999999998 O
0.6250000000000000 0.8749999999999999 0.1250000000000001 O
0.6249999999999999 0.8749999999999999 0.3750000000000001 O
0.6249999999999999 0.8750000000000000 0.6250000000000001 O
0.6249999999999999 0.8750000000000000 0.8750000000000001 O
0.8750000000000000 0.1250000000000000 0.1250000000000000 O
0.8749999999999999 0.1250000000000001 0.3750000000000000 O
0.8750000000000000 0.1249999999999999 0.6250000000000000 O
0.8750000000000000 0.1250000000000000 0.8750000000000000 O
0.8750000000000000 0.3749999999999999 0.1250000000000000 O
0.8749999999999999 0.3750000000000000 0.3750000000000001 O
0.8750000000000000 0.3749999999999999 0.6250000000000000 O
0.8750000000000000 0.3749999999999999 0.8750000000000000 O
0.8750000000000000 0.6250000000000000 0.1250000000000000 O
0.8750000000000000 0.6250000000000001 0.3749999999999999 O
0.8750000000000001 0.6250000000000000 0.6250000000000000 O
0.8750000000000002 0.6249999999999999 0.8750000000000000 O
0.8750000000000000 0.8750000000000000 0.1250000000000000 O
0.8750000000000000 0.8750000000000000 0.3750000000000000 O
0.8750000000000001 0.8750000000000000 0.6250000000000000 O
0.8750000000000000 0.8750000000000000 0.8750000000000000 O
PREC = Accurate
EDIFF = 1.E-08
ISMEAR = 0
SIGMA = 0.05
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.
# finite differences
IBRION = 6
POTIM = 0.015
Click to see supercell_2x2x2/KPOINTS!
K points
0
Gamma
3 3 3
0 0 0
Click to see supercell_3x3x3/KPOINTS!
K points
0
Gamma
2 2 2
0 0 0
Click to see supercell_4x4x4/KPOINTS!
K points
0
Gamma
2 2 2
0 0 0
PAW_PBE Mg 13Apr2007 PAW_PBE O 08Apr2002
5.3 Calculation¶
To compute the force constants using the finite differences approach we set IBRION=6 in the INCAR file and run VASP for the different supercells. This can be done using the following script
cd $TUTORIALS/phonon/e05_*
for folder in supercell_2x2x2 supercell_3x3x3 supercell_4x4x4
do
cd $folder
mpirun -np 4 vasp_std
cd ..
done
Note that the calculations for the 3×3×3 and 4×4×4 supercells take quite a while to run. In the next sections, we will use the calculations based on the 2×2×2 supercell keeping in mind that this is clearly insufficient to have an accurate phonon dispersion that can be compared with experimental data. You might choose to launch the 3×3×3 and 4×4×4 supercells now and leave them running on a separate terminal while you proceed with the rest of the tutorial. In the end, you might re-run the calculations of sections 6 and 8 using the force constants computed on larger supercells.
5.4 Questions¶
- What is the main difference between creating a supercell to compute the force constants for bulk material and a 2D material?
6.1 Task¶
Compute the phonon dispersion and phonon DOS for MgO.
Using the force constants computed for a 2×2×2 supercell in the previous run, we will now compute the phonon dispersion and phonon density of states. As outlined in the theory page about phonons, this is done by building the dynamical matrices at arbitrary q points by Fourier interpolating the force-constants from the supercell to the primitive cell.
6.2 Input¶
To compute the phonon dispersion, we require the force-constants on a supercell that we computed in the previous step.
These were written to the vaspout.h5 file.
When the LPHON_READ_FORCE_CONSTANTS is present in the INCAR file, VASP will read the force constants from the vaspin.h5 file, plot the phonon dispersion on the q points specified in the QPOINTS file and exit.
To do so copy the $TUTORIALS/phonon/e05_*/supercell_2x2x2/vaspout.h5 file generated in the previous run to $TUTORIALS/phonon/e06_mgo-dispersion/supercell_2x2x2/dispersion/vaspin.h5.
cp $TUTORIALS/phonon/e05_*/supercell_2x2x2/vaspout.h5 $TUTORIALS/phonon/e06_mgo-dispersion/supercell_2x2x2/dispersion/vaspin.h5
The additional input files are listed below
Mg8 O8
1.0
0.0000000000000000 4.2458414324158662 4.2458414324158662
4.2458414324158662 0.0000000000000000 4.2458414324158662
4.2458414324158662 4.2458414324158662 0.0000000000000000
Mg O
8 8
direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 Mg
0.0000000000000000 0.0000000000000000 0.5000000000000000 Mg
0.0000000000000000 0.5000000000000000 0.0000000000000000 Mg
0.0000000000000000 0.5000000000000000 0.5000000000000000 Mg
0.5000000000000000 0.0000000000000000 0.0000000000000000 Mg
0.5000000000000000 0.0000000000000000 0.5000000000000000 Mg
0.5000000000000000 0.5000000000000000 0.0000000000000000 Mg
0.5000000000000000 0.5000000000000000 0.5000000000000000 Mg
0.2500000000000000 0.2500000000000000 0.2500000000000000 O
0.2500000000000000 0.2500000000000000 0.7500000000000000 O
0.2500000000000000 0.7500000000000000 0.2500000000000000 O
0.2500000000000000 0.7500000000000000 0.7500000000000000 O
0.7500000000000000 0.2500000000000000 0.2500000000000000 O
0.7500000000000000 0.2499999999999999 0.7500000000000000 O
0.7500000000000000 0.7500000000000000 0.2499999999999999 O
0.7500000000000001 0.7500000000000000 0.7500000000000000 O
#computation of the phonon dispersion
PHON_NWRITE = -2
LPHON_DISPERSION = .TRUE.
LPHON_POLAR=.FALSE.
PHON_DOS=1
LPHON_READ_FORCE_CONSTANTS=.TRUE.
k-points for bandstructure using seekpath GAMMA-X-L-W-W_2-K-U
50
line
reciprocal
0.00000000 0.00000000 0.00000000 GAMMA
0.50000000 0.00000000 0.50000000 X
0.50000000 0.00000000 0.50000000 X
0.50000000 0.25000000 0.75000000 W
0.50000000 0.25000000 0.75000000 W
0.37500000 0.37500000 0.75000000 K
0.37500000 0.37500000 0.75000000 K
0.00000000 0.00000000 0.00000000 GAMMA
0.00000000 0.00000000 0.00000000 GAMMA
0.50000000 0.50000000 0.50000000 L
0.50000000 0.50000000 0.50000000 L
0.62500000 0.25000000 0.62500000 U
0.62500000 0.25000000 0.62500000 U
0.75000000 0.25000000 0.50000000 W_2
0.75000000 0.25000000 0.50000000 W_2
0.50000000 0.50000000 0.50000000 L
0.50000000 0.50000000 0.50000000 L
0.37500000 0.37500000 0.75000000 K
PAW_PBE Mg 13Apr2007
PAW_PBE O 08Apr2002
6.3 Calculation¶
Now, use VASP to generate the phonon dispersion:
cd $TUTORIALS/phonon/e06_mgo-dispersion/supercell_2x2x2/dispersion
mpirun -np 1 vasp_std
Plot the phonon dispersion using py4vasp!
import py4vasp
mycalc = py4vasp.Calculation.from_path("./e06_mgo-dispersion/supercell_2x2x2/dispersion")
mycalc.phonon_band.plot(selection="Mg O")
To check the convergence of the phonon DOS with respect to the q-point sampling, copy the force-constants to the dos folder:
cp $TUTORIALS/phonon/e05_*/supercell_2x2x2/vaspout.h5 $TUTORIALS/phonon/e06_mgo-dispersion/supercell_2x2x2/dos/vaspin.h5
and execute the following code in the folder $TUTORIALS/phonons/e06_mgo-dispersion/supercell_2x2x2/dos
cd $TUTORIALS/phonon/e06_*/supercell_2x2x2/dos
for i in 10 20 30 40
do
echo "QPOINTS" > QPOINTS
echo "0" >> QPOINTS
echo "Gamma" >> QPOINTS
echo "$i $i $i" >> QPOINTS
vasp_std
cp vaspout.h5 vaspout_$i.h5
done
This will compute the phonon density of states by interpolating the phonon frequencies in the q points specified by the different QPOINTS files with meshes of 10×0x×0, 20x×0x×0, 30x30x×0 ×nd 40x40x40 and save the results in different vaspout_*.h5 files.
Now, we can plot the phonon DOS on these different q-point meshes and compare the results using py4vasp:
import py4vasp
from plotly.subplots import make_subplots
from plotly.express.colors import sample_colorscale
calculations = [10,20,30,40]
c = sample_colorscale('viridis',[.8/(len(calculations)-1)*i for i in range(len(calculations))])
#create plotly figure
final_fig = make_subplots()
graph_list = []
for i,label in enumerate(calculations):
mycalc = py4vasp.Calculation.from_file(f"./e06_mgo-dispersion/supercell_2x2x2/dos/vaspout_{label}.h5")
graph = mycalc.phonon_dos.plot()
graph.series[0].name=label
graph_list.append(graph)
if i == 0:
final_graph = graph
else:
final_graph.series = final_graph.series + graph.series
# final_graph
from py4vasp import plot
data = list((g.series[0].x, g.series[0].y, g.series[0].name) for g in graph_list)
plot(*data)
What conclusion can you take from the plots above?
Click to see an answer!
We computed the phonon density of states (DOS) by interpolating the phonon frequencies on increasingly dense q-point grids. Beyond a certain density of q points, the phonon DOS is well described, and increasing it further does not change the DOS. There are still visible differences between the DOS computed on a q-point mesh of 10×10×10 and 20×20×20. In contrast, the differences between the 30×30×30 and 40×40×40 q-point mesh are negligible. In either case, the interpolation of the phonon modes is rather quick with the computation on 40×40×40 q points taking 5 seconds with force constants from a 2x2x2 supercell.
Now, look at the atom-decomposed phonon DOS using py4vasp!
import py4vasp
mycalc = py4vasp.Calculation.from_file("./e06_mgo-dispersion/supercell_2x2x2/dos/vaspout_40.h5")
mycalc.phonon_dos.plot(selection="Mg O")
Looking at the atom decomposed phonon density of states, we see that the Mg atom has the most weight in the phonon modes at low phonon frequencies while the O atom has most of the weight at high phonon frequencies.
What do you think is the reason for this?
The O atom has a lower mass than the Mg atom which yields higher phonon frequencies in the modes in which it is oscillating. The other factor determining the magnitude of the phonon frequencies are the force constants whose strength depends on the strength of the chemical bonds being deformed by the phonon vibration mode.
6.4 Questions¶
- The phonon modes are computed exactly at q points commensurate with the supercell size. What type of interpolation is used to obtain the phonon modes and frequencies for points not commensurate with the supercell size?
- What are the two main quantities determining the magnitude of the phonon frequency for a given phonon mode?
7.1 Task¶
Compute the static dielectric tensor and Born effective charges for MgO.
For semiconductors or insulators, the electronic screening of the ions is incomplete which leads to long-range (LR) interatomic force constants. To compute them explicitly would require infinitely large supercell calculations. For practical calculations, a finite size truncation is needed which leads to Gibbs oscillations in the phonon dispersion. Fortunately, this long-range behavior can be modeled by looking at the analytic form of the ion-ion contribution to the total energy.
For that, follow the approach outlined in Ref. 1 and start by splitting the second-order force constants into short-range and long-range parts, $$ \Phi_{I\alpha J\beta} = \Phi_{I\alpha J\beta}^\text{SR} + \Phi_{I\alpha J\beta}^\text{LR} $$ with the long-range part being obtained from the analytic derivative of the long-range part of the ion-ion contribution to the total energy $E_\text{ion-ion}$. This contribution is typically evaluated using an Ewald-sum technique, in which we separate the ion-ion contribution to the total energy into two parts: the real-space term captures the short-range interaction and reciprocal-space term captures the long-range part $E_\text{ion-ion} = E^\text{SR}_\text{ion-ion} + E^\text{LR}_\text{ion-ion}$. An Ewald parameter $\lambda$ governs the separation and represents a truncation length.
This leads to the following analytical expression for the long-range interatomic force constants, $$ \Phi_{I\alpha J\beta}^\text{LR} = \frac{4\pi e^2}{\Omega_0} \sum_\mathbf{G} \frac{(\mathbf{G} \cdot \mathbf{Z}^*_{I\alpha})(\mathbf{G} \cdot \mathbf{Z}^*_{J\beta})} {\mathbf{G} \cdot \epsilon^\infty \cdot \mathbf{G}} e^{i\mathbf{G} \cdot (\mathbf{R}_J-\mathbf{R}_I)} \exp\Bigl[\frac{-\mathbf{G} \cdot \epsilon^\infty \cdot \mathbf{G}}{4\lambda^2}\Bigr] $$
with $\epsilon^\infty$ the clamped ion dielectric tensor, $\mathbf{Z}^*_{I\alpha}$ the Born effective charges, $\alpha$ the Ewald parameter which is chosen such that the contributions from $\exp[-\mathbf{G} \cdot \epsilon^\infty \cdot \mathbf{G} / (4\lambda^2)]$ are negligible within a certain $\mathbf{G}$ vector cutoff sphere PHON_G_CUTOFF.
- X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997).
MgO
1.0
0.000000000000000 2.122920716207933 2.122920716207933
2.122920716207933 0.000000000000000 2.122920716207933
2.122920716207933 2.122920716207933 0.000000000000000
Mg O
1 1
Direct
0.0 0.0 0.0
0.5 0.5 0.5
PREC = Accurate
EDIFF = 1.E-08
ISMEAR = 0
SIGMA = 0.05
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.
# Compute dielectric tensor and Born effective charges
LEPSILON = .TRUE.
Automatically generated mesh
0
Gamma
6 6 6
0 0 0
PAW_PBE Mg 13Apr2007
PAW_PBE O 08Apr2002
7.3 Calculation¶
Run the VASP calculation using the following command
cd $TUTORIALS/phonon/e07_mgo_dielectric_becs
mpirun -np 4 vasp_std
You can now read the dielectric tensor and born effective charges from the vaspout.h5 file using py4vasp.
import py4vasp
calc = py4vasp.Calculation.from_path("./e07_mgo_dielectric_becs")
calc.born_effective_charge
calc.dielectric_tensor
This information now needs to be present in the run where we compute the phonon dispersion. To do so we need to specify the PHON_DIELECTRIC and PHON_BORN_CHARGES tags which contain the dielectric tensor and Born effective charges respectively. These tags will be used in the INCAR file in the next section to perform the special treatment of the long-range force constants.
The Python code below reads the information from the vaspout.h5 file and creates the PHON_DIELECTRIC and PHON_BORN_CHARGES tags.
def formatVal(v):
return f"{v:18.12f}"
bornCharges = calc.born_effective_charge.to_dict()["charge_tensors"]
dielectricTensor = calc.dielectric_tensor.to_dict()["clamped_ion"]
print("PHON_DIELECTRIC = \\")
strings = []
for i in range(3):
strings.append(" ".join(map(formatVal, dielectricTensor[:,i])))
print("\\\n".join(strings))
print("PHON_BORN_CHARGES = \\")
strings = []
for b in bornCharges:
for i in range(3):
strings.append(" ".join(map(formatVal, b[:,i])))
print("\\\n".join(strings))
7.4 Questions¶
- What is the origin of the long-range interaction of the force constants?
- What quantities are required to describe these long-range force constants?
8.1 Task¶
Compute the long-range part of the interatomic force constants for polar materials using the static dielectric tensor and Born effective charges for MgO.
The separation of the force constants outlined in the previous section allows us to separate the dynamical matrix into short- and long-range parts $$ D_{I\alpha J\beta} (\mathbf{q}) = D^\text{SR}_{I\alpha J\beta} (\mathbf{q}) + D^\text{LR}_{I\alpha J\beta} (\mathbf{q}), $$ with the long-range part of the dynamical matrix
$$ D^\text{LR}_{i\alpha j\beta} (\mathbf{q}) = \frac{4\pi e^2}{\Omega_0} \sum_\mathbf{G} \sum_{l} \frac{\big[ (\mathbf{G}+\mathbf{q}) \cdot \mathbf{Z}^*_{i\alpha} \big] \big[ (\mathbf{G}+\mathbf{q}) \cdot \mathbf{Z}^*_{j\beta} \big]} {(\mathbf{G}+\mathbf{q}) \cdot \epsilon^\infty \cdot (\mathbf{G}+\mathbf{q})} e^{i(\mathbf{q}+\mathbf{G}) \cdot (\mathbf{l}+\mathbf{r}_i-\mathbf{r}_j)} e^\frac{-(\mathbf{G}+\mathbf{q}) \cdot \epsilon^\infty \cdot (\mathbf{G}+\mathbf{q})}{4\lambda^2}. $$
The equations above give us the practical method for computing the phonon dynamical matrices including the long-range force constants using a moderately sized supercell calculation with the steps:
- Compute $\Phi_{I\alpha J\beta}$ using a finite size supercell
- Compute $\Phi_{I\alpha J\beta}^\text{SR}=\Phi_{I\alpha J\beta}-\Phi^\text{LR}_{I\alpha J\beta}$
- Compute $D^\text{SR}_{i\alpha i\beta}(\mathbf{q})$ using $\Phi_{I\alpha J\beta}^\text{SR}$ (Fourier interpolation)
- Compute $D^\text{LR}_{i\alpha i\beta}(\mathbf{q})$ in the unit cell and add to $D^\text{SR}_{i\alpha i\beta}(\mathbf{q})$
The treatment is done automatically inside VASP using the LPHON_POLAR=.TRUE. tag and specifying the dielectric tensor with PHON_DIELECTRIC and the Born effective charges with PHON_BORN_CHARGES.
Plot the phonon dispersion and the density of states now including the special treatment of the long-range dipole-dipole interaction. Compare with the plots generated in section 6
8.2 Input¶
Copy the vasoput.h5 to vaspin.h5 using the following commands
cp $TUTORIALS/phonon/e05_*/supercell_2x2x2/vaspout.h5 $TUTORIALS/phonon/e08_mgo-dispersion-polar/supercell_2x2x2/dispersion/vaspin.h5
cp $TUTORIALS/phonon/e05_*/supercell_2x2x2/vaspout.h5 $TUTORIALS/phonon/e08_mgo-dispersion-polar/supercell_2x2x2/dos/vaspin.h5
8.3 Calculation¶
Run VASP using the following command
cd $TUTORIALS/phonon/e08_mgo-dispersion-polar/supercell_2x2x2/dispersion
vasp_std
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e08_mgo-dispersion-polar/supercell_2x2x2/dispersion")
mycalc.phonon_band.plot(selection="Mg O")
Compute the phonon DOS for different q-point meshes (same as in section 6) using the following command
cd $TUTORIALS/phonon/e08_mgo-dispersion-polar/supercell_2x2x2/dos
bash run.sh
and look at the DOS for the densest q-point mesh
8.4 Questions¶
- What is the main change in the phonon dispersion when including the special treatment of the long-range force constants?
As a bonus task repeat the calculations of sections 6 and 8 starting from force constants computed in the 3×3×3 and 4×4×4 supercells and compare the results. This is instructive to verify how the special treatment of the long-range dipole-dipole interaction leads to a smoother interpolation of the phonon dispersion.