Part 1: Graphene¶
By the end of this tutorial, you will be able to:
- list some limitations of the static-lattice approximation
- perform a relaxation of the ionic position and cell shape at fixed volume for 2D materials
- set the verbosity of the OUTCAR file
- visualize the crystal structure using py4vasp
1.1 Task¶
Relax the ionic positions and cell shape of graphene at fixed volume.
The standard ab-initio calculation uses a static lattice model, where the ions constitute a fixed, rigid, immobile periodic structure. That is, the electronic and ionic degrees of freedom are separated using the Born-Oppenheimer approximation. To find a stable structure, we typically perform the following sequence of steps
- chooses an initial structure
- treats the electrons quantum mechanically
- updates the ionic positions based on the ab-initio forces (Hellmann—Feyman theorem).
The static lattice model describes many phenomena, but there are considerable limitations when lattice dynamics play a dominant role. For instance, phonons contribute more to the specific heat than electrons for many materials. Furthermore, the static lattice model cannot describe thermal expansion, thermal conductivity, piezoelectricity, melting, the transmission of sound, the temperature-dependent intensity of X-ray diffraction, certain aspects of neutron scattering, and some optical and dielectric properties. An example for the latter is crystals that show resonances in their reflectivity at energies far below the bandgap.
1.2 Input¶
The input files to run this example are prepared at $TUTORIALS/phonon/e01_static-lattice. Check them out!
Graphene
2.458
1.0000000000000000 0.0000000000000000 0.0000000000000000
-0.5000000000000000 0.8660254037844386 0.0000000000000000
-0.0000000000000000 -0.0000000000000000 3.2546786004882016
C
2
Selective dynamics
Direct
0.3333333333333333 0.6666666666666666 0.0000000000000000 T T F
0.6666666666666666 0.3333333333333333 0.0000000000000000 T T F
SYSTEM = graphene
NWRITE = 3
ENCUT = 400
# electronic
PREC = Accurate
EDIFF = 1e-8
ISMEAR = -1
SIGMA = 0.2
# ionic
IBRION = 2
NSW = 100
ISIF = 4
k mesh
0
Gamma
12 12 1
0 0 0
PAW C_s 04May1998
The POSCAR file defines the graphene structure; a single layer of carbon atoms in a Honeycomb lattice. An 8-Å vacuum layer suppresses the interaction between periodic replica in the z direction. The selective dynamics blocks relaxations in z direction and does not impact the resulting graphene after the relaxation.
The INCAR file specifies tags regarding the electronic and ionic minimization. The electronic and ionic degrees of freedom are updated iteratively:
- For a given structure, the Kohn-Sham equations define the electronic ground state.
- The gradient of the electronic energy (Hellmann—Feyman theorem) yields the forces and stress acting on the structure.
- The IBRION tag determines the optimization algorithm for the structure (here: conjugate gradient) resulting in a new structure.
- VASP repeats step 1 to 3 until satisfying the convergence criterion (EDIFFG).
Note that EDIFF must be small enough to avoid a spurious convergence. If EDIFF is too large, the electronic density may not be converged well enough to obtain precise forces and stress. As a result, VASP may struggle to converge the structure or potentially end up with an incorrect solution. Therefore, it is a good idea to check whether the forces and stress really vanish using a tighter EDIFF cutoff.
What does ISIF=4 imply? Do we need to worry about Pulay stress in this calculation? Read the VASP Wiki to find the answer!
Click to see the answer!
Even though the volume is fixed, the lattice parameters are allowed to change. This may induce Pulay stress, because the number of G vectors satisfying the energy cutoff ENCUT may change. Take care of Pulay stress by restarting after a couple of ionic steps especially if large changes to the unit cell shape occurred. Alternatively, increase ENCUT to reduce the impact of the Pulay stress.
The verbosity of the OUTCAR file is controlled by the NWRITE tag. Here, we expect a rather short calculation including both electronic and ionic steps. Therefore, it is a good opportunity to increase verbosity and observe the additional output to the OUTCAR file.
The KPOINTS file defines a Γ-centered equidistant $\mathbf{k}$ mesh. Monolayers do not need a subdivision of the $\mathbf{k}$ mesh normal to the surface direction. Additional $\mathbf{k}$ points in the direction of the surface normal would merely describe the interaction between the periodic images of the graphene layer.
1.3 Calculation¶
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/phonon/e01_*
mpirun -np 2 vasp_std
Visualize the relaxed graphene structure using py4vasp!
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e01_static-lattice")
mycalc.structure.plot()
To see the typical graphene layers, plot a 3x3x1 supercell using py4vasp!
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e01_static-lattice")
mycalc.structure.plot(supercell=(3,3,1))
Open the OUTCAR file and find information about the forces computed at each ionic iteration! Can you identify additional information due to the NWRITE tag? What is the force on each ionic site in the first and in the final iteration step? What is the convergence criterion? What is the value of the EDIFFG tag in this calculation?
In the next example, we continue with the structure in the POSCAR file because it is symmetrized.
By the end of this tutorial, you will be able to:
- create supercells for phonon calculations
- state the definition and compute the force constants using finite differences
- judge when to set LREAL=Auto
- explain the use of symmetry in the context of the finite differences algorithm
2.1 Task¶
Compute the force constants for graphene
We consider the Born-Oppenheimer energy surface $E(\{\mathbf{R}\})$. This ground-state energy of the electronic system depends on the position $\mathbf{R}_{I} = \mathbf{R}^0_{I} + \mathbf{u}_{I}$. Here, $\mathbf{u}_{I}$ is a the displacement of ion $I$ with respect to the ground state. Expanding the total energy $E$ around the unperturbed system with total energy $E^0$ yields $$\tag{1.1} E(\{\mathbf{R}\})= E(\{\mathbf{R}^0\}) - \sum_{I\alpha} F_{I\alpha} (\{\mathbf{R}^0\}) u_{I\alpha}+ \sum_{I\alpha J\beta} \Phi_{I\alpha J\beta} (\{\mathbf{R}^0\}) u_{I\alpha} u_{J\beta} + \mathcal{O}(\mathbf{R}^3), $$ where $\alpha$ and $\beta$ run over the Cartesian axis. The harmonic approximation truncates the expansion at the second derivative of the total energy with ionic displacements. This approximation is reasonable when the displacements are small with respect to the equilibrium position.
The name and sign convention of the expansion coefficients relate to their physical interpretation. The first-order coefficient is the force $F$ and the second-order one the Hessian or force-constant matrix $\Phi$: $$\tag{1.2} F_{I\alpha} (\{\mathbf{R}^0\}) = - \left. \frac{\partial E(\{\mathbf{R}\})}{\partial R_{I\alpha}} \right|_{\mathbf{R} =\mathbf{R^0}} ,\quad \Phi_{I\alpha J\beta} (\{\mathbf{R}^0\}) = \left. \frac{\partial E(\{\mathbf{R}\})}{\partial R_{I\alpha} \partial R_{J\beta}} \right|_{\mathbf{R} =\mathbf{R^0}} = - \left. \frac{\partial F_{I\alpha}(\{\mathbf{R}\})}{\partial R_{J\beta}} \right|_{\mathbf{R} =\mathbf{R^0}} . $$ The forces are computed using the Hellmann–Feynman theorem.
To capture phonon modes with wavelengths larger than the primitive unit cell, VASP requires a calculation of the appropriate supercell. Mind: Converge all phonon calculations with respect to the size of the supercell. When increasing the supercell size by a factor, reduce the $\mathbf{k}$-point mesh by the same factor so that the sampling density remains invariant.
2.2 Input¶
The input files to run this example are prepared at $TUTORIALS/phonon/e02_fin-diff-force-constants.
The POSCAR file contains a previously relaxed primitive cell of graphene.
In the following python code, we generate 3x3, 4x4, 5x5, and 6x6 supercells of graphene replicated in-plane.
Run them and inspect the different generated POSCAR files.
from pymatgen.core import Structure
from pymatgen.io.vasp import Poscar
my_struc = Structure.from_file("./e02_fin-diff-force-constants/POSCAR")
# make a 3x3x1 supercell
my_struc.make_supercell([3,3,1])
# write supercell to POSCAR format with specified filename
Poscar(my_struc).write_file(filename="./e02_fin-diff-force-constants/supercell_3x3x1/POSCAR",significant_figures=16)
# A POSCAR file will appear at $TUTORIALS/MD/e02_fin-diff-force-constants/supercell_3x3x1/
# You may need to refresh the file browser to see it
my_struc = Structure.from_file("./e02_fin-diff-force-constants/POSCAR")
my_struc.make_supercell([4,4,1])
Poscar(my_struc).write_file(filename="./e02_fin-diff-force-constants/supercell_4x4x1/POSCAR",significant_figures=16)
my_struc = Structure.from_file("./e02_fin-diff-force-constants/POSCAR")
my_struc.make_supercell([5,5,1])
Poscar(my_struc).write_file(filename="./e02_fin-diff-force-constants/supercell_5x5x1/POSCAR",significant_figures=16)
my_struc = Structure.from_file("./e02_fin-diff-force-constants/POSCAR")
my_struc.make_supercell([6,6,1])
Poscar(my_struc).write_file(filename="./e02_fin-diff-force-constants/supercell_6x6x1/POSCAR",significant_figures=16)
Click to see supercell_3x3x1/POSCAR!
C18
1.0
7.3521657209830806 0.0000000000000000 0.0000000000000000
-3.6760828604915403 6.3671622872044793 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.0000000000000000
C
18
direct
0.1111111111111133 0.2222222222222200 0.0000000000000000 C
0.1111111111111133 0.5555555555555532 0.0000000000000000 C
0.1111111111111133 0.8888888888888866 0.0000000000000000 C
0.4444444444444466 0.2222222222222200 0.0000000000000000 C
0.4444444444444466 0.5555555555555532 0.0000000000000000 C
0.4444444444444466 0.8888888888888866 0.0000000000000000 C
0.7777777777777801 0.2222222222222200 0.0000000000000000 C
0.7777777777777799 0.5555555555555532 0.0000000000000000 C
0.7777777777777799 0.8888888888888866 0.0000000000000000 C
0.2222222222222200 0.1111111111111133 0.0000000000000000 C
0.2222222222222200 0.4444444444444466 0.0000000000000000 C
0.2222222222222199 0.7777777777777799 0.0000000000000000 C
0.5555555555555532 0.1111111111111133 0.0000000000000000 C
0.5555555555555534 0.4444444444444466 0.0000000000000000 C
0.5555555555555534 0.7777777777777799 0.0000000000000000 C
0.8888888888888866 0.1111111111111133 0.0000000000000000 C
0.8888888888888866 0.4444444444444466 0.0000000000000000 C
0.8888888888888866 0.7777777777777799 0.0000000000000000 C
Click to see supercell_4x4x1/POSCAR!
C32
1.0
9.8028876279774408 0.0000000000000000 0.0000000000000000
-4.9014438139887204 8.4895497162726397 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.0000000000000000
C
32
direct
0.0833333333333350 0.1666666666666650 0.0000000000000000 C
0.0833333333333350 0.4166666666666650 0.0000000000000000 C
0.0833333333333350 0.6666666666666650 0.0000000000000000 C
0.0833333333333349 0.9166666666666649 0.0000000000000000 C
0.3333333333333350 0.1666666666666650 0.0000000000000000 C
0.3333333333333350 0.4166666666666650 0.0000000000000000 C
0.3333333333333350 0.6666666666666650 0.0000000000000000 C
0.3333333333333349 0.9166666666666649 0.0000000000000000 C
0.5833333333333350 0.1666666666666650 0.0000000000000000 C
0.5833333333333350 0.4166666666666650 0.0000000000000000 C
0.5833333333333350 0.6666666666666650 0.0000000000000000 C
0.5833333333333349 0.9166666666666649 0.0000000000000000 C
0.8333333333333350 0.1666666666666650 0.0000000000000000 C
0.8333333333333350 0.4166666666666650 0.0000000000000000 C
0.8333333333333350 0.6666666666666650 0.0000000000000000 C
0.8333333333333349 0.9166666666666649 0.0000000000000000 C
0.1666666666666650 0.0833333333333350 0.0000000000000000 C
0.1666666666666650 0.3333333333333350 0.0000000000000000 C
0.1666666666666650 0.5833333333333349 0.0000000000000000 C
0.1666666666666649 0.8333333333333348 0.0000000000000000 C
0.4166666666666650 0.0833333333333350 0.0000000000000000 C
0.4166666666666650 0.3333333333333350 0.0000000000000000 C
0.4166666666666650 0.5833333333333349 0.0000000000000000 C
0.4166666666666649 0.8333333333333348 0.0000000000000000 C
0.6666666666666650 0.0833333333333350 0.0000000000000000 C
0.6666666666666650 0.3333333333333350 0.0000000000000000 C
0.6666666666666650 0.5833333333333349 0.0000000000000000 C
0.6666666666666649 0.8333333333333348 0.0000000000000000 C
0.9166666666666650 0.0833333333333350 0.0000000000000000 C
0.9166666666666650 0.3333333333333350 0.0000000000000000 C
0.9166666666666650 0.5833333333333349 0.0000000000000000 C
0.9166666666666649 0.8333333333333348 0.0000000000000000 C
Click to see supercell_5x5x1/POSCAR!
C50
1.0
12.2536095349718011 0.0000000000000000 0.0000000000000000
-6.1268047674859005 10.6119371453408000 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.0000000000000000
C
50
direct
0.0666666666666680 0.1333333333333320 0.0000000000000000 C
0.0666666666666680 0.3333333333333320 0.0000000000000000 C
0.0666666666666680 0.5333333333333320 0.0000000000000000 C
0.0666666666666680 0.7333333333333321 0.0000000000000000 C
0.0666666666666680 0.9333333333333320 0.0000000000000000 C
0.2666666666666680 0.1333333333333320 0.0000000000000000 C
0.2666666666666680 0.3333333333333320 0.0000000000000000 C
0.2666666666666680 0.5333333333333320 0.0000000000000000 C
0.2666666666666680 0.7333333333333321 0.0000000000000000 C
0.2666666666666679 0.9333333333333320 0.0000000000000000 C
0.4666666666666680 0.1333333333333320 0.0000000000000000 C
0.4666666666666680 0.3333333333333320 0.0000000000000000 C
0.4666666666666680 0.5333333333333320 0.0000000000000000 C
0.4666666666666680 0.7333333333333321 0.0000000000000000 C
0.4666666666666680 0.9333333333333320 0.0000000000000000 C
0.6666666666666681 0.1333333333333320 0.0000000000000000 C
0.6666666666666681 0.3333333333333320 0.0000000000000000 C
0.6666666666666681 0.5333333333333320 0.0000000000000000 C
0.6666666666666681 0.7333333333333321 0.0000000000000000 C
0.6666666666666681 0.9333333333333320 0.0000000000000000 C
0.8666666666666680 0.1333333333333320 0.0000000000000000 C
0.8666666666666680 0.3333333333333320 0.0000000000000000 C
0.8666666666666680 0.5333333333333320 0.0000000000000000 C
0.8666666666666680 0.7333333333333321 0.0000000000000000 C
0.8666666666666680 0.9333333333333320 0.0000000000000000 C
0.1333333333333320 0.0666666666666680 0.0000000000000000 C
0.1333333333333320 0.2666666666666680 0.0000000000000000 C
0.1333333333333320 0.4666666666666680 0.0000000000000000 C
0.1333333333333320 0.6666666666666681 0.0000000000000000 C
0.1333333333333319 0.8666666666666679 0.0000000000000000 C
0.3333333333333320 0.0666666666666680 0.0000000000000000 C
0.3333333333333319 0.2666666666666680 0.0000000000000000 C
0.3333333333333319 0.4666666666666680 0.0000000000000000 C
0.3333333333333320 0.6666666666666681 0.0000000000000000 C
0.3333333333333319 0.8666666666666679 0.0000000000000000 C
0.5333333333333320 0.0666666666666680 0.0000000000000000 C
0.5333333333333320 0.2666666666666680 0.0000000000000000 C
0.5333333333333320 0.4666666666666680 0.0000000000000000 C
0.5333333333333320 0.6666666666666681 0.0000000000000000 C
0.5333333333333319 0.8666666666666679 0.0000000000000000 C
0.7333333333333321 0.0666666666666680 0.0000000000000000 C
0.7333333333333321 0.2666666666666680 0.0000000000000000 C
0.7333333333333321 0.4666666666666680 0.0000000000000000 C
0.7333333333333321 0.6666666666666681 0.0000000000000000 C
0.7333333333333321 0.8666666666666679 0.0000000000000000 C
0.9333333333333319 0.0666666666666680 0.0000000000000000 C
0.9333333333333319 0.2666666666666680 0.0000000000000000 C
0.9333333333333319 0.4666666666666680 0.0000000000000000 C
0.9333333333333319 0.6666666666666681 0.0000000000000000 C
0.9333333333333319 0.8666666666666679 0.0000000000000000 C
Click to see supercell_6x6x1/POSCAR!
C72
1.0
14.7043314419661613 0.0000000000000000 0.0000000000000000
-7.3521657209830806 12.7343245744089586 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.0000000000000000
C
72
direct
0.0555555555555567 0.1111111111111100 0.0000000000000000 C
0.0555555555555566 0.2777777777777766 0.0000000000000000 C
0.0555555555555566 0.4444444444444433 0.0000000000000000 C
0.0555555555555567 0.6111111111111101 0.0000000000000000 C
0.0555555555555566 0.7777777777777765 0.0000000000000000 C
0.0555555555555567 0.9444444444444433 0.0000000000000000 C
0.2222222222222233 0.1111111111111100 0.0000000000000000 C
0.2222222222222233 0.2777777777777766 0.0000000000000000 C
0.2222222222222233 0.4444444444444433 0.0000000000000000 C
0.2222222222222233 0.6111111111111101 0.0000000000000000 C
0.2222222222222232 0.7777777777777765 0.0000000000000000 C
0.2222222222222233 0.9444444444444433 0.0000000000000000 C
0.3888888888888901 0.1111111111111100 0.0000000000000000 C
0.3888888888888899 0.2777777777777766 0.0000000000000000 C
0.3888888888888899 0.4444444444444433 0.0000000000000000 C
0.3888888888888900 0.6111111111111101 0.0000000000000000 C
0.3888888888888899 0.7777777777777765 0.0000000000000000 C
0.3888888888888900 0.9444444444444433 0.0000000000000000 C
0.5555555555555567 0.1111111111111100 0.0000000000000000 C
0.5555555555555567 0.2777777777777766 0.0000000000000000 C
0.5555555555555567 0.4444444444444433 0.0000000000000000 C
0.5555555555555567 0.6111111111111101 0.0000000000000000 C
0.5555555555555566 0.7777777777777765 0.0000000000000000 C
0.5555555555555567 0.9444444444444433 0.0000000000000000 C
0.7222222222222234 0.1111111111111100 0.0000000000000000 C
0.7222222222222233 0.2777777777777766 0.0000000000000000 C
0.7222222222222233 0.4444444444444433 0.0000000000000000 C
0.7222222222222234 0.6111111111111101 0.0000000000000000 C
0.7222222222222232 0.7777777777777765 0.0000000000000000 C
0.7222222222222234 0.9444444444444433 0.0000000000000000 C
0.8888888888888899 0.1111111111111100 0.0000000000000000 C
0.8888888888888898 0.2777777777777766 0.0000000000000000 C
0.8888888888888899 0.4444444444444433 0.0000000000000000 C
0.8888888888888899 0.6111111111111101 0.0000000000000000 C
0.8888888888888897 0.7777777777777765 0.0000000000000000 C
0.8888888888888899 0.9444444444444433 0.0000000000000000 C
0.1111111111111100 0.0555555555555567 0.0000000000000000 C
0.1111111111111100 0.2222222222222233 0.0000000000000000 C
0.1111111111111100 0.3888888888888899 0.0000000000000000 C
0.1111111111111100 0.5555555555555567 0.0000000000000000 C
0.1111111111111099 0.7222222222222232 0.0000000000000000 C
0.1111111111111100 0.8888888888888899 0.0000000000000000 C
0.2777777777777766 0.0555555555555567 0.0000000000000000 C
0.2777777777777767 0.2222222222222233 0.0000000000000000 C
0.2777777777777767 0.3888888888888899 0.0000000000000000 C
0.2777777777777767 0.5555555555555567 0.0000000000000000 C
0.2777777777777766 0.7222222222222232 0.0000000000000000 C
0.2777777777777767 0.8888888888888899 0.0000000000000000 C
0.4444444444444433 0.0555555555555567 0.0000000000000000 C
0.4444444444444433 0.2222222222222233 0.0000000000000000 C
0.4444444444444433 0.3888888888888899 0.0000000000000000 C
0.4444444444444433 0.5555555555555567 0.0000000000000000 C
0.4444444444444433 0.7222222222222232 0.0000000000000000 C
0.4444444444444433 0.8888888888888899 0.0000000000000000 C
0.6111111111111101 0.0555555555555567 0.0000000000000000 C
0.6111111111111101 0.2222222222222233 0.0000000000000000 C
0.6111111111111099 0.3888888888888899 0.0000000000000000 C
0.6111111111111101 0.5555555555555567 0.0000000000000000 C
0.6111111111111099 0.7222222222222232 0.0000000000000000 C
0.6111111111111101 0.8888888888888899 0.0000000000000000 C
0.7777777777777767 0.0555555555555567 0.0000000000000000 C
0.7777777777777766 0.2222222222222233 0.0000000000000000 C
0.7777777777777767 0.3888888888888899 0.0000000000000000 C
0.7777777777777767 0.5555555555555567 0.0000000000000000 C
0.7777777777777766 0.7222222222222232 0.0000000000000000 C
0.7777777777777767 0.8888888888888899 0.0000000000000000 C
0.9444444444444432 0.0555555555555567 0.0000000000000000 C
0.9444444444444433 0.2222222222222233 0.0000000000000000 C
0.9444444444444432 0.3888888888888899 0.0000000000000000 C
0.9444444444444432 0.5555555555555567 0.0000000000000000 C
0.9444444444444431 0.7222222222222232 0.0000000000000000 C
0.9444444444444433 0.8888888888888899 0.0000000000000000 C
SYSTEM = graphene
ENCUT = 400
# electronic
PREC = Accurate
NELMIN = 5
EDIFF = 1e-8
ISMEAR = -1
SIGMA = 0.2
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.
# ionic (finite differences)
IBRION = 6
POTIM = 0.015
Click to see supercell_3x3x1/KPOINTS!
K points
0
Gamma
4 4 1
0 0 0
Click to see supercell_4x4x1/KPOINTS!
K points
0
Gamma
3 3 1
0 0 0
Click to see supercell_5x5x1/KPOINTS!
K points
0
Gamma
2 2 1
0 0 0
Click to see supercell_6x6x1/KPOINTS!
K points
0
Gamma
2 2 1
0 0 0
PAW C_s 04May1998
In the INCAR file, IBRION=6 triggers a computation of the force constants and phonon modes using the finite differences approach.
2.3 Calculation¶
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/phonon/e02_*
for d in 3x3x1 4x4x1 5x5x1 6x6x1
do
cd supercell_$d
mpirun -np 4 vasp_std
cd ..
done
You will see an advice:
-----------------------------------------------------------------------------
| |
| ----> ADVICE to this user running VASP <---- |
| |
| You have a (more or less) 'large supercell' and for larger cells it |
| might be more efficient to use real-space projection operators. |
| Therefore, try LREAL= Auto in the INCAR file. |
| Mind: For very accurate calculation, you might also keep the |
| reciprocal projection scheme (i.e. LREAL=.FALSE.). |
| |
-----------------------------------------------------------------------------
The calculation of the force constants takes some time for the larger 5x5x1 and 6x6x1 supercells (5 and 12 minutes respectively). You can proceed with the tutorial while the larger calculations are running.
In the meantime let us understand the projection operators and the meaning of the LREAL tag. Read the VASP Wiki to answer the following questions: What is aliasing? The evaluation of which part of the pseudo potential requires performing a projection? What is projected onto what? Which tag should be set in conjuction with LREAL=Auto? How is the appropriate setting for that tag determined? Should we use LREAL=Auto or LREAL=False here? Why?
Plot the supercells using py4vasp!
import py4vasp
supercell3 = py4vasp.Calculation.from_path("./e02_fin-diff-force-constants/supercell_3x3x1")
supercell3.structure[:].plot()
py4vasp visualizes the ground state structure along with 4 perturbed ones used for the finite-difference algorithm. When you press iterate over the structure in a loop, recognize that one atom wiggles. Structure 0 is the unperturbed structure. VASP exploits the symmetry of the system, to extract the force constants of the entire system using only these small ionic displacements. Will the larger supercell have more degrees of freedom? Execute the next cell!
supercell4 = py4vasp.Calculation.from_path("./e02_fin-diff-force-constants/supercell_4x4x1")
supercell4.structure[:].plot()
There are also 5 structures: the equilibrium structure and 2 ionic displacements for each of the 2 degrees of freedom. Why are there 2 degrees of freedom and not 3 $\times$ the number of ions in the unit cell?
Click to see a hint!
Symmetry operations map the carbon atoms onto each other. For the larger cell, a simple translation of the crystal maps all atoms into the primitive cell with 2 atoms. But there are also symmetry operations in the primitive cell to consider.
What is the length of the lattice vectors? Which angle do they enclose? What other symmetries might be taken into account?
Click to see the answer!
Graphene has a hexagonal structure and, thus, the in-plane lattice vectors have the same length. They enclose an angle of $60^\circ$ with a 6-fold rotational symmetry around the third lattice vector. This point-group operation reduces the number of degrees of freedom from 6 to 4, because for each site there is only one in-plane degree of freedom. The inversion symmetry reduces the number of degrees of freedom by two, because only one of the sites needs to be displaced.
If you use IBRION=5 instead of 6, the finite difference algorithm will not exploit the symmetry of the system.
supercell5 = py4vasp.Calculation.from_path("./e02_fin-diff-force-constants/supercell_5x5x1")
supercell5.structure[:].plot()
supercell6 = py4vasp.Calculation.from_path("./e02_fin-diff-force-constants/supercell_6x6x1")
supercell6.structure[:].plot()
The force constant matrix can be retrieved using py4vasp. Run the following command! How many entries do you find for the different supercells?
supercell3.force_constant.print()
supercell4.force_constant.print()
VASP computes the force constants in Cartesian coordinates. In the case of graphene and other 2D materials, the in-plane force constants do not couple to the out-of-plane ones. This manifests as vanishing force constants for all off-diagonal elements involving in-plane ($x$ and $y$) and out-of-plane ($z$) components.
Let us plot the out-of-plane force constants as a function of the distance between different nuclei. With the following function, we select one atom index and evaluate the Cartesian distance to all the other atoms. We use the distance to the periodic image whenever that is smaller. Sorting the distances by magnitude allows plotting the corresponding out-of-plane force constants.
import numpy as np
from py4vasp import Calculation, plot
atom_index = 0
def out_of_plane_force_constants(supercell):
calculation = Calculation.from_path(f"./e02_fin-diff-force-constants/supercell_{supercell}")
structure = calculation.structure[0].read()
# evaluate minimal distance for a structure with periodic boundaries
distances = structure["positions"][atom_index] - structure["positions"]
minimal_distances = distances - np.rint(distances)
# convert to absolute distance
cartesian_distances = minimal_distances @ structure["lattice_vectors"]
absolute_distances = np.linalg.norm(cartesian_distances, axis=1)
# sort the distances
order = np.argsort(absolute_distances)
absolute_distances = absolute_distances[order]
# get out-of-plane force constants
slice_ = lambda i: slice(3 * i, 3 * (i + 1))
ii = slice_(atom_index)
force_constants = calculation.force_constant.read()["force_constants"]
force_constants = np.array([force_constants[ii, slice_(j)][2, 2] for j in order])
# skip the atom interacting with itself
return absolute_distances[1:], force_constants[1:], supercell
plot(
out_of_plane_force_constants("3x3x1"),
xlabel='distance (Å)',
ylabel='Force-constant (eV/Ų)',
)
As you can see the force constants oscillate between negative and positive values when considering neighbors at different distances.
What does it mean to have negative force-constants? (Click to see the answer)
The second order force constants are the second derivative of the total energy with respect to two ionic positions. The sign of the force constants determines whether the total energy increases or decreases when these two ions move relative to each other. This begs the question whether negative force constants mean that the structure is meta-stable. Here, meta-stable means a structure with vanishing forces but decreasing total energy upon distortion. However, the appearance of some negative force constants is not sufficient for a meta-stable structure. While the displacement of one pair of ions may lower the total energy, other positive force constants coupling different ion pairs may counterbalance them yielding an overall stable structure.
For the in-plane force constants, we define a local reference frame based on the relative positions of the relevant ion pair. We rotate the 3x3 force constant matrix ($\mathbf{\Phi}$) for the ion pair $(I, J)$ to a new reference the coordinate vectors point along the longitudinal and the transverse direction as well as normal to the monolayer. Here the longitudinal direction $\hat{\mathbf{r}}_{\mathrm L}$ points from one atom to the other and the transverse direction $\hat{\mathbf{r}}_{\mathrm T}$ is perpendicular to it in the plane. The matrix $\mathbf{R} = (\hat{\mathbf{r}}_{\mathrm{L}}, \hat{\mathbf{r}}_{\mathrm T}, \hat{\mathbf{z}})$ and rotates the force-constants matrix $\tilde{\mathbf{\Phi}} = \mathbf{R}^T \mathbf{\Phi} \mathbf{R}$. We plot the longitudinal and transverse components of $\tilde{\mathbf{\Phi}}$ as a function of the distance between the atoms.
import numpy as np
from py4vasp import Calculation, plot
def in_plane_force_constants(calculation, atom_index):
structure = calculation.structure[0].read()
# evaluate minimal distance in period boundaries
distances = structure["positions"][atom_index] - structure["positions"]
minimal_distances = distances - np.rint(distances)
# convert to absolute distance
cartesian_distances = minimal_distances @ structure["lattice_vectors"]
absolute_distances = np.linalg.norm(cartesian_distances, axis=1)
# sort the distances (skip the atom interacting with itself)
order = np.argsort(absolute_distances)[1:]
absolute_distances = absolute_distances[order]
# construct rotation matrix
r_l = (cartesian_distances[order].T / absolute_distances).T
zz = [0,0,1]
r_t = np.cross(r_l, zz)
rotations = np.array([[l, t, zz] for l, t in zip(r_l, r_t)])
# split force constants into longitudinal and transversal
force_constants = calculation.force_constant.read()["force_constants"]
slice_ = lambda i: slice(3 * i, 3 * (i + 1))
ii = slice_(atom_index)
longitudinal = []
transversal = []
for R, j in zip(rotations, order):
Φ = force_constants[ii, slice_(j)]
Φt = R.T @ Φ @ R
longitudinal.append(Φt[0,0])
transversal.append(Φt[1,1])
return (absolute_distances,longitudinal,"longitudinal"), (absolute_distances,transversal,"transversal")
atom_index = 0
my_calc3 = Calculation.from_path("./e02_fin-diff-force-constants/supercell_3x3x1")
plot(
*in_plane_force_constants(my_calc3,atom_index),
xlabel='distance (Å)',
ylabel='force-constant (eV/Ų)',
)
Try plotting the force constants for different supercell sizes.
my_calc4 = Calculation.from_path("./e02_fin-diff-force-constants/supercell_4x4x1")
plot(
*in_plane_force_constants(my_calc4,atom_index),
xlabel='distance (Å)',
ylabel='force-constant (eV/Ų)',
)
my_calc5 = Calculation.from_path("./e02_fin-diff-force-constants/supercell_5x5x1")
plot(
*in_plane_force_constants(my_calc5,atom_index),
xlabel='distance (Å)',
ylabel='force-constant (eV/Ų)',
)
What additional information do you obtain when increasing the supercell size?
When increasing the supercell size the force constants for a larger number of atom pairs are computed.
We can increase the accuracy of the phonon dispersion by including force-constants on increasingly large supercells (i.e considering more distant interacting ions). Because the magnitude of these force-constants decreases for pairs of atoms that are further apart we know that the phonon dispersion will eventually converge with reasonably sized supercells. There are a few notable exceptions to this that we should point out:
- When ions with different charges move, they form dipoles which interact with a slowly decaying force constants (we will tackle this topic in the next tutorial).
- For metals and semi-metals the nuclei movement leads to changes in the Fermi level which leads to changes total energy. This interaction is often long-ranged and is known as a Kohn anomally.
2.4 Questions¶
- What is the difference between IBRION=5 and 6?
- What is the definition of the force constants?
- What does the LREAL tag control?
- Why is the size of the supercell important for a precise phonon dispersion calculation?
By the end of this tutorial, you will be able to:
- compute and plot the phonon dispersion
3.1 Task¶
Compute the phonon dispersion of graphene.
In the harmonic approximation, the ionic Hamiltonian yields $$ H = \frac{1}{2} \sum_{I\alpha} M_I \dot{u}^2_{I\alpha} + \frac{1}{2} \sum_{I\alpha J\beta} \Phi_{I\alpha J\beta} u_{I\alpha} u_{J\beta}, $$ which leads to the following equation of motion: $$ M_I \ddot{u}_{I\alpha} = - \Phi_{I\alpha J\beta} u_{J\beta} $$
When we consider solutions in the form of plane waves traveling parallel to the wave vector $\mathbf{q}$, i.e., $$ \mathbf{u}^\mu_{I\alpha}(\mathbf{q},t) = \frac{1}{2} \frac{1}{\sqrt{M_I}} \left\{ A^\mu(\mathbf{q}) \varepsilon^{\mu }_{I\alpha}(\mathbf{q}) e^{ i [\mathbf{q} \cdot \mathbf{\mathbf{R}}_I -\omega_\mu(\mathbf{q})t]}+ A^{\mu*}(\mathbf{q}) \varepsilon^{\mu*}_{I\alpha}(\mathbf{q}) e^{-i [\mathbf{q} \cdot \mathbf{\mathbf{R}}_I -\omega_\mu(\mathbf{q})t]} \right\} $$ determining the phonon modes $\varepsilon^{\mu }_{I\alpha}(\mathbf{q})$ and vibrational frequencies $\omega^\mu(\mathbf{q})$ amounts to solving: $$\tag{2} \sum_{J\beta} \frac{1}{\sqrt{M_I M_J}} \Phi_{I\alpha J\beta} e^{i\mathbf{q} \cdot (\mathbf{R}_J-\mathbf{R}_I)} \varepsilon^{\mu }_{J\beta}(\mathbf{q}) = \omega^\mu(\mathbf{q})^2 \varepsilon^{\mu }_{I\alpha}(\mathbf{q}) $$ where $\Phi_{I\alpha J\beta}$ is the Hessian with respect to the ionic site $I, J$ with mass $M_{I/J}$ and direction $\alpha,\beta=x, y, z$. $A(\mathbf q)$ is the phonon amplitude and $\varepsilon^{\mu }_{I\alpha}(\mathbf q)$ is the phonon normal modes of vibration. With $N$ ions in the supercell, the dimension of the dynamical matrix $\frac{1}{\sqrt{M_I M_J}} \Phi_{I\alpha J\beta} e^{i\mathbf{q} \cdot (\mathbf{R}_J-\mathbf{R}_I)}$ is thus $3N$, which yields a set of $3N$ linear homogeneous equations for each wave vector $\mathbf q$.
Solving the eigenvalue problem yields eigenvectors $\varepsilon^{\mu }_{I\alpha}(\mathbf q)$ and eigenvalues $\omega^\mu(\mathbf{q})$.
We can write the positions of the atoms in the supercell $\mathbf{R}_I$ in terms of integer multiples of the lattice vectors of the unit cell $\mathbf{L}_i$ such that $\mathbf{R}_I = \mathbf{L}_i + \mathbf{r}_i$ with $\mathbf{r}_i$ being the position of the ion in the unit cell. This allows us to compute the phonons in the unit cell using the following equation $$ \sum_{L_i,L_j} \sum_{i\beta} \frac{1}{\sqrt{M_i M_j}} \Phi_{I\alpha J\beta} e^{i\mathbf{q} \cdot (\mathbf{L}_j-\mathbf{L}_i+\mathbf{r}_j-\mathbf{r}_i)} \varepsilon^{\mu }_{j\beta}(\mathbf{q}) = \omega^\mu(\mathbf{q})^2 \varepsilon^{\mu }_{i\alpha}(\mathbf{q}) $$ in this case, the dynamical matrix has dimension $3n$ with $n$ the number of atoms in the unit cell.
The phonon frequencies in the primitive cell are commensurate for $\mathbf{q}$ points that satisfy $$ \mathbf q = \frac{s_1}{n_1} \mathbf{b}_1 + \frac{s_2}{n_2} \mathbf{b}_2 + \frac{s_3}{n_3} \mathbf{b}_3 $$ where $s_i = 1...n_i$, with $n_i$ the number of repetitions of the unit cell in each direction, $\mathbf q$ is restricted to the first Brillouin zone in reciprocal space spanned by $\mathbf{b}_1$, $\mathbf{b}_2$, and $\mathbf{b}_3$.
3.2 Input¶
The input files to run this example are prepared at $TUTORIALS/phonon/e03_dispersion. Check them out!
Click to see supercell_3x3x1/POSCAR!
C18
1.0
7.3521657209830806 0.0000000000000000 0.0000000000000000
-3.6760828604915403 6.3671622872044793 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.0000000000000000
C
18
direct
0.1111111111111133 0.2222222222222200 0.0000000000000000 C
0.1111111111111133 0.5555555555555532 0.0000000000000000 C
0.1111111111111133 0.8888888888888866 0.0000000000000000 C
0.4444444444444466 0.2222222222222200 0.0000000000000000 C
0.4444444444444466 0.5555555555555532 0.0000000000000000 C
0.4444444444444466 0.8888888888888866 0.0000000000000000 C
0.7777777777777801 0.2222222222222200 0.0000000000000000 C
0.7777777777777799 0.5555555555555532 0.0000000000000000 C
0.7777777777777799 0.8888888888888866 0.0000000000000000 C
0.2222222222222200 0.1111111111111133 0.0000000000000000 C
0.2222222222222200 0.4444444444444466 0.0000000000000000 C
0.2222222222222199 0.7777777777777799 0.0000000000000000 C
0.5555555555555532 0.1111111111111133 0.0000000000000000 C
0.5555555555555534 0.4444444444444466 0.0000000000000000 C
0.5555555555555534 0.7777777777777799 0.0000000000000000 C
0.8888888888888866 0.1111111111111133 0.0000000000000000 C
0.8888888888888866 0.4444444444444466 0.0000000000000000 C
0.8888888888888866 0.7777777777777799 0.0000000000000000 C
Click to see supercell_4x4x1/POSCAR!
C32
1.0
9.8028876279774408 0.0000000000000000 0.0000000000000000
-4.9014438139887204 8.4895497162726397 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.0000000000000000
C
32
direct
0.0833333333333350 0.1666666666666650 0.0000000000000000 C
0.0833333333333350 0.4166666666666650 0.0000000000000000 C
0.0833333333333350 0.6666666666666650 0.0000000000000000 C
0.0833333333333349 0.9166666666666649 0.0000000000000000 C
0.3333333333333350 0.1666666666666650 0.0000000000000000 C
0.3333333333333350 0.4166666666666650 0.0000000000000000 C
0.3333333333333350 0.6666666666666650 0.0000000000000000 C
0.3333333333333349 0.9166666666666649 0.0000000000000000 C
0.5833333333333350 0.1666666666666650 0.0000000000000000 C
0.5833333333333350 0.4166666666666650 0.0000000000000000 C
0.5833333333333350 0.6666666666666650 0.0000000000000000 C
0.5833333333333349 0.9166666666666649 0.0000000000000000 C
0.8333333333333350 0.1666666666666650 0.0000000000000000 C
0.8333333333333350 0.4166666666666650 0.0000000000000000 C
0.8333333333333350 0.6666666666666650 0.0000000000000000 C
0.8333333333333349 0.9166666666666649 0.0000000000000000 C
0.1666666666666650 0.0833333333333350 0.0000000000000000 C
0.1666666666666650 0.3333333333333350 0.0000000000000000 C
0.1666666666666650 0.5833333333333349 0.0000000000000000 C
0.1666666666666649 0.8333333333333348 0.0000000000000000 C
0.4166666666666650 0.0833333333333350 0.0000000000000000 C
0.4166666666666650 0.3333333333333350 0.0000000000000000 C
0.4166666666666650 0.5833333333333349 0.0000000000000000 C
0.4166666666666649 0.8333333333333348 0.0000000000000000 C
0.6666666666666650 0.0833333333333350 0.0000000000000000 C
0.6666666666666650 0.3333333333333350 0.0000000000000000 C
0.6666666666666650 0.5833333333333349 0.0000000000000000 C
0.6666666666666649 0.8333333333333348 0.0000000000000000 C
0.9166666666666650 0.0833333333333350 0.0000000000000000 C
0.9166666666666650 0.3333333333333350 0.0000000000000000 C
0.9166666666666650 0.5833333333333349 0.0000000000000000 C
0.9166666666666649 0.8333333333333348 0.0000000000000000 C
Click to see supercell_5x5x1/POSCAR!
C50
1.0
12.2536095349718011 0.0000000000000000 0.0000000000000000
-6.1268047674859005 10.6119371453408000 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.0000000000000000
C
50
direct
0.0666666666666680 0.1333333333333320 0.0000000000000000 C
0.0666666666666680 0.3333333333333320 0.0000000000000000 C
0.0666666666666680 0.5333333333333320 0.0000000000000000 C
0.0666666666666680 0.7333333333333321 0.0000000000000000 C
0.0666666666666680 0.9333333333333320 0.0000000000000000 C
0.2666666666666680 0.1333333333333320 0.0000000000000000 C
0.2666666666666680 0.3333333333333320 0.0000000000000000 C
0.2666666666666680 0.5333333333333320 0.0000000000000000 C
0.2666666666666680 0.7333333333333321 0.0000000000000000 C
0.2666666666666679 0.9333333333333320 0.0000000000000000 C
0.4666666666666680 0.1333333333333320 0.0000000000000000 C
0.4666666666666680 0.3333333333333320 0.0000000000000000 C
0.4666666666666680 0.5333333333333320 0.0000000000000000 C
0.4666666666666680 0.7333333333333321 0.0000000000000000 C
0.4666666666666680 0.9333333333333320 0.0000000000000000 C
0.6666666666666681 0.1333333333333320 0.0000000000000000 C
0.6666666666666681 0.3333333333333320 0.0000000000000000 C
0.6666666666666681 0.5333333333333320 0.0000000000000000 C
0.6666666666666681 0.7333333333333321 0.0000000000000000 C
0.6666666666666681 0.9333333333333320 0.0000000000000000 C
0.8666666666666680 0.1333333333333320 0.0000000000000000 C
0.8666666666666680 0.3333333333333320 0.0000000000000000 C
0.8666666666666680 0.5333333333333320 0.0000000000000000 C
0.8666666666666680 0.7333333333333321 0.0000000000000000 C
0.8666666666666680 0.9333333333333320 0.0000000000000000 C
0.1333333333333320 0.0666666666666680 0.0000000000000000 C
0.1333333333333320 0.2666666666666680 0.0000000000000000 C
0.1333333333333320 0.4666666666666680 0.0000000000000000 C
0.1333333333333320 0.6666666666666681 0.0000000000000000 C
0.1333333333333319 0.8666666666666679 0.0000000000000000 C
0.3333333333333320 0.0666666666666680 0.0000000000000000 C
0.3333333333333319 0.2666666666666680 0.0000000000000000 C
0.3333333333333319 0.4666666666666680 0.0000000000000000 C
0.3333333333333320 0.6666666666666681 0.0000000000000000 C
0.3333333333333319 0.8666666666666679 0.0000000000000000 C
0.5333333333333320 0.0666666666666680 0.0000000000000000 C
0.5333333333333320 0.2666666666666680 0.0000000000000000 C
0.5333333333333320 0.4666666666666680 0.0000000000000000 C
0.5333333333333320 0.6666666666666681 0.0000000000000000 C
0.5333333333333319 0.8666666666666679 0.0000000000000000 C
0.7333333333333321 0.0666666666666680 0.0000000000000000 C
0.7333333333333321 0.2666666666666680 0.0000000000000000 C
0.7333333333333321 0.4666666666666680 0.0000000000000000 C
0.7333333333333321 0.6666666666666681 0.0000000000000000 C
0.7333333333333321 0.8666666666666679 0.0000000000000000 C
0.9333333333333319 0.0666666666666680 0.0000000000000000 C
0.9333333333333319 0.2666666666666680 0.0000000000000000 C
0.9333333333333319 0.4666666666666680 0.0000000000000000 C
0.9333333333333319 0.6666666666666681 0.0000000000000000 C
0.9333333333333319 0.8666666666666679 0.0000000000000000 C
Click to see supercell_6x6x1/POSCAR!
C72
1.0
14.7043314419661613 0.0000000000000000 0.0000000000000000
-7.3521657209830806 12.7343245744089586 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.0000000000000000
C
72
direct
0.0555555555555567 0.1111111111111100 0.0000000000000000 C
0.0555555555555566 0.2777777777777766 0.0000000000000000 C
0.0555555555555566 0.4444444444444433 0.0000000000000000 C
0.0555555555555567 0.6111111111111101 0.0000000000000000 C
0.0555555555555566 0.7777777777777765 0.0000000000000000 C
0.0555555555555567 0.9444444444444433 0.0000000000000000 C
0.2222222222222233 0.1111111111111100 0.0000000000000000 C
0.2222222222222233 0.2777777777777766 0.0000000000000000 C
0.2222222222222233 0.4444444444444433 0.0000000000000000 C
0.2222222222222233 0.6111111111111101 0.0000000000000000 C
0.2222222222222232 0.7777777777777765 0.0000000000000000 C
0.2222222222222233 0.9444444444444433 0.0000000000000000 C
0.3888888888888901 0.1111111111111100 0.0000000000000000 C
0.3888888888888899 0.2777777777777766 0.0000000000000000 C
0.3888888888888899 0.4444444444444433 0.0000000000000000 C
0.3888888888888900 0.6111111111111101 0.0000000000000000 C
0.3888888888888899 0.7777777777777765 0.0000000000000000 C
0.3888888888888900 0.9444444444444433 0.0000000000000000 C
0.5555555555555567 0.1111111111111100 0.0000000000000000 C
0.5555555555555567 0.2777777777777766 0.0000000000000000 C
0.5555555555555567 0.4444444444444433 0.0000000000000000 C
0.5555555555555567 0.6111111111111101 0.0000000000000000 C
0.5555555555555566 0.7777777777777765 0.0000000000000000 C
0.5555555555555567 0.9444444444444433 0.0000000000000000 C
0.7222222222222234 0.1111111111111100 0.0000000000000000 C
0.7222222222222233 0.2777777777777766 0.0000000000000000 C
0.7222222222222233 0.4444444444444433 0.0000000000000000 C
0.7222222222222234 0.6111111111111101 0.0000000000000000 C
0.7222222222222232 0.7777777777777765 0.0000000000000000 C
0.7222222222222234 0.9444444444444433 0.0000000000000000 C
0.8888888888888899 0.1111111111111100 0.0000000000000000 C
0.8888888888888898 0.2777777777777766 0.0000000000000000 C
0.8888888888888899 0.4444444444444433 0.0000000000000000 C
0.8888888888888899 0.6111111111111101 0.0000000000000000 C
0.8888888888888897 0.7777777777777765 0.0000000000000000 C
0.8888888888888899 0.9444444444444433 0.0000000000000000 C
0.1111111111111100 0.0555555555555567 0.0000000000000000 C
0.1111111111111100 0.2222222222222233 0.0000000000000000 C
0.1111111111111100 0.3888888888888899 0.0000000000000000 C
0.1111111111111100 0.5555555555555567 0.0000000000000000 C
0.1111111111111099 0.7222222222222232 0.0000000000000000 C
0.1111111111111100 0.8888888888888899 0.0000000000000000 C
0.2777777777777766 0.0555555555555567 0.0000000000000000 C
0.2777777777777767 0.2222222222222233 0.0000000000000000 C
0.2777777777777767 0.3888888888888899 0.0000000000000000 C
0.2777777777777767 0.5555555555555567 0.0000000000000000 C
0.2777777777777766 0.7222222222222232 0.0000000000000000 C
0.2777777777777767 0.8888888888888899 0.0000000000000000 C
0.4444444444444433 0.0555555555555567 0.0000000000000000 C
0.4444444444444433 0.2222222222222233 0.0000000000000000 C
0.4444444444444433 0.3888888888888899 0.0000000000000000 C
0.4444444444444433 0.5555555555555567 0.0000000000000000 C
0.4444444444444433 0.7222222222222232 0.0000000000000000 C
0.4444444444444433 0.8888888888888899 0.0000000000000000 C
0.6111111111111101 0.0555555555555567 0.0000000000000000 C
0.6111111111111101 0.2222222222222233 0.0000000000000000 C
0.6111111111111099 0.3888888888888899 0.0000000000000000 C
0.6111111111111101 0.5555555555555567 0.0000000000000000 C
0.6111111111111099 0.7222222222222232 0.0000000000000000 C
0.6111111111111101 0.8888888888888899 0.0000000000000000 C
0.7777777777777767 0.0555555555555567 0.0000000000000000 C
0.7777777777777766 0.2222222222222233 0.0000000000000000 C
0.7777777777777767 0.3888888888888899 0.0000000000000000 C
0.7777777777777767 0.5555555555555567 0.0000000000000000 C
0.7777777777777766 0.7222222222222232 0.0000000000000000 C
0.7777777777777767 0.8888888888888899 0.0000000000000000 C
0.9444444444444432 0.0555555555555567 0.0000000000000000 C
0.9444444444444433 0.2222222222222233 0.0000000000000000 C
0.9444444444444432 0.3888888888888899 0.0000000000000000 C
0.9444444444444432 0.5555555555555567 0.0000000000000000 C
0.9444444444444431 0.7222222222222232 0.0000000000000000 C
0.9444444444444433 0.8888888888888899 0.0000000000000000 C```
</details>
#computation of the phonon dispersion
PHON_NWRITE = -2
LPHON_DISPERSION = .TRUE.
LPHON_POLAR=.FALSE.
LPHON_READ_FORCE_CONSTANTS=.TRUE.
Contains the force constants in a format readable by VASP
wave vector q
50
line
reciprocal
0.00000000 0.00000000 0.00000000 GAMMA
0.50000000 0.00000000 0.00000000 M
0.50000000 0.00000000 0.00000000 M
0.33333333 0.33333333 0.00000000 K
0.33333333 0.33333333 0.00000000 K
0.00000000 0.00000000 0.00000000 GAMMA
Click to see supercell_3x3x1/KPOINTS!
K points
0
Gamma
4 4 1
0 0 0
Click to see supercell_4x4x1/KPOINTS!
K points
0
Gamma
3 3 1
0 0 0
Click to see supercell_5x5x1/KPOINTS!
K points
0
Gamma
2 2 1
0 0 0
Click to see supercell_6x6x1/KPOINTS!
K points
0
Gamma
2 2 1
0 0 0
PAW C_s 04May1998
4.00000000000000
Check the meaning of all tags used in the INCAR file and make comments in the file!
3.3 Calculation¶
Open a terminal, navigate to this example's directory, copy the vaspout.h5 and POSCAR file of the previous example to vaspin.h5 and run VASP by entering the following:
cd $TUTORIALS/phonon/e03_*
for d in 3x3x1 4x4x1 5x5x1 6x6x1
do
cd supercell_$d
cp ../../e02_*/supercell_$d/POSCAR POSCAR
cp ../../e02_*/supercell_$d/vaspout.h5 vaspin.h5
mpirun -np 1 vasp_std
cd ..
done
Plot the phonon dispersion using py4vasp!
import py4vasp
supercell3 = py4vasp.Calculation.from_path( "./e03_dispersion/supercell_3x3x1")
supercell4 = py4vasp.Calculation.from_path( "./e03_dispersion/supercell_4x4x1")
supercell5 = py4vasp.Calculation.from_path( "./e03_dispersion/supercell_5x5x1")
supercell6 = py4vasp.Calculation.from_path( "./e03_dispersion/supercell_6x6x1")
supercell3.phonon_band.plot()
supercell4.phonon_band.plot()
supercell5.phonon_band.plot()
supercell6.phonon_band.plot()
With the following code you can get the phonon dispersions for all the different supercell sizes on the same plot
import py4vasp
from plotly.subplots import make_subplots
calculations = ['3x3','4x4','5x5','6x6']
#create plotly figure
final_fig = make_subplots()
for i,label in enumerate(calculations):
fig = py4vasp.Calculation.from_path(f"./e03_dispersion/supercell_{label}x1").phonon_band.to_plotly()
for d in fig.data:
d.name = label
d['line'].color=None # reset the color scheme
final_fig.add_traces(fig.data)
final_fig.update_layout(fig.layout)
final_fig.update_layout(autosize=False,width=800,showlegend=True,margin={"l":0,"r":0,"t":0,"b":0})
What are the differences between the different phonon dispersions?
Near the Γ point, a saddle point appears when the supercell size is too small. These slightly imaginary modes might appear again for larger supercell sizes but tend to disappear when the supercell size is large enough. The convergence of these quadratic 'soft' phonon modes in 2D materials is particularly difficult and a lot of care should be taken when performing these calculations.
Negative phonon frequencies indicate of instability in the structure. While the forces were zero in the equilibrium structure, the distortion according to a certain phonon mode lowers the total energy. This implies that the structure is metastable (given the approximations and limitations of the current approach: 0 K and harmonic approximation). Sometimes these instabilities are signatures of an inaccurate calculation or an artifact of the interpolation of the phonon dispersion using a truncated set of real space force constants (finite supercell size). Sometimes they have physical meaning and distorting according to that phonon mode and performing ionic relaxation leads to a structure with different symmetry and lower energy.
In the present calculation (and in most cases) these negative modes near the Γ point are indicative of a problem with the calculation. This can be insufficient supercell size or too low ENCUT.
Look at the phonon frequencies at the Γ point for the different supercell sizes.
Why does the 5x5 supercell yield optical phonon frequencies different from all the other calculations?
While for the 3x3, 4x4, and 6x6 supercells we kept the k point sampling consistent (i.e. decreasing with increasing supercell size) for the 5x5 supercell this is not possible so we used a coarser k point mesh. The difference between the frequencies at Γ for 5x5 and the other supercells will be small when the k-point sampling is dense enough.
As we have mentioned before, at the q points commensurate with the supercell the phonon frequencies are computed exactly. All the other points are obtained through interpolation. To illustrate this better, we show how increasing the size of the supercell corresponds to sampling the Brillouin zone with a Γ centered k point mesh. Shown below are the plots of the reciprocal Brillouin zone with the k points generated from a regular sampling of the Brillouin zone corresponding to a certain supercell size.
For a given supercell if we interpolate the dynamical matrix along a path in reciprocal space that passes on the points in blue and diagonalize it, we will obtain the exact phonon frequencies (given the current approximations: 0 K, harmonic approximation, DFT, etc..) only on those blue points. Everywhere else the results are obtained through interpolation.
What is the minimum supercell size we need to compute the phonon frequencies at the M point?
2x2
What is the minimum supercell size we need to compute the phonon frequencies at the K point?
3x3
Now you can interactively visualize the phonon modes. This is done using the phononwebsite.
Executing the following code will use the phononwebsite to visualize the phonon modes. First a .json file will be created from the vasput.h5 file.
#create json file from hdf5
from phononweb.vaspphonon import VaspPhonon
vp = VaspPhonon("graphene",filename="e03_dispersion/supercell_4x4x1/vaspout.h5")
vp.write_json()
This will create a graphene.json file in the current working directory. Download this file and visualize it using the Choose-file option on https://henriquemiranda.github.io/phononwebsite/phonon.html.
On the left panel, you can see the structural information, choose the number of repetitions of the unit cell to be displayed, change the camera, change the amplitude of the oscillations or display vectors in the directions where the atoms are moving.
The middle panel shows the structure oscillating according to some phonon modes. You can move it around by dragging it with the mouse or zoom in/out by scrolling.
The right panel shows the phonon dispersion. By clicking on any point in it you can select which phonon mode should be animated in the central panel.
3.4 Questions¶
- How far does the wavelength of $\mathbf q=0$ phonons extend in the crystal?
- What is phonon dispersion?
- How are the force constants and the phonon frequencies connected?