Phonons with Density-Functional Perturbation Theory

Getting ready

- A short reminder of the relevant theory can be found here.
- The following command must have been executed:
module load qe/5.0

- You may download the entire exercise file
examples_phon.tgz
in a directory of your choice. Uncompress and unpack
the file and enter in the resulting directory:
tar -zxvf examples_phon.tgz cd examples_phon

Phonons at Gamma point (**q**=0)

- Calculate the self-consistent charge density and Kohn-Sham
orbitals for Si. Sample input file: si.scf.in.
Edit "pseudo_dir" and "outdir" so that they point respectively to the
directory containing the pseudopotential (you will need
Si.pz-vbc.UPF) and to the scratch directory.
pw.x < si.scf.in > si.scf.out

You will get the usual self-consistent results for Si, and the needed data written to directory "outdir"/"prefix".save - Calculate phonons at wavevector
**q**=(0,0,0). Sample input file: si.phG.in. The values of "outdir" and "prefix" must be the same as in the preceding scf step. The value of the mass will overwrite the value given in the scf calculation. The line after the namelist is the q-point to be computed (gamma).ph.x < si.phG.in > si.phG.out

The code should perform 3*Nat, Nat=Number of atoms, linear response calculations, one per atomic displacement. For efficiency reason, the 3*Nat atomic displacements are grouped into "irreps" (irreproducible representations), based on the symmetry of the system. At the end the dynamical matrix is computed, written to file si.dynG [i.e. the variable fildyn in "si.phG.in"], diagonalized. Note the 3 degenerate non zero frequencies and the 3 frequencies close to zero. The latter should be exactly zero due to Acoustic Sum Rule (ASR)), i.e. traslational symmetry. Notice that the degeneracy of the modes [ 3 + 3 ] corresponds to the dimensionality of the irreps: it is a property related to symmetry.

dynmat.x &input fildyn='si.dynG' /You may also impose the Acoustic Sum Rule, as in e.g.:

dynmat.x &input fildyn='si.dynG', asr='simple' /Notice the change. File "dynmat.axsf" contains normal modes in a format that can be visualized by XCrySDen: they appear as forces on atoms. Use

xcrysden --axsf dynmat.axsfto visualize it, show forces (f) and change force options (shift-f) for a better visualization. Do they look reasonable? what did you expect?

Phonons at a generic **q**-point

- Run the calculation of the self-consistent charge density and Kohn-Sham orbitals for Si,if you haven't already, or if the data files have been lost for some reason (same as above)
- Calculate phonons at wavevector
**q**=(1,0,0) 2pi/a. Copy "si.phG.in" into "si.phX.in" and edit it: modify the title if you wish, the "fildyn" variable ("si.dynG" -> "si.dynX") and the last line, defining the q-point:0.0 0.0 0.0 ===> 1.0 0.0 0.0

Sample file: si.phX.in.ph.x < si.phX.in > si.phX.out

In the first part of the output, a non-scf calculation is performed. The symmetry of the small group of**q**is assumed instead of the crystal symmetry, so the number of**k**-points increases. Note the presence of**k**+**q**vectors intercalated between the**k**vectors.

The linear-response calculation is performed in the second part of the output. Notice that the number and degeneracies of irreps are different from those at**q**=0: they depend on symmetry, i.e. on the small group of**q**. The file si.dynX contains now three matrices, i.e. the three equivalent X points.

grep omega si.dyn* si.dynG: omega( 1) = 0.063598 [THz] = 2.121419 [cm-1] si.dynG: omega( 2) = 0.063598 [THz] = 2.121419 [cm-1] si.dynG: omega( 3) = 0.063598 [THz] = 2.121419 [cm-1] si.dynG: omega( 4) = 15.475124 [THz] = 516.197992 [cm-1] si.dynG: omega( 5) = 15.475124 [THz] = 516.197992 [cm-1] si.dynG: omega( 6) = 15.475124 [THz] = 516.197992 [cm-1] si.dynL: omega( 1) = 3.281779 [THz] = 109.469103 [cm-1] si.dynL: omega( 2) = 3.281779 [THz] = 109.469103 [cm-1] si.dynL: omega( 3) = 11.303519 [THz] = 377.047327 [cm-1] si.dynL: omega( 4) = 12.548573 [THz] = 418.578134 [cm-1] si.dynL: omega( 5) = 14.786637 [THz] = 493.232384 [cm-1] si.dynL: omega( 6) = 14.786637 [THz] = 493.232384 [cm-1] si.dynX: omega( 1) = 4.317843 [THz] = 144.028682 [cm-1] si.dynX: omega( 2) = 4.317843 [THz] = 144.028682 [cm-1] si.dynX: omega( 3) = 12.400213 [THz] = 413.629325 [cm-1] si.dynX: omega( 4) = 12.400213 [THz] = 413.629325 [cm-1] si.dynX: omega( 5) = 13.962671 [THz] = 465.747664 [cm-1] si.dynX: omega( 6) = 13.962671 [THz] = 465.747664 [cm-1]

Electric Fields and Acoustic Sum Rules

Let us consider a little bit more in detail the Gamma case. In insulating materials the system can sustain macroscopic electric fields and these can be coupled with vibrations. The physical properties that describe these properties are the dielectric tensor ("epsilon") and the effective charges ("Z*"). Note that in metals there are no macroscopic electric fields: the dielectric constant is infinite. Edit "si.phG.in" and add variable epsil=.true. to the namelist. This can be done ONLY ifph.x < si.phG.in > si.phG.outNote that an additional linear-response calculation for electric fields in 3 independent directions is performed, dielectric constant and effective charges are calculated and stored into the "si.dynG" file. Note that effective charges are tensors: they are neither diagonal nor even symmetric in general! In high-symmetry systems like this one they reduce to a scalar, though. Due to translational symmetry, ASR requires that acoustic modes have zero frequencies AND the sum of effective charges over all atoms is zero: \sum_i Z*_i = 0. In Si, where the two atoms in the cell are equivalent, this means Z*=0. This is never exactly verified in realistic calculations because:

- insufficiently accurate scf threshold (in pw.x and/or ph.x part).
Usually pw.x takes very little time w.r.t. ph.x, so
**do not**use poor (large) thresholds in pw.x in order to spare some time: you gain very little and you may get lousy results. - XC energy is computed on a real-space grid that assumes an explicit origin. More problematic for GGA than for LDA.
**k**-point sampling is not accurate enough (for Z* and epsilon): macroscopic electric fields arise from the limit**q**->0, which is known to require a dense grid of k-points.

K_POINTS automatic 6 6 6 1 1 1Run again pw.x+ph.x, save the output for later comparison:

pw.x < si.scf.in-28k > si.scf.out-28k ph.x < si.phG.in > si.phG.out-28kDo the same for the 2 chadi-cohen / monkhorst-pack points :

K_POINTS automatic 2 2 2 1 1 1for the 6 monkhorst-pack points :

K_POINTS automatic 3 3 3 1 1 1and for the 10 chadi-cohen / monkhorst-pack points (what we have used up to now)

K_POINTS automatic 4 4 4 1 1 1Notice that

- ASR for effective charges improves by improving BZ sampling.
- ASR for acoustic mode improves but does not go to zero => k-point sampling is not an issue for ASR for acoustic modes.
- Optical modes are practically unaffected.

Phonons in Polar Materials

In polar materials, macroscopic electric fields are present in the long-wavelength limit and induce the so-called LO-TO splitting. This can be computed from the knowledge of epsilon, Z* and the propagation directionATOMIC_SPECIES Al 26.98 Al.pz-vbc.UPF As 74.92 As.pz-bhs.UPFand alat = 10.6 (theoretical lattice parameter for LDA and those pseudopotentials). Sample file: alas.scf.in. You will need the pseudopotentials for Al, Al.pz-vbc.UPF, and for As, As.pz-bhs.UPF. Run the scf calculation as usual:

pw.x < alas.scf.in > alas.scf.outCreate an input for phonon at Gamma for AlAs (including electric fields: "epsil=.true.") Sample file: alas.phG.in.

ph.x < alas.phG.in > alas.phG.outNotice that the optical phonon frequencies are 3-fold degenerate and there is no TO-LO splitting, because the non-analytic term that produces the TO-LO splitting has to be added to the dynamical matrix. This is done by auxiliary program "dynmat.x". Create an input (alas.dynmat.in for instance) for dynmat.x, specifying the additional variables

&input fildyn='alas.dynG', asr='simple' q(1)=1.0, q(2)=0.0, q(3)=0.0 /Run the dynmat.x code:

dynmat.x < alas.dynmat.inThe code will take note of the presence of Z* and epsilon in the dynamical matrix file and will add to the dynamical matrix the nonanalytic part corresponding to the given

In cubic system, such as AlAs, the frequencies (NOT the eigenvectors)
are independent on the **q** direction. In lower symmetries both
eigenvalues and eigenvectors depend on the phonon propagation direction,
and modes cannot be distinguished between purely transverse (i.e.
**q** perpendicular to the displacement pattern) and purely
longitudinal (i.e. **q** parallel to the displacement pattern).