Phonons with Density-Functional Perturbation Theory
module load qe/5.0
tar -zxvf examples_phon.tgz cd examples_phon
Phonons at Gamma point (q=0)All phonon calculations start from the self-consistent charge density and Kohn-Sham orbitals, calculated for equilibrium atomic positions. For q=0 phonons in Si, you just need to run the phonon code afterwards.
pw.x < si.scf.in > si.scf.outYou will get the usual self-consistent results for Si, and the needed data written to directory "outdir"/"prefix".save
ph.x < si.phG.in > si.phG.outThe 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-pointIn order to perform a phonon calculation at a generic (nonzero) q-point, an intermediate non-scf calculation must be run to (re)construct electronic states at k and k+q. This could be performed as a separate step, but it can more simply be done by the phonon code. In the following, the X point, q = (1,0,0) 2pi/a, is used. The steps are
0.0 0.0 0.0 ===> 1.0 0.0 0.0Sample file: si.phX.in.
ph.x < si.phX.in > si.phX.outIn 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.
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 RulesLet 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 if q=(0.0 0.0 0.0), otherwise the code complains. Run the phonon calculation at Gamma again:
ph.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:
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
Phonons in Polar MaterialsIn 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 direction q. We consider AlAs, a polar semiconductor with the zincblende structure (i.e. like Si but one different atom per sublattices; Si is nonpolar so there is no LO-TO splitting). Create an input for scf calculation in AlAs using the following data:
ATOMIC_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 q=>0. The 3-fold degenerate optical modes are now split by the TO-LO splitting. You should get two degenerate TO frequencies around 10.71 THz, one LO frequency at 11.79 THz. Note that asr='simple' (or 'three-dim', a more sophisticated version) makes the frequency of the three lowest acoustic modes at q=0 to vanish, but has no effect on the other modes.
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).