# Examples of calculations with ZG.x¶

Author: Marios Zacharias

Note

It is perfectly reasonable to obtain ZG displacements that differ from those reported in this tutorial. This difference arises from the gauge freedom of the vibrational modes: modes obtained by diagonalizing the dynamical matrix can differ by a phase factor (or a unitary matrix in case of degeneracy) if calculated on different systems. The final results are not affected by these differences.

## Running ZG.x for evaluating temperature-dependent properties with Quantum Espresso¶

The tutorial is based on the **Special Displacement Method** described in the papers Phys. Rev. Research 2, 013357 (2020) and Phys. Rev. B 94, 075125 (2016). Tutorials on anharmonic lattice dynamics are based on the paper Phys. Rev. B 108, 035155 (2023). It is advisable to check them for the interpretation of the results and description of the main equations related to this tutorial.

The results in this example are obtained with Quantum_ESPRESSO v7.2 and the ZG.x code in EPW.

For this tutorial purposes set the following environment variable in your shell and scripts:

```
QE=/path_to_q-e/ # this is the path to q-e
```

Executables in `$QE/EPW/ZG/`

directory are compiled together with `make epw`

. To use local executables go to `$QE/EPW/ZG/src/local`

and type `./compile_ifort.sh`

.

In this tutorial we will show how to obtain all input and output files. We also provide inputs as a reference or to speed up the process. Most of the temperature-dependent calculations will be performed for *T = 0* K; one can repeat the steps using a different temperature.

Exercise 1: Generation of the ZG configuration and DFT-ZG calculation.

Exercise 2: Temperature-dependent band structures.

Exercise 3: Anharmonic lattice dynamics.

Exercise 4: Temperature-dependent optical spectra.

Exercise 5: Multiphonon diffuse scattering intensity.

Exercise 6: Phonon unfolding.

A complete list of tutorials and files are available in this link.

## Exercise 1¶

The tutorial tarball is located inside the `$QE/EPW/ZG/`

directory. Untar the tutorial tarball, go to `exercise1`

, and create directory `workdir`

:

```
tar -xvf tutorial.tar.gz; cd tutorial/exercise1; mkdir workdir
```

In this exercise we will generate the ZG configuration of silicon in a \(3\times3\times3\) supercell for the temperature *T = 0* K and run your first DFT-ZG calculation. In the following, all steps for obtaining the phonons of silicon are provided, but to speed up the process you can skip the first four steps. These steps are standard Quantum Espresso runs for obtaining the interatomic force constants. We note that for generating successfully a ZG configuration you need to make sure that the phonon dispersion is as you would expect (compare with literature). If there exist negative frequencies (soft modes), the code excludes them. If the system is anharmonic one should consider to evaluate anharmonic phonons using A-SDM as in Exercise 3.

Run a self-consistent calculation for silicon in the `workdir`

. For converged results, the energy cutoff **ecutwfc** should be 30 Ry:

```
cd workdir; cp ../inputs/si.scf.in .; cp ../inputs/Si.pz-vbc.UPF .
mpirun -np 4 $QE/bin/pw.x -nk 4 < si.scf.in > si.scf.out
```

Run a `ph.x`

calculation on a homogeneous \(4\times4\times4\) **q**-point grid using the input:

```
-- si.ph.in
Phonons of Silicon
&inputph
amass(1) = 28.0855,
prefix = 'si'
outdir = './'
ldisp = .true.
fildyn = 'si.dyn'
tr2_ph = 1.0d-12
nq1 = 4, nq2 = 4, nq3 = 4
/
```

which can be found from the `inputs`

directory:

```
cp ../inputs/si.ph.in .
mpirun -np 4 $QE/bin/ph.x -nk 4 < si.ph.in > si.ph.out
```

This will generate 8 **si.dynX** output files containing the dynamical matrix calculated for each irreducible **q**-point. The list of irreducible **q**-points is written in the **si.dyn0** file.

Run a `q2r.x`

calculation to obtain the interatomic force constants (IFC) file using the input:

```
-- q2r.in
&input
fildyn='si.dyn', flfrc = 'si.444.fc'
/
```

i.e.:

```
cp ../inputs/q2r.in .
mpirun -np 1 $QE/bin/q2r.x < q2r.in > q2r.out
```

This will generate **si.444.fc** which contains the IFCs. As we will see below, this is the only external input file necessary for running `ZG.x`

.

Run a `matdyn.x`

calculation to check the phonon dispersion using the input file:

```
-- matdyn.in
&input
asr='simple', amass(1)=28.0855,
flfrc='si.444.fc', fldyn='si.dyn.mat', flfrq='si.freq', fleig='si.dyn.eig',
q_in_cryst_coord = .false., q_in_band_form = .true.
/
9
0.00 0.00 0.00 100
0.75 0.75 0.00 1
0.25 1.00 0.25 100
0.00 1.00 0.00 100
0.00 0.00 0.00 100
0.50 0.50 0.50 100
0.00 1.00 0.00 100
0.50 1.00 0.00 100
0.50 0.50 0.50 100
```

and

```
cp ../inputs/matdyn.in .
mpirun -np 1 $QE/bin/matdyn.x < matdyn.in > matdyn.out
$QE/bin/plotband.x
Input file > si.freq
Reading 6 bands at 702 k-points
Range: 0.0000 509.7631eV Emin, Emax, [firstk, lastk] > 0 600
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 0.0000
high-symmetry point: 0.7500 0.7500 0.0000 x coordinate 1.0607
...
high-symmetry point: 0.5000 0.5000 0.5000 x coordinate 5.3534
output file (gnuplot/xmgr) > si_ph_dispersion.xmgr
bands in gnuplot/xmgr format written to file si_ph_dispersion.xmgr
output file (ps) >
stopping ...
```

`matdyn.x`

will generate **si.freq** which contains the phonon frequencies along the specified path in the Brillouin zone. Then we combine **si.freq** with `plotband.x`

to obtain **si_ph_dispersion.xmgr** in gnuplot/xmgr format for the phonon dispersion. To plot the phonon dispersion, as shown in Figure 1, type:

```
cp si_ph_dispersion.xmgr ../gnuplot/.
cd ../gnuplot/; gnuplot gp_coms.p; evince ph_dispersion_444.eps
```

Is the phonon dispersion the one you would expect? One can verify the results with literature data. Once you have obtained the IFC file (e.g. here **si.444.fc**) and verified the correctness of your phonons, you are ready to perform a `ZG.x`

calculation. If you have skipped the four steps above, you can copy the IFC file from `inputs`

into `workdir`

:

```
cd ../workdir; cp ../inputs/si.444.fc .
```

Run a ZG calculation using the input:

```
-- ZG_333.in
&input
flfrc='si.444.fc',
asr='simple', amass(1)=28.0855, atm_zg(1) = 'Si',
T = 0.00,
dim1 = 3, dim2 = 3, dim3 = 3
compute_error = .true., synch = .true., error_thresh = 0.035, niters = 15000
incl_qA = .false.
/
```

Details about all input flags are provided in https://docs.epw-code.org/doc/InputsZG.html. The **q**-point grid (**nq1** * **nq2** * **nq3**) used to obtain the dynamical matrices should not be necessarily the same with the supercell size (**dim1** * **dim2** * **dim3**). In particular, **dimX** is independent of **nqX**, since `ZG.x`

takes advantage of Fourier interpolation as implemented in `matdyn.x`

; thus any size of ZG configuration can be generated from **si.444.fc**.

```
cp ../inputs/ZG_333.in .
mpirun -np 4 $QE/bin/ZG.x -nk 4 < ZG_333.in > ZG_333.out
```

The outputs from this `ZG.x`

run are: (a) **ZG-configuration_0.00K.dat** and **equil_pos.dat** which contain the ZG (quantum nuclei at 0 K) and equilibrium (classical nuclei at 0 K) coordinates in angstroms, and (b) **ZG-scf_333_0.00K.in** and **equil-scf_333.in** which are the corresponding scf files for performing calculations in a \(3\times3\times3\) supercell.

Run a self-consistent calculation for silicon with the nuclei clamped at their equilibrium and ZG coordinates (for *T* = 0.00 K) in the supercell:

```
mpirun -np 54 $QE/pw.x -nk 3 < equil-scf_333.in > equil-scf_333.out
mpirun -np 56 $QE/pw.x -nk 8 < ZG-scf_333_0.00K.in > ZG-scf_333_0.00K.out
```

The input file `equil-scf_333.in`

is also provided here:

```
&control equil-scf_333.in
calculation = 'scf'
restart_mode = 'from_scratch',
prefix = 'equil-si',
pseudo_dir = './',
outdir='./'
/
&system
ibrav = 0
nat = 54
ntyp = 1
ecutwfc = 20.00
occupations = 'fixed'
smearing = 'gaussian'
degauss = 0.0D+00
/
&electrons
diagonalization = 'david'
mixing_mode= 'plain'
mixing_beta = 0.70
conv_thr = 0.1D-06
/
ATOMIC_SPECIES
Si 28.086 Si.pz-vbc.UPF
K_POINTS automatic
2 2 2 0 0 0
CELL_PARAMETERS {angstrom}
-8.09641133 0.00000000 8.09641133
0.00000000 8.09641133 8.09641133
-8.09641133 8.09641133 0.00000000
ATOMIC_POSITIONS (angstrom)
Si 0.00000000 0.00000000 0.00000000
Si -1.34940189 1.34940189 1.34940189
Si -2.69880378 0.00000000 2.69880378
...
```

Note

Cell parameters, number of atoms, **k**-grid, and atomic coordinates have been modified automatically based on the supercell dimensions. The code will always generate the lattice information in angstroms.

The file **ZG-scf_333_0.00K.in** is also provided here:

```
&control ZG-scf_333_0.00K.in
calculation = 'scf'
restart_mode = 'from_scratch'
prefix = 'ZG-0.00K-si'
pseudo_dir = './'
outdir = './'
/
&system
ibrav = 0
nat = 54
ntyp = 1
ecutwfc = 20.00
occupations = 'fixed'
smearing = 'gaussian'
degauss = 0.0D+00
/
&electrons
diagonalization = 'david'
mixing_mode= 'plain'
mixing_beta = 0.70
conv_thr = 0.1D-06
/
ATOMIC_SPECIES
Si 28.085 Si.pz-vbc.UPF
K_POINTS automatic
2 2 2 0 0 0
CELL_PARAMETERS (angstrom)
-8.09641133 0.00000000 8.09641133
0.00000000 8.09641133 8.09641133
-8.09641133 8.09641133 0.00000000
ATOMIC_POSITIONS (angstrom)
Si -0.02610870 -0.00524574 0.02749582
Si -1.36329713 1.31619849 1.42481245
Si -2.75218400 -0.00549830 2.71681282
...
```

Note

Same information with `equil-scf 333.in`

, except for the prefix and atomic positions. In general, you should add `nosym = .true.`

, since symmetries do not apply after ZG displacements. The information printed by the code for constructing these scf files is basic; so for more complex cases one can modify the files based on their calculation needs.

Let us now discuss some apsects about the generation of the **ZG-configuration** file. Type:

```
grep -3 "Optimum" ZG_333.out | tail -4
```

This prints:

```
Optimum configuration found !
Sum of diagonal terms per q-point: 0.419087
Error and niter index: 0.032411 7403
```

This information is printed by `ZG.x`

when the optimum configuration is found. Line 3 gives the sum of all diagonal terms per q-point, representing the denominator of Eq. (54) in Phys. Rev. Research 2, 013357, (2020). The first entry in line 4 represents the value of the minimization function E(S_{q,ν}) given by Eq. (54) in Phys. Rev. Research 2, 013357, (2020). The integer in line 4 represents the number of attempts required to achieve E(S_{q,ν}) smaller than the value specified for `error_thresh`

flag. If this number exceeds the integer specified by `niters`

flag, then `ZG.x`

will stop without printing the ZG-configuration file. A general rule is: the larger the supercell, the fewer attempts are required to achieve E(S_{q,ν}) < `error_thresh`

. The reason is: the order of choosing the unique set of signs, assigned to the **q**-points, is less important as we approach the thermodynamic limit.

Now type:

```
head -3 ZG-configuration_0.00K.dat
```

This prints:

```
Temperature is: 0.00 K
Atomic positions 54
Si -0.02610870 -0.00524574 0.02749582
```

Lines 1 and 2 give the temperature for which ZG displacements are generated and the total number of atomic coordinates, respectively. In line 3, we have the first ZG atomic coordinate. It is
perfectly reasonable to find different ZG coordinates, since the modes obtained by diagonalizing the dynamical matrix can differ by a phase factor (or a unitary matrix in case of degeneracy) if the processor, or compiler, or libraries have changed. The eigenvalues, of course, should remain the same in all cases. The flag `synch = .true.`

should apply a smooth Berry connection and align the sign of the modes with respect to a reference mode, but the sign of this reference depends on the machine. Degeneracy is not taken into account. The best way to check the validity of your configuration is by comparing the anisotropic mean-squared displacement tensor with the exact values, both obtained at the end of the **ZG_333.out**. To this aim type:

```
grep -12 "Anisotropic" ZG_333.out | tail -13
```

This prints:

```
Anisotropic mean-squared displacement tensor vs exact values (Ang^2):
Atom 1
Si 0.002243 0.002332 0.002322
Exact values
Si 0.002352 0.002352 0.002352
Atom 2
Si 0.002529 0.002407 0.002277
Exact values
Si 0.002352 0.002352 0.002352
off-diagonal terms
Si 0.000189 0.001294 0.000610
Si 0.000652 -0.000851 0.000027
```

Reducing `error_thresh`

brings the anisotropic displacement tensor closer to the exact values.

Now check whether the calculations have finished and type:

```
grep ! equil-scf_333.out
grep ! ZG-scf_333_0.00K.out
```

This prints:

```
equil-scf_333.out:! total energy = -427.83892824 Ry
ZG-scf_333_0.00K.out:! total energy = -427.72018072 Ry
```

These are the total energies (Kohn-Sham potential energy surfaces) of the equilibrium and ZG structures. What is the reason of the difference between these two values? The difference is ΔE = 0.1187 Ry and is due to the vibrational potential energy. To check this type:

```
grep -5 "Potent" ZG_333.out
```

This prints:

```
==============================================
Total vibrational energy: 0.24020789 Ry
Potential energy: 0.12010394 Ry
Kinetic energy: 0.12010394 Ry
Vibrational free energy: 0.24020789 Ry
```

Our computed *ΔE* is indeed almost equal to half the total vibrational energy, since a DFT calculation will not account for the vibrational kinetic energy contribution. The value 0.119 deviates from 0.12, since (i) we exclude from our calculations the **q**-points belonging in set *A* and (ii) due to the error coming by using small ZG configurations. As we use larger supercells this small deviation reduces to zero and *ΔE = Total vibrational energy / 2*. The code also prints the vibrational free energy evaluated using Eq. (A5) of Phys. Rev. B 108, 035155 (2023). The information for the vibrational free energy is important for checking the convergence of the trial free energy in A-SDM calculations (Exercise 3).

For the next exercises we will only need the charge-density files from **equil-si.save** and **ZG-0.00K-si.save**. Thus remove:

```
rm -r *wfc* si.save _ph0/ equil-si.save/*wfc* ZG-0.00K-si.save/*wfc*
```

## Exercise 2¶

In this exercise we will learn how to calculate the phonon-induced band gap renormalization and temperature-dependent band structures using the example of silicon and the math:3times3times3 ZG supercell. To obtain temperature-dependent band structures we will employ the band structure unfolding (BSU) technique with plane-waves as basis sets. The theory of BSU with plane-waves is described in Phys. Rev. B 85, 085201 (2012). To speed up the process one can skip the following two steps.

First go to the directory `exercise2`

and copy the following input files in your `workdir`

:

```
cd ../../exercise2/; mkdir workdir; cd workdir
cp ../inputs/Si.pz-vbc.UPF .; cp ../inputs/si.scf.in .
cp ../inputs/si.bands.in .; cp ../inputs/bands.in .
```

In file **si.bands.in** we set `nbnd = 5`

to include one empty band. We also sample the Γ-X path using 50 **k**-points in Cartesian coordinates (units of 2 pi/a).

Run a standard band structure calculation in the unit-cell of silicon along Γ-X using 50 **k**-points:

```
mpirun -np 4 $QE/bin/pw.x -nk 4 < si.scf.in > si.scf.out
mpirun -np 4 $QE/bin/pw.x -nk 4 < si.bands.in > si.bands.out
mpirun -np 4 $QE/bin/bands.x -nk 4 < bands.in > bands.out
```

Plot the band structure as shown in Figure 2 using the file **bands.dat.gnu**:

```
cp bands.dat.gnu ../gnuplot/; cd ../gnuplot/; gnuplot gp_coms.p
evince band_structure.eps
```

Can you find the k-coordinates of the valence band maximum (VBM) and conduction band minimum (CBM)? Inspecting file **bands.dat**, those are:

```
0.000000 0.000000 0.000000
-5.811 6.253 6.253 6.253 8.817
...
0.000000 0.840000 0.000000
-2.783 -0.278 3.440 3.440 6.743
```

We can see that the VBM lies at the Γ-point and is threefold degenerated with energy E_VBM = 6.253 eV. The CBM lies at ~0.84 Γ-X with energy E_CBM = 6.743 eV. Therefore, the DFT-LDA band gap with these settings is E_G = E_CBM - E_VBM = 0.49 eV.

Run a non self-consistent calculation for two **K**-points that define the VBM and CBM in the equilibrium and ZG supercell structures. The selection rule is to find the reciprocal lattice vector of the supercell **G** that maps **K** onto **k**, i.e. **K** = **k** - **G**. The trivial choice is **G** = **Γ** and thus set **K** = **k**.

To prepare the calculation proceed as follows:

```
cd ../workdir/
cp ../../exercise1/workdir/equil-scf_333.in equil-nscf_333.in
cp ../../exercise1/workdir/ZG-scf_333_0.00K.in ZG-nscf_333_0.00K.in
cp -r ../../exercise1/workdir/equil-si.save/ .
cp -r ../../exercise1/workdir/ZG-0.00K-si.save/ .
```

In **equil-nscf_333.in** and **ZG-nscf_333_0.00K.in** apply the following changes:

```
1. calculation = 'nscf'
2. nbnd = 135 # i.e. include one empty band / unit-cell
3. Replace the automatic K-grid:
K_POINTS automatic
2 2 2 0 0 0
with the two K-points of the VBM and CBM:
K_POINTS crystal
2
0.000000 0.000000 0.000000 1
0.000000 1.260000 1.260000 1
4. For ZG input file add nosym = .true. after nbnd = 135.
```

Since the coordinates are given in units of the reciprocal lattice parameters, we obtain the **K**-coordinates by multiplying the **k**-coordinates of VBM and CBM states with the dimensions of the supercell, i.e. if **k** *=* [ *x y z*] then **K** *=* [ \(m \times x n \times y p \times z\) ], where integers \(m, n, p\) define an \(m \times n \times p\) supercell.

If you did not complete Exercise 1 and the steps above, then copy the files from the `inputs`

directory:

```
cd ../workdir/; cp ../inputs/*nscf_333*in .; cp -r ../inputs/*-si.save/ .
```

For example the ZG input file should look like this:

```
&control ZG-nscf_333_0.00K.in
calculation = 'nscf'
restart_mode = 'from_scratch'
prefix = 'ZG-0.00K-si'
pseudo_dir = './'
outdir = './'
/
&system
ibrav = 0
nat = 54
ntyp = 1
ecutwfc = 20.00
nbnd = 135
nosym = .true.
occupations = 'fixed'
smearing = 'gaussian'
degauss = 0.0D+00
/
&electrons
diagonalization = 'david'
mixing_mode= 'plain'
mixing_beta = 0.70
conv_thr = 0.1D-06
/
ATOMIC_SPECIES
Si 28.085 Si.pz-vbc.UPF
K_POINTS crystal
2
0.000000 0.000000 0.000000 1
0.000000 1.260000 1.260000 1
CELL_PARAMETERS (angstrom)
-8.09641133 0.00000000 8.09641133
0.00000000 8.09641133 8.09641133
-8.09641133 8.09641133 0.00000000
ATOMIC_POSITIONS (angstrom)
...
```

Then run:

```
mpirun -n 28 $QE/pw.x -nk 2 < equil-nscf_333.in > equil-nscf_333.out
mpirun -n 28 $QE/pw.x -nk 2 < ZG-nscf_333_0.00K.in > ZG-nscf_333_0.00K.out
```

The calculations should take around 5 secs each. When the equilibrium calculation is finished, check that you obtain the VBM and CBM energies you would expect. To this aim type:

```
grep highest equil-nscf_333.out
```

It prints

```
highest occupied, lowest unoccupied level (ev): 6.2534 6.7435
```

Indeed we obtained the correct VBM and CBM energies.

To find the band gap renormalization from the ZG calculation type:

```
grep -41 "End of band" ZG-nscf_333_0.00K.out | tail -40
```

Can you find the energy of the VBM and CBM? You should find that there is a splitting of the originally threefold degenerated VBM of about 10 meV. The degeneracy splitting results from the ZG displacements in a finite size supercell and should reduce to zero for larger ZG simulation cells. This is a consequence of the harmonic approximation, where the symmetries of the system should be maintained upon thermal averaging (see thermal ellipsoids in Figure 3).

We note that to minimize the effect of band degeneracy splitting we exclude the modes in set *A*. In Phys. Rev. Research 2, 013357, (2020) we show that these modes dominate VBM splitting (see the spectral function in Figure 3), giving ΔE > 80 meV for a \(4\times4\times4\) ZG supercell. Hence, for finite size systems with degenerated band edges, it is suggested to keep `incl_qA = .false.`

. To deal with this finite size supercell artefact, we determine the VBM renormalization by taking the average of the three bands. Our analysis yields: E_VBM = (6.2740 + 6.2863 + 6.2920)/3 = 6.2841 eV. Thus the band gap at T = 0 K is E(0 K) = 6.7250 - 6.2841 = 0.4409 eV, and we can determine a zero-point renormalization E_ZPR = 0.49 - 0.4409 = 49 meV. This value is in good agreement with literature data [check, for example, J. Chem. Phys. 143, 102813 (2015) and New J. Phys. 20, 123008 (2018)].

In general, for obtaining reliable values for the zero-point renormalization, you should always check convergence with respect to the ZG supercell size (see the right plot in Figure 3).

Before going to the next step remove unnecessary files:

```
rm -r *wfc* *save
```

Now prepare a band structure unfolding calculation using the \(3\times3\times3\) ZG configuration. Create a new directory and copy the necessary input files:

```
mkdir band_structure_unfolding/; cd band_structure_unfolding/
cp ../../inputs/bands_unfold.in bands.in; cp ../../inputs/Si.pz-vbc.UPF .
cp -r ../../../exercise1/workdir/ZG-0.00K-si.save/ .
cp ../ZG-nscf_333_0.00K.in ZG-bands_333_0.00K.in
```

In **ZG-bands_333_0.00K.in** apply the following modifications:

```
1. calculation = 'bands'
2. Replace the following lines containing information about the K-points:
K_POINTS crystal
2
0.000000 0.000000 0.000000 1
0.000000 1.260000 1.260000 1
with
K_POINTS crystal_b
2
0.000000 0.000000 0.000000 35
0.000000 1.500000 1.500000 1
```

Here we choose to sample the Gamma-X path with 36 equally-spaced **K**-points. If you did not complete the previous step, then copy:

```
cp ../../inputs/ZG-bands_333_0.00K.in .
```

The file **ZG-bands_333_0.00K.in** should look like this:

```
&control ZG-bands_333_0.00K.in
calculation = 'bands'
restart_mode = 'from_scratch'
prefix = 'ZG-0.00K-si'
pseudo_dir = './'
outdir = './'
/
&system
ibrav = 0
nat = 54
ntyp = 1
ecutwfc = 20.00
nbnd = 135
nosym = .true.
occupations = 'fixed'
smearing = 'gaussian'
degauss = 0.0D+00
/
&electrons
diagonalization = 'david'
mixing_mode= 'plain'
mixing_beta = 0.70
conv_thr = 0.1D-06
/
ATOMIC_SPECIES
Si 28.085 Si.pz-vbc.UPF
K_POINTS crystal_b
2
0.000000 0.000000 0.000000 35
0.000000 1.500000 1.500000 1
CELL_PARAMETERS {angstrom}
-8.09641133 0.00000000 8.09641133
0.00000000 8.09641133 8.09641133
-8.09641133 8.09641133 0.00000000
ATOMIC_POSITIONS {angstrom}
...
```

We also show here the file **bands.in**:

```
-- bands.in
&bands
prefix = 'ZG-0.00K-si'
outdir = './'
filband = 'bands.dat'
lsym = .false.
dim1 = 3,
dim2 = 3,
dim3 = 3
/
```

The input variables in **bands.in** are the same as in a standard `bands.x`

calculation. The only difference is the flags `dimX`

used to specify the dimensions of the ZG supercell.

Now run a band structure unfolding calculation:

```
mpirun -n 56 $QE/pw.x -nk 28 < ZG-bands_333_0.00K.in > ZG-bands_333_0.00K.out
mpirun -n 28 $QE/bands_unfold.x < bands.in > bands.out
```

Note that parallelization over **k**-points in `bands_unfold.x`

is not implemented. The bands calculation should take around 50 secs and the unfolding around 1 sec. Now read the outputs of **bands01.dat** and **spectral_weights01.dat**. Can you guess what’s the meaning of the spectral weights for each band? The files **bands01.dat** and **spectral_weights01.dat**, containing the band energies, E_m (**K**, T), and spectral weights, P_m ( **K**, **k**, T), for all **K**-points, respectively. Now you have all the ingredients to evaluate the spectral function given by Eq. (61) of Phys. Rev. Research 2, 013357, (2020).

Now calculate the spectral function using the energies and spectral weights. To this aim convert first **bands01.dat** and **spectral_weights01.dat** into a different format using `plot_bands.x`

of QE:

```
cp ../../../inputs/energies.in .
cp ../../../inputs/spectral_weights.in .
$QE/bin/plotband.x spectral_weights01.dat < spectral_weights.in > spw.out
$QE/bin/plotband.x bands01.dat < energies.in > energies.out
sed -i '/^\s*$/d' spectral_weights.dat # remove empty lines
sed -i '/^\s*$/d' energies.dat # remove empty lines
paste energies.dat spectral_weights.dat > tmp
awk '{print $1,$2,$4}' tmp > energies_weights.dat; rm tmp
```

You have just generated **energies_weights.dat** file in a similar format to a “.gnu” file where the first column has the momentum-axis values, the second the energies, and the third the spectral weights. Now proceed with the calculation of the spectral function:

```
cp ../../../inputs/energies_weights.dat . # if you missed the step above
cp ../../../inputs/pp_spctrlfn.in .
```

The file `pp_spctrlfn.in`

takes the following input variables:

```
-- pp_spctrlfn.in
&input
flin = 'energies_weights.dat'
steps = 4725,
ksteps = 200,
esteps = 200,
kmin = 0,
kmax = 1.0,
emin = 3.00
emax = 9.00
flspfn = 'spectral_function.dat'
/
```

Here, (i) `steps`

is the number of rows in the input file `flin`

, (ii) `ksteps`

and `esteps`

define the resolution of the spectral function along the momentum-axis and energy-axis, (iii) `kmin`

- `kmax`

and `emin`

- `emax`

define the momentum and energy windows, and (iv) in `flspfn`

we provide the name of the output file.

Now calculate the spectral function using `pp_spctrlfn.x`

and plot:

```
mpirun -n 4 $QE/bin/pp_spctrlfn.x -nk 4 < pp_spctrlfn.in > pp_spctrlfn.out
cp spectral_function.dat ../../../gnuplot/. ; cd ../../../gnuplot/
cp ../inputs/bands.dat.gnu .
gnuplot gp_pm3d.p; evince bsu_si_0K.eps
```

We also copied **bands.dat.gnu** for comparing the spectral function with the band structure calculated with the atoms at their equilibrium positions (white dots in Figure 4). In the right part of Figure 4, we also reproduced the converged calculation from Phys. Rev. Research 2, 013357, (2020), employing an \(8\times8\times8\) ZG supercell, more **K**-points, and empty bands.

## Exercise 3a¶

In this exercise we will learn how to calculate temperature-dependent anharmonic phonons of Zr at *T=1188* K using \(2\times2\times2\) ZG supercells and the anharmonic special displacement method (A-SDM). A-SDM combines the special displacement method with the self-consistent phonon (SCP) theory; details for the theory related to this exercise please refer to Phys. Rev. B 108, 035155 (2023) and D. Hooton, LI. a new treatment of anharmonicity in lattice thermodynamics: I, London Edinburgh Philos. Mag. J. Sci. 46, 422 (1955). The main aim is to find the effective matrix of interatomic force constants (IFCs) *C* that best describes anharmonicity in the potential energy surface *U*. For this reason we minimize the trial free energy with respect to the matrix of IFCs which requires, essentially, to solve self-consistently the following equation:

In the following we are going to show how to solve the above equation self-consistently by performing a series of finite difference calculations.

First go to the directory `exercise3`

and copy the following input files in your `workdir`

:

```
cd ../../exercise3/; mkdir workdir; cd workdir
cp ../inputs/Zr.upf ../inputs/scf.in ../inputs/ph.in .
cp ../inputs/q2r.in ../inputs/matdyn.in .
```

Below we provide the **scf.in** file:

```
&control scf.in
calculation='scf'
restart_mode='from_scratch'
prefix = 'scf'
outdir='./'
pseudo_dir='./'
verbosity='high'
/
&system
ibrav=3
nat=1
ntyp=1
ecutwfc= 30
A = 3.51695
occupations = 'smearing'
smearing='gaussian',
degauss=0.005d0
/
&electrons
diagonalization='david'
mixing_mode='plain'
mixing_beta=0.7
conv_thr=1e-7
/
ATOMIC_SPECIES
Zr 91.224 Zr.upf
K_POINTS {automatic}
6 6 6 0 0 0
ATOMIC_POSITIONS (crystal)
Zr 0.0 0.0 0.0
```

Note that Gaussian smearing is used for the occupations since Zr is a metal. For speed up we use a small energy cutoff of 30 Ry.

Now run a self-consistent calculation for Zr in the `workdir`

:

```
mpirun -np 56 $QE/pw.x -nk 14 < scf.in > scf.out
```

Run a `ph.x`

calculation on a homogeneous \(3\times3\times3\) **q**-point grid using:

```
mpirun -np 56 $QE/ph.x -nk 14 < ph.in > ph.out
```

The **ph.in** is provided below:

```
-- ph.in
&inputph
amass(1)= 91.224,
prefix = 'scf'
outdir = './'
verbosity = 'high'
ldisp = .true.
fildyn = 'Zr.dyn'
tr2_ph = 1.0d-14
nq1 = 3
nq2 = 3
nq3 = 3
/
```

This will generate 4 **Zr.dynX** output files containing the dynamical matrix calculated for each irreducible **q**-point. The list of irreducible **q**-points is written in the **Zr.dyn0** file. Now run a `q2r.x`

calculation to obtain the IFCs file using:

```
mpirun -np 1 $QE/q2r.x < q2r.in > q2r.out
```

The **q2r.in** is also provided:

```
-- q2r.in
&input
fildyn='Zr.dyn', zasr='crystal', flfrc = 'Zr.333.fc'
/
```

This will generate **Zr.333.fc** which contains the harmonic IFCs.

Now run a `matdyn.x`

calculation to obtain the harmonic phonon dispersion using:

```
mpirun -np 1 $QE/matdyn.x < matdyn.in > matdyn.out
```

The **matdyn.in** is provided:

```
-- matdyn.in
&input
&input
asr='crystal', amass(1)= 91.224,
flfrc='Zr.333.fc', fldyn=' ', flfrq='Zr_harm.freq', fleig=' ', q_in_cryst_coord = .true.
q_in_band_form = .true.
/
5
0.00 0.00 0.00 100
0.50 0.50 0.50 100
0.75 0.25 -0.25 100
0.00 0.00 0.00 100
0.50 0.50 0.00 1
```

Then obtain the phonon dispersion using:

```
$QE/plotband.x
Input file > Zr_harm.freq
Reading 3 bands at 401 k-points
Range: -131.7185 166.3237eV Emin, Emax, [firstk, lastk] > -200 200
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 0.0000
high-symmetry point: 0.0000 0.0000 1.0000 x coordinate 1.0000
high-symmetry point: 0.5000 0.5000 0.5000 x coordinate 1.8660
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 2.7321
high-symmetry point: 0.0000 0.5000 0.5000 x coordinate 3.4392
output file (gnuplot/xmgr) > Zr_harm_ph_dispersion.xmgr
output file (ps) >
stopping ...
```

`matdyn.x`

will generate **Zr_harm.freq** which contains the harmonic phonon frequencies. Then we combine **Zr_harm.freq** with `plotband.x`

to obtain **Zr_harm_ph_dispersion.xmgr** in gnuplot/xmgr format for the harmonic phonon dispersion. To plot the phonon dispersion type:

```
cp Zr_harm_ph_dispersion.xmgr ../gnuplot/.
cd ../gnuplot/; gnuplot gp_coms_harm.p; evince Zr_harm_ph_dispersion.pdf
```

Is the matrix of IFCs positive definite (i.e. all phonon frequencies are positive)? Unlike Si phonon dispersion in Exercise 1 , the answer is no and the phonon dispersion exhibits severe instabilities. Here the harmonic approximation is not enough since anharmonic effects are important in Zr. To treat anharmonicity we are going to use the A-SDM.

The starting point of A-SDM for obtaining self-consistent phonons is a crucial step towards convergence and stability. This is because the SCP theory is designed to work with positive definite IFCs matrices. Although one could employ the harmonic phonons as a starting point (and neglect imaginary phonons), we strongly suggest the use of stable phonons as a starting point. For systems with a double well potential like Zr, well-defined dispersions can be obtained by exploring the ground state polymorphous structure, i.e a low-symmetry structure whose energy corresponds to a minimum of the potential energy surface [see also Phys. Rev. B 101, 155137 (2020)].

In the following we are going to show how to obtain the ground state polymorphous network of Zr and its stable phonons. If you have skipped the steps above, you can copy the harmonic IFCs file from `inputs`

into `workdir`

and run an A_ZG calculation:

```
cd ../workdir; cp ../inputs/Zr.333.fc .
cp ../inputs/ZG_1.in .
$QE/ZG.x < ZG_1.in > ZG_1.out
```

The input is provided here:

```
-- ZG_1.in
&input
asr='crystal', amass(1)= 91.224, flfrc='Zr.333.fc'
flscf = 'scf.in'
atm_zg(1) = 'Zr',
T = 0
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.2
ASDM = .true.
/
&A_ZG
poly = .true.
poly_fd_forces = .false.
/
```

Note

The **q**-point grid **nq1** * **nq2** * **nq3** used to obtain the initial dynamical matrices should not be necessarily the same with the supercell size (**dim1** * **dim2** * **dim3**). In particular, `dimX`

is independent of `nqX`

, since `ZG.x`

takes advantage of Fourier interpolation as implemented in `matdyn.x`

; thus any size of ZG configuration can be generated from **Zr.333.fc**.

The input format of `&input`

list in **ZG_1.in** is the same as we have seen in Exercise 1. We remark that here we are using a \(2\times2\times2\) supercell and therefore `incl_qA = .true.`

. The temperature is set initially to `T = 0`

, since we want to explore the polymorphous geometry. Line 7 contains the flag `ASDM`

which enables the A-SDM procedure and allows the code to read the input variables in the list `&A_ZG`

. At this step we ask the code by turning `poly = .true.`

to generate the input file for performing a coordinates’ optimization of the structure that will lead to the polymorphous geometry. Apart from the standard `ZG.x`

output files described in Exercise 1, the code now prints **ZG-relax.in**. This file contains displaced coordinates of the atoms away from their high-symmetry positions (that define the monomorphous network) in a \(2\times2\times2\) supercell using ZG displacements.

The file **ZG-relax.in** is provided here:

```
&control ZG-relax.in
calculation = relax
restart_mode = 'from_scratch'
prefix = 'ZG-relax'
pseudo_dir = './'
outdir = './'
tprnfor=.true.
disk_io = 'none'
forc_conv_thr = 1.0D-5
nstep = 300
/
&system
ibrav = 0
nat = 8
ntyp = 1
ecutwfc = 30.00
occupations = 'smearing'
smearing = 'gaussian'
degauss = 0.5D-02
/
&electrons
diagonalization = 'david'
mixing_mode= 'plain'
mixing_beta = 0.70
conv_thr = 0.1D-06
/
&ions
ion_dynamics = 'bfgs'
/
ATOMIC_SPECIES
Zr 91.224 Zr.upf
K_POINTS automatic
3 3 3 0 0 0
CELL_PARAMETERS (angstrom)
3.51695000 3.51695000 3.51695000
-3.51695000 3.51695000 3.51695000
-3.51695000 -3.51695000 3.51695000
ATOMIC_POSITIONS (angstrom)
...
```

Note

Cell paramaters, number of atoms, **k**-grid, and atomic coordinates have been modified automatically based on the supercell dimensions. The code will always generate the lattice information in angstroms.

Now run nuclei coordinates optimization for Zr with the nuclei clamped at their ZG harmonic coordinates (for *T = 0.00* K) in the \(2\times2\times2\) supercell.

```
mpirun -np 56 $QE/pw.x -nk 14 < ZG-relax.in > ZG-relax.out
```

The calculation should take around 3 minutes. Meanwhile you can check how the total energy evolves by typing `grep ! ZG-relax.out`

. Can you compare the total energy of the relaxed polymorphous structure with the one of the monomorphous structure (unit cell) calculation in **scf.out**? The polymorpous network is more stable by 27 meV/f.u.

Now run an A_ZG calculation using the same input file:

```
cp ../inputs/ZG_1.in ZG_2.in
$QE/ZG.x < ZG_2.in > ZG_2.out
```

The code now will search for the output file **ZG-relax.out**, read the relaxed final coordinates, and apply finite differences. This will generate **ZG-scf_poly_iter_00_XXXX.in** files corresponding to the scf inputs for evaluating the IFCs of the polymorphous network by finite differences. The default displacement on each atom is 0.01058 Å (i.e. 0.02 Bohr). One can change the default value by specifying a value in Å, using the flag `fd_displ`

under the namelist `&A_ZG`

.

Run scf calculations for the files **ZG-scf_poly_iter_00_XXXX.in** and save all scf outputs in the folder **fd_forces**:

```
file=ZG-scf_poly_iter_00
mkdir fd_forces/
for i in {0001..0048}; do
mkdir $i
cd $i
cp ../Zr.upf .
mv ../"$file"_"$i".in .
mpirun -np 56 $QE/pw.x -nk 14 < "$file"_"$i".in > "$file"_"$i".out
mv "$file"_"$i".out ../fd_forces/.
rm Zr.upf
cd ../
done
```

Each scf run should take around 3 seconds and the total calculation around 6 minutes. Now run an A_ZG calculation using the same input file but turn the flag `poly_fd_forces = .true.`

.

```
cp ../inputs/ZG_3.in ZG_3.in
$QE/ZG.x < ZG_3.in > ZG_3.out
```

**ZG_3.in** is also provided:

```
-- ZG_3.in
&input
asr='crystal', amass(1)= 91.224, flfrc='Zr.333.fc'
flscf = 'scf.in'
atm_zg(1) = 'Zr',
T = 0
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.3
ASDM = .true.
/
&A_ZG
poly = .true.
poly_fd_forces = .true.
incl_epsil = .false.
/
```

The code now will read all the forces from the files **ZG-scf_poly_iter_00_XXXX.out** in **fd_forces** and evaluate the IFCs of the polymorphous structure. Apart from the standard `ZG.x`

output files, the code now prints: (i) **FORCE_CONSTANTS_poly_iter_00** which contains the IFCs in the same format as phonopy, (ii) **FORCE_CONSTANTS_sym_poly_iter_00** which contains the symmetrized IFCs (i.e. after applying all symmetry operations of the crystal), and (iii) **poly_iter_00.fc** which contains the symmetrized IFCs in the same format as `q2r.x`

output, ready to be used for the next A-SDM iteration. You can check the symmetries found by the code using `grep "Sym." ZG_3.out`

or read **ZG_3.out** directly. For systems where the dielectric constant and Born effective charges are not zero, one may include this information in the **XXXX.fc** file by setting `incl_epsil = .true.`

. Note that for including the non-analytic contribution to the dynamical matrix one needs to set `na_ifc = .true.`

for the next ZG iteration run `&input`

list. This is because the IFCs are evaluated by finite differences instead of DFPT. The same flag is also necessary for obtaining the phonon dispersion with `matdyn.x`

; if not included the calculation will lead to erroneous results. An example is given for the case of PbTe in Exercise 3b

Run a `matdyn.x`

calculation to obtain the harmonic phonons of the Zr polymorphous network:

```
cp ../inputs/matdyn_poly.in .
mpirun -np 1 $QE/matdyn.x < matdyn_poly.in > matdyn_poly.out
$QE/plotband.x
```

`matdyn.x`

will generate **Zr_poly.freq** which contains the harmonic phonon frequencies of the Zr polymorphous network. Then combine **Zr_poly.freq** with `plotband.x`

to create the phonon dispersion file **Zr_poly_ph_dispersion.xmgr** as before. To plot the new phonon dispersion type:

```
cp Zr_poly_ph_dispersion.xmgr ../gnuplot/.
cd ../gnuplot/; gnuplot gp_coms_poly.p; evince Zr_poly_ph_dispersion.pdf
```

In this case, the matrix of IFCs is positive definite (i.e. all phonon frequencies are positive) and we are ready to account for anharmonicity via A-SDM using as a starting point the symmetrized IFCs of the Zr polymorphous network.

Proceed to run an A_ZG calculation using:

```
cd ../workdir; cp ../inputs/ZG_4.in ZG_4.in
$QE/ZG.x < ZG_4.in > ZG_4.out
```

The input file should look like this:

```
-- ZG_4.in
&input
asr='crystal', amass(1)= 91.224, flfrc='poly_iter_00.fc'
flscf = 'scf.in'
atm_zg(1) = 'Zr',
T = 1188
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.3
ASDM = .true.
fd = .true.
/
&A_ZG
iter_idx = 1
apply_fd = .true.
/
```

This is essentially the first iteration indicated by `iter_idx = 1`

. The file of IFCs is now `flfrc='poly_iter_00.fc'`

and the target temperature is *T= 1188* K. We also use `fd = .true.`

to specify that IFCs have been calculated using finite differences. `ZG.x`

will generate a ZG configuration on which the method of finite differences is applied after setting the flag `apply_fd = .true.`

. The scf files are **ZG-scf_1188.00K_iter_01_XXXX.in** to be used for the evaluation of the temperature-dependent IFCs of the configuration generated using ZG displacements for *T= 1188* K.

Run scf calculations for the files **ZG-scf_1188.00K_iter_01_XXXX.in** and save all scf outputs in the folder **fd_forces**. The script below is as before but with a new name for `file`

variable.

```
file=ZG-scf_1188.00K_iter_01
mkdir fd_forces/
for i in {0001..0048}; do
mkdir $i
cd $i
cp ../Zr.upf .
mv ../"$file"_"$i".in .
mpirun -np 56 $QE/pw.x -nk 14 < "$file"_"$i".in > "$file"_"$i".out
mv "$file"_"$i".out ../fd_forces/.
rm Zr.upf
cd ../
done
```

Each scf run should take around 3 seconds and the total calculation around 6 minutes.

Read the calculated forces with an A_ZG calculation using:

```
cp ../inputs/ZG_5.in ZG_5.in
$QE/ZG.x < ZG_5.in > ZG_5.out
```

The input file reads:

```
-- ZG_5.in
&input
asr='crystal', amass(1)= 91.224, flfrc='poly_iter_00.fc'
flscf = 'scf.in', atm_zg(1) = 'Zr',
T = 1188
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.3
ASDM = .true.
/
&A_ZG
iter_idx0 = 0, iter_idx = 1
apply_fd = .false.
read_fd_forces = .true.
mixing = .true.
incl_epsil = .false.
/
```

The flag `apply_fd`

is set to .false. since `apply_fd`

and `read_fd_forces`

cannot be both true. With the flag `read_fd_forces = .true.`

the code reads all the forces from files **ZG-scf_1188.00K_iter_01_XXXX.out** in the folder **fd_forces** and evaluates the IFCs of the ZG configuration. If one scf calculation is not successful, then the code will complain and print the error:

```
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Error in routine ZG (XXXX):
forces in file are missing ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
```

Here XXXX will be replaced by the code with the file index of the unsuccessful calculation. If all calculations run successfully, the code evaluates the IFCs using the finite differences formula and prints the following new files related to the first iteration: (i) **FORCE_CONSTANTS_1188.00K_iter_01**, (ii) **FORCE_CONSTANTS_sym_1188.00K_iter_01** and (iii) **1188.00K_iter_01.fc** which contains the symmetrized IFCs in the same `q2r.x`

output format, ready to be used for the next A-SDM iteration. Here, we set the flag `mixing = .true.`

, so that iterative mixing of the IFCs is performed. The input variables `iter_idx0`

and `iter_idx`

denote the initial and final index of iterations, respectively, for which mixing of IFCs is performed.

Now run `matdyn.x`

to obtain the temperature-dependent anharmonic phonons of the first iteration:

```
cp ../inputs/matdyn_01.in .
mpirun -np 1 $QE/matdyn.x < matdyn_01.in > matdyn_01.out
$QE/plotband.x
```

`matdyn.x`

will generate **Zr_01.freq** which contains the temperature-dependent anharmonic phonon frequencies of Zr at the first iteration. Then combine **Zr_01.freq** with `plotband.x`

to create the phonon dispersion file **Zr_01_ph_dispersion.xmgr** as before. To plot the new phonon dispersion and compare with the previous results type:

```
cp Zr_01_ph_dispersion.xmgr ../gnuplot/.
cd ../gnuplot/; gnuplot gp_coms_iter_01.p; evince Zr_01_ph_dispersion.pdf
```

As we can see the result from the first iteration follows a similar trend with the one calculated for the polymorphous network. In the following, we proceed to the second iteration and re-evaluate the phonon dispersion of Zr at *T= 1188* K starting from the IFCs of the previous iteration.

Run an A_ZG calculation using:

```
cd ../workdir; cp ../inputs/ZG_6.in ZG_6.in
$QE/ZG.x < ZG_6.in > ZG_6.out
```

The input file should look like this:

```
-- ZG_6.in
&input
asr='crystal', amass(1)= 91.224, flfrc='1188.00K_iter_01.fc'
flscf = 'scf.in'
atm_zg(1) = 'Zr',
T = 1188
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.35
ASDM = .true.
fd = .true.
/
&A_ZG
iter_idx = 2
apply_fd = .true.
/
```

This will generate a new ZG configuration and the scf files **ZG-scf_1188.00K_iter_02_XXXX.in** to be used for the evaluation of the temperature-dependent IFCs of the second iteration. In a similar spirit, run scf calculations for files **ZG-scf_1188.00K_iter_02_XXXX.in** and save all scf outputs in the folder **fd_forces**.

```
file=ZG-scf_1188.00K_iter_02
mkdir fd_forces/
for i in {0001..0048}; do
mkdir $i
cd $i
cp ../Zr.upf .
mv ../"$file"_"$i".in .
mpirun -np 56 $QE/pw.x -nk 14 < "$file"_"$i".in > "$file"_"$i".out
mv "$file"_"$i".out ../fd_forces/.
rm Zr.upf
cd ../
done
```

Now read the forces of iteration 2 using:

```
cp ../inputs/ZG_7.in ZG_7.in
$QE/ZG.x < ZG_7.in > ZG_7.out
```

The input file should look like this:

```
-- ZG_7.in
&input
asr='crystal', amass(1)= 91.224, flfrc='1188.00K_iter_01.fc'
flscf = 'scf.in', atm_zg(1) = 'Zr',
T = 1188
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.3
ASDM = .true.
fd = .true.
/
&A_ZG
iter_idx0 = 0, iter_idx = 2
apply_fd = .false.
read_fd_forces = .true.
mixing = .true.
incl_epsil = .false.
/
```

The code will read all the forces from the files **ZG-scf_1188.00K_iter_02_XXXX.in** in **fd_forces** and evaluate the IFCs of the new ZG configuration. The code now prints all files related to the second iteration: (i) **FORCE_CONSTANTS_1188.00K_iter_02**, (ii) **FORCE_CONSTANTS_sym_1188.00K_iter_02** and (iii) **1188.00K_iter_02.fc**. Mixing is performed between the current and all previous IFCs matrices as specified by the input variables `iter_idx0`

and `iter_idx`

. We found that convergence is achieved faster for the current system. Instead, one can perform mixing between consecutive iterations (current and previous iteration) by choosing appropriately `iter_idx0`

and `iter_idx`

. In any case the converged result should be similar.

Now run again a `matdyn.x`

calculation to obtain the temperature-dependent anharmonic phonons at the second iteration:

```
cp ../inputs/matdyn_02.in .
mpirun -np 1 $QE/matdyn.x < matdyn_02.in > matdyn_02.out
$QE/plotband.x
```

`matdyn.x`

will generate **Zr_02.freq** which contains the temperature-dependent anharmonic phonon frequencies of Zr at the second iteration. Then combine **Zr_02.freq** with `plotband.x`

to create the phonon dispersion file **Zr_02_ph_dispersion.xmgr** as before. To plot the new phonon dispersion and compare with the previous results type:

```
cp Zr_02_ph_dispersion.xmgr ../gnuplot/.
cd ../gnuplot/; gnuplot gp_coms_iter_02.p; evince Zr_02_ph_dispersion.pdf
```

As we can see the result at the second iteration is almost converged. One can repeat further iterative steps to check convergence. As a homework repeat the previous steps for the third and fourth iterations. The results should look like the following:

We can readily see that convergence is obtained at the third iteration; note that the result from the first iteration is still in excellent agreement with the converged result. As a homework, also calculate the phonons of Zr at *T=1188* K using \(3\times3\times3\) supercells.

**HINT**: Start the A-SDM procedure from the IFCs of Zr at *T=1188* K obtained at the third iteration with \(2\times2\times2\) supercells (i.e. using **1188.00K_iter_03.fc**). The result is provided in the plot below.

For the convergence of the phonon dispersion when larger supercells are used, please refer to the results in Phys. Rev. B 108, 035155 (2023). Therein, results are obtained using more accurate settings, i.e. denser **k**-grid and higher kinetic energy cutoff.

To this end we also discuss four final important points:

A quicker and more accurate way than checking the phonon dispersion is to plot the Frobenius norm (\(||.||\)) of the leading component of the IFCs matrix (\(\tilde{C}_{1 1 \alpha, 1 1 \alpha} (T)\)). This information is found from the first \(3\times3\) matrix in

FORCE_CONSTANTS_sym_1188.00K_iter_XX, where XX indicates the iteration. We note that convergence in IFCs guarantees essentially minimization of the system’s free energy. The result for a few more iterations is shown below. “iteration -1” corresponds to the harmonic \(||\tilde{C}_{1 1 \alpha, 1 1 \alpha}||\) and the difference is taken with respect to \(|| \tilde{C}_{1 1 \alpha, 1 1 \alpha} (T) ||\) of the last iteration. Note also that IFCs unit has been converted from Ry/Bohr \(^2\) to eV/ Å \(^2\).

A full self-consistent phonon theory also minimizes the trial free energy with respect to nuclear coordinates. One can use

`ZG.x`

which implements the Newton-Raphson method and find the nuclear coordinates for which \(\big< \partial U / \partial \tau_{\kappa \alpha} \big>_T = 0\). This requires to solve self-consistently the equation:

using at each iterative step the calculated \(\tilde{C}_{p \kappa \alpha, p \kappa \alpha} (T)\). In order to account for this feature in A-SDM we need to set the flag `update_equil = .true.`

when `read_fd_forces = .true.`

. If we use this flag for the case of bcc Zr (point group symmetry \(m\bar{3}m\)), then the code will give \(\tau^{\rm new}_{p \kappa \alpha} (T) = \tau_{p \kappa \alpha} (T)\) (i.e. no change of the nuclear coordinates) at every iterative step. This comes as no surprise since symmetries are enforced and atoms in bcc Zr occupy only special Wyckoff positions. However, for systems containing atoms at general Wyckoff positions, the equation above will give new thermal equilibrium positions at each temperature that contribute to the minimization of the free energy.

Obtaining the polymorphous structure is not a necessary step. We found that this is a very useful step for obtaining positive definite initial IFCs matrix for systems with a double well potential. If the harmonic approximation gives a positive definite initial IFCs matrix, then one could start the iterative procedure without obtaining the polymorphous network. See for example Exercise 3b.

After computing temperature-dependent anharmonic phonons using A-SDM, one could then do a ZG calculation using the new configuration and evaluate anharmonic electron-phonon properties.

## Exercise 3b¶

In this exercise we will learn how to calculate temperature-dependent anharmonic phonons of PbTe at *T=300* K using \(2\times2\times2\) ZG supercells and the A-SDM. Unlike Zr, seen in Exercise 3a, PbTe is a polar material and a correction to the dynamical matrix arising from long-range dipole-dipole interactions is essential. Since the A-SDM relies on the computation of IFCs by finite differences, we employ the mixed space approach described in J. Phys.: Condens. Matter 22 202201 (2010), as implemented in `matdyn.x`

. The key flags here to be used are `fd = .true.`

, `na_ifc = .true.`

, and `incl_epsil = .true.`

. Furthermore, in this example, the harmonic IFCs matrix is positive definite and therefore the computation of the polymorphous structure is not required to initialize the A-SDM procedure.

First go to the directory `exercise3`

and create `workdirb`

:

```
cd ../../exercise3/; mkdir workdirb; cd workdirb
```

Below we provide the **scf.in** file:
::
– scf.in

- &control
calculation=’scf’ restart_mode=’from_scratch’ prefix = ‘scf’ outdir=’./’ pseudo_dir=’/users/zacharia/nc-sr-04_pbesol_standard_upf/’ tstress=.true. tprnfor=.true.

/ &system

ibrav=2 space_group = 225 nat=2 ntyp=2 ecutwfc= 100 A = 6.4440

/ &electrons

diagonalization=’david’ mixing_mode=’plain’ mixing_beta=0.7 conv_thr=1e-9

/ ATOMIC_SPECIES Pb 207.2 Pb.upf Te 127.60 Te.upf K_POINTS {automatic} 8 8 8 0 0 0 ATOMIC_POSITIONS (crystal_sg)

Pb 4b Te 4a

Note

Here we provide the space group of the crystal, as given in the International Tables of Crystallography A (ITA), and the atomic coordinates in Wyckoff positions. This is an alternative way to specify the geometry of the system in Quantum Espresso which ensures that the symmetry of the crystal is respected.

Here we use the PBEsol approximation; the pseudopotential files can be downloaded from the Pseudo-dojo Library.

We also provide the **ph.in** file:

```
PbTe ph.in
&inputph
prefix = 'scf',
amass(1) = 207.2
amass(2) = 127.60
outdir = './'
verbosity = 'high'
ldisp = .true.
fildyn = 'PbTe.dyn'
tr2_ph = 1.0d-14
nq1 = 2
nq2 = 2
nq3 = 2
/
```

**q2r.in** file:

```
-- q2r.in
&input
fildyn='PbTe.dyn', zasr='crystal', flfrc = 'PbTe.222.fc'
/
```

**matdyn.x** file

```
-- matdyn.in
&input
asr='crystal',
flfrc='PbTe.222.fc', fldyn=' ', flfrq='PbTe_harm.freq', fleig=' ', q_in_cryst_coord = .true.
q_in_band_form = .true.
/
3
0 0.5 0.5 100
0 0.0 0.0 100
0 0.5 0.5 100
```

Now run a self-consistent and a DFPT calculation for PbTe and obtain the IFCs and harmonic phonon frequencies in the `workdirb`

:

```
mpirun -np 58 $QE/pw.x -nk 29 < scf.in > scf.out
mpirun -np 58 $QE/ph.x -nk 29 < ph.in > ph.out
mpirun -np 1 $QE/q2r.x < q2r.in > q2r.out
mpirun -np 1 $QE/matdyn.x < matdyn.in > matdyn.out
```

After the calculation is completed use `plotband.x`

of Quantum Espresso to obtain the phonon dispersion data as an a xmgr file.

```
$QE/plotband.x
Input file > PbTe_harm.freq
Reading 6 bands at 201 k-points
Range: 0.0000 107.8059eV Emin, Emax, [firstk, lastk] > 0 110
high-symmetry point: 0.0000 1.0000 0.0000 x coordinate 0.0000
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 1.0000
high-symmetry point: 0.0000 1.0000 0.0000 x coordinate 2.0000
output file (gnuplot/xmgr) > PbTe_harm_ph_dispersion.xmgr
output file (ps) >
stopping ...
```

Now plot the phonon dispersion which should look like the one shown in **Figure 12**.

The matrix of IFCs is positive definite and can be used for initializing the A-SDM procedure. Perform now a `ZG.x`

calculation to generate a \(2\times2\times2\) ZG configuration on which the method of finite differences is applied after setting the flag `apply_fd = .true.`

. The scf files are **ZG-scf_300.00K_iter_01_XXXX.in** to be used for the evaluation of the temperature-dependent IFCs of the configuration generated using ZG displacements for *T= 300* K.

Now run:

```
$QE/ZG.x < ZG_1.in > ZG_1.out
```

**ZG_1.in** is provided here:

```
-- ZG_1.in
&input
asr='crystal',
flscf = 'scf.in'
flfrc='PbTe.222.fc'
atm_zg(1) = 'Pb',
atm_zg(2) = 'Te',
T = 300
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.4
ASDM = .true.
/
&A_ZG
iter_idx = 1
apply_fd = .true.
/
```

Run scf calculations for the files **ZG-scf_300.00K_iter_01_XXXX.in** and save all scf outputs in the folder **fd_forces**:

```
file=ZG-scf_300.00K_iter_01
mkdir fd_forces/
for i in {0001..0096}; do
mkdir $i
cd $i
cp ../Pb.upf .
cp ../Te.upf .
mv ../"$file"_"$i".in .
mpirun -np 36 $QE/pw.x -nk 18 < "$file"_"$i".in > "$file"_"$i".out
mv "$file"_"$i".out ../fd_forces/.
rm Pb.upf Te.upf
cd ../
done
```

Read the calculated forces with an A_ZG calculation using:

```
$QE/ZG.x < ZG_2.in > ZG_2.out
```

The input file reads:

```
-- ZG_2.in
&input
asr='crystal',
flscf = 'scf.in'
flfrc='PbTe.222.fc'
atm_zg(1) = 'Pb',
atm_zg(2) = 'Te',
T = 300
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.6
ASDM = .true.
/
&A_ZG
iter_idx = 1
apply_fd = .false.
read_fd_forces = .true.
mixing = .false.
incl_epsil = .true.
/
```

The flag `apply_fd`

is set to .false. since `apply_fd`

and `read_fd_forces`

cannot be both true. With the flag `read_fd_forces = .true.`

the code reads all the forces from files **ZG-scf_300.00K_iter_01_XXXX.out** in the folder **fd_forces** and evaluates the IFCs of the ZG configuration. If all calculations run successfully, the code evaluates the IFCs using the finite differences formula and prints the following new files related to the first iteration: (i) **FORCE_CONSTANTS_300.00K_iter_01**, (ii) **FORCE_CONSTANTS_sym_300.00K_iter_01** and (iii) **300.00K_iter_01.fc** which contains the symmetrized IFCs in the same `q2r.x`

output format, ready to be used for the next A-SDM iteration. Here, we set the flag `incl_epsil = .true.`

, to write the information of the dielectric constant and Born Effective charges in **300.00K_iter_01.fc**. Those should look like this:

```
T
31.235827868894 0.000000000000 -0.000000000000
0.000000000000 31.235827868894 0.000000000000
-0.000000000000 0.000000000000 31.235827868893
1
5.7365775 0.0000000 0.0000000
0.0000000 5.7365775 -0.0000000
-0.0000000 -0.0000000 5.7365775
2
-5.7365775 -0.0000000 -0.0000000
-0.0000000 -5.7365775 0.0000000
-0.0000000 0.0000000 -5.7365775
```

Now run `matdyn.x`

to obtain the temperature-dependent anharmonic phonons of the first iteration:

```
mpirun -np 1 $QE/matdyn.x < matdyn_01.in > matdyn_01.out
$QE/plotband.x
```

**matdyn_01.in** is provided here:

```
-- matdyn_01.in
&input
asr='crystal',
flfrc='300.00K_iter_01.fc', fldyn=' ', flfrq='PbTe_01.freq', fleig=' ', q_in_cryst_coord = .true.
q_in_band_form = .true., fd = .true., na_ifc =.true.
/
3
0 0.5 0.5 100
0 0.0 0.0 100
0 0.5 0.5 100
```

Note

Here we use `fd = .true.`

and `na_ifc = .true.`

since finite differences are used for computing the A-SDM IFCs.

`matdyn.x`

will generate **PbTe_01.freq** which contains the temperature-dependent anharmonic phonon frequencies of PbTe at the first iteration. Then **PbTe_01.freq** is combined with `plotband.x`

to create the phonon dispersion file **PbTe_01_ph_dispersion.xmgr** as before. Now plot the new phonon dispersion and compare with the previous results as shown in **Figure 13**.

Now proceed to the second and third iterations using the same procedure as before for calculating the forces and also apply mixing. The input files are provided below:

**ZG_3.in**:

```
-- ZG_3.in
&input
asr='crystal',
flscf = 'scf.in'
flfrc='300.00K_iter_01.fc'
atm_zg(1) = 'Pb',
atm_zg(2) = 'Te',
T = 300
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.6
ASDM = .true.
na_ifc = .true., fd =.true.
/
&A_ZG
iter_idx = 2
apply_fd = .true.
/
```

**ZG_4.in**:

```
-- ZG_4.in
&input
asr='crystal',
flscf = 'scf.in'
flfrc='300.00K_iter_01.fc'
atm_zg(1) = 'Pb',
atm_zg(2) = 'Te',
T = 300
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.6
ASDM = .true.
na_ifc = .true., fd =.true.
/
&A_ZG
iter_idx0 =1, iter_idx = 2
apply_fd = .false.
read_fd_forces = .true.
mixing = .true.
incl_epsil = .true.
/
```

**matdyn_02.in**

```
-- matdyn_02.in
&input
asr='crystal',
flfrc='300.00K_iter_02.fc', fldyn=' ', flfrq='PbTe_02.freq', fleig=' ', q_in_cryst_coord = .true.
q_in_band_form = .true., na_ifc =.true., fd =.true.
/
3
0 0.5 0.5 100
0 0.0 0.0 100
0 0.5 0.5 100
```

**ZG_5.in**:

```
-- ZG_5.in
&input
asr='crystal',
flscf = 'scf.in'
flfrc='300.00K_iter_02.fc'
atm_zg(1) = 'Pb',
atm_zg(2) = 'Te',
T = 300
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.6
ASDM = .true.
na_ifc = .true., fd =.true.
/
&A_ZG
iter_idx = 3
apply_fd = .true.
/
```

**ZG_6.in**:

```
-- ZG_6.in
&input
asr='crystal',
flscf = 'scf.in'
flfrc='300.00K_iter_02.fc'
atm_zg(1) = 'Pb',
atm_zg(2) = 'Te',
T = 300
dim1 = 2, dim2 =2, dim3 =2, incl_qA = .true.
compute_error = .true., synch = .true., error_thresh = 0.6
ASDM = .true.
na_ifc = .true., fd =.true.
/
&A_ZG
iter_idx0 =1, iter_idx = 3
apply_fd = .false.
read_fd_forces = .true.
mixing = .true.
incl_epsil = .true.
/
```

**matdyn_03.in**

```
-- matdyn_03.in
&input
asr='crystal',
flfrc='300.00K_iter_03.fc', fldyn=' ', flfrq='PbTe_03.freq', fleig=' ', q_in_cryst_coord = .true.
q_in_band_form = .true., na_ifc =.true., fd =.true.
/
3
0 0.5 0.5 100
0 0.0 0.0 100
0 0.5 0.5 100
```

After you finish iteration 2 and 3 plot all phonon dispersions and check for convergence as shown in **Figure 14**:

Can you repeat now the calculation for *T=500* K and *T=700* K (first iteration only)? The result should look like as in **Figure 15**:

## Exercise 4¶

In this exercise we will learn how to calculate the phonon-assisted optical spectra using the example of silicon and the \(3\times3\times3\) ZG supercell for *T=0* K. To obtain temperature-dependent optical spectra we will evaluate the optical matrix elements in the independent-particle approximation using a routine very similar to `epsilon.x`

of Quantum Espresso. We emphasize that this methodology can be extended to higher-level theories of optical absorption, like RPA, BSE, etc … The present approach of calculating phonon-assisted optical spectra is based on Phys. Rev. B 94, 075125, (2016).

The calculation of phonon-assisted optical spectra using ZG configurations has some drawbacks and advantages when compared to the perturbative methodology implemented in `epw.x`

.

First go to the directory `exercise4`

, create your working directories and copy the input files:

```
cd ../../exercise4/; mkdir workdir; cd workdir; mkdir equil; mkdir ZG
cd equil
cp ../../inputs/equil/K_list.in .
cp ../../inputs/equil/epsilon.in .
cp ../../inputs/equil/epsi_av.sh .
cp -r ../../../exercise1/workdir/equil-si.save/ .
cp ../../../exercise2/workdir/equil-nscf_333.in .
```

The file **K_list.in** contains the crystal coordinates of 200 randomly generated **K**-points. We use random **K**-points, instead of a uniform grid, in order to speed up convergence. In **equil-nscf_333.in** apply the following changes:

```
1. Add nosym = .true. below nbnd = 135.
2. Delete the information for the K-points:
K_POINTS crystal
2
0.000000 0.000000 0.000000 1
0.000000 1.260000 1.260000 1
3. At the bottom of the file add the following:
K_POINTS crystal
30
4. Copy and paste the first 30 K-points from file K_list.in.
```

Note

Although we use the equilibrium structure we set `nosym = .true.`

since random **K**-points are employed.

If you did not complete Exercise 1 and Exercise 2, you can copy:

```
cp ../../inputs/equil/equil-nscf_333.in .
cp -r ../../../exercise2/inputs/equil-si.save/ .
```

The file **equil-nscf_333.in** should look like this:

```
&control equil-nscf_333.in
calculation = 'nscf'
restart_mode = 'from_scratch'
prefix = 'equil-si'
pseudo_dir = './'
outdir = './'
/
&system
ibrav = 0
nat = 54
ntyp = 1
ecutwfc = 20.00
nbnd = 135
nosym = .true.
occupations = 'fixed'
smearing = 'gaussian'
degauss = 0.0D+00
/
&electrons
diagonalization = 'david'
mixing_mode= 'plain'
mixing_beta = 0.70
conv_thr = 0.1D-06
/
ATOMIC_SPECIES
Si 28.085 Si.pz-vbc.UPF
CELL_PARAMETERS {angstrom}
-8.09641133 0.00000000 8.09641133
0.00000000 8.09641133 8.09641133
-8.09641133 8.09641133 0.00000000
ATOMIC_POSITIONS {angstrom}
Si 0.00000000 0.00000000 0.00000000
Si -1.34940189 1.34940189 1.34940189
Si -2.69880378 0.00000000 2.69880378
...
K_POINTS crystal
30
0.236650 0.022946 0.917156 1.0
0.392728 0.655175 0.471114 1.0
...
```

We also show here the file **epsilon.in**:

```
-- epsilon.in
&inputpp
outdir = './',
prefix = 'equil-si',
calculation = 'eps'
/
&energy_grid
smeartype = 'gauss',
intersmear = 0.03d0,
wmax = 4.5d0,
wmin = 0.2d0,
nw = 600,
shift = 0.0d0,
/
```

The input variables in **epsilon.in** are the same as in a standard `epsilon.x`

calculation.

Run an optical spectrum calculation for silicon using the \(3\times3\times3\) equilibrium supercell and 30 random **K**-points. The aim is to calculate the imaginary part of the dielectric function [\(\varepsilon_2(\omega)\)] using Eq. (2) of Phys. Rev. B 94, 075125, (2016). To calculate the dielectric function we use `epsilon_Gaus.x`

, which is a slightly modified routine of `epsilon.x`

, that replaces the delta function with a Gaussian instead of a Lorentzian function. We choose to apply Gaussian broadening in order to minimize numerical artefacts at the absorption onset.

We calculate \(\varepsilon_2(\omega)\) using:

```
mpirun -n 56 $QE/pw.x -nk 28 < equil-nscf_333.in > equil-nscf_333.out
mpirun -n 28 $QE/epsilon_Gaus.x < epsilon.in > epsilon.out
```

The full calculation should take around 40 secs. The output files are (i) **epsi_equil-si.dat** containing \(\varepsilon_2(\omega)\), (ii) **epsr_equil-si.dat** containing the real part of the dielectric function \(\varepsilon_1(\omega)\), and (iii) **eels_equil-si.dat** containing the electron energy loss spectrum; all calculated for the equilibrium structure. In each file the first column is the energy grid, and the rest three represent the observable along the three Cartesian directions.

Take the isotropic average of the dielectric function and copy the output file in the `gnuplot`

directory using:

```
awk '{print $1,($2+$3+$4)/3}' epsi_equil-si.dat > epsi_si_333_equil_30_av.dat
cp epsi_si_333_equil_30_av.dat ../../gnuplot/.
```

Run an optical spectra calculation for silicon using the \(3\times3\times3\) ZG supercell and 30 random **K**-points. We will essentially repeat all steps performed for calculating the optical spectra of the \(3\times3\times3\) equilibrium supercell. First go to the `ZG`

directory and copy the input files:

```
cd ../ZG/
cp ../../inputs/ZG/K_list.in .
cp ../../inputs/ZG/epsilon.in .
cp ../../inputs/ZG/epsi_av.sh .
cp -r ../../../exercise1/workdir/ZG-0.00K-si.save/ .
cp ../../../exercise2/workdir/ZG-nscf_333_0.00K.in .
```

In **ZG-nscf_333_0.00K.in** apply the following changes:

```
1. Delete the information for the K-points:
K_POINTS crystal
2
0.000000 0.000000 0.000000 1
0.000000 1.260000 1.260000 1
2. At the bottom of the file add the following:
K_POINTS crystal
30
3. Copy and paste the first 30 K-points from file K_list.in.
```

If you did not complete Exercise 1 and Exercise 2 you can copy **ZG-0.00K.save** and
**ZG-nscf_333_0.00K.in** from `inputs`

directory, i.e.:

```
cp ../../inputs/ZG/ZG-nscf_333_0.00K.in ; cp -r ../../../exercise2/inputs/ZG-0.00K-si.save/ .
```

The file **ZG-nscf_333_0.00K.in** should look like this:

```
&control ZG-nscf_333_0.00K.in
calculation = 'nscf'
restart_mode = 'from_scratch'
prefix = 'ZG-0.00K-si'
pseudo_dir = './'
outdir = './'
/
&system
ibrav = 0
nat = 54
ntyp = 1
ecutwfc = 20.00
nbnd = 135
nosym = .true.
occupations = 'fixed'
smearing = 'gaussian'
degauss = 0.0D+00
/
&electrons
diagonalization = 'david'
mixing_mode= 'plain'
mixing_beta = 0.70
conv_thr = 0.1D-06
/
ATOMIC_SPECIES
Si 28.085 Si.pz-vbc.UPF
CELL_PARAMETERS {angstrom}
-8.09641133 0.00000000 8.09641133
0.00000000 8.09641133 8.09641133
-8.09641133 8.09641133 0.00000000
ATOMIC_POSITIONS {angstrom}
Si -0.02610870 -0.00524574 0.02749582
Si -1.36329713 1.31619849 1.42481245
Si -2.75218400 -0.00549830 2.71681282
...
K_POINTS crystal
30
0.236650 0.022946 0.917156 1.0
0.392728 0.655175 0.471114 1.0
...
```

**epsilon.in** is the same with the one used before, but with `prefix = 'ZG-0.00K-si'`

. Now run the following to calculate \(\varepsilon_2(\omega)\) at *T=0.00* K:

```
mpirun -n 56 $QE/pw.x -nk 28 < ZG-nscf_333_0.00K.in > ZG-nscf_333_0.00K.out
mpirun -n 28 $QE/epsilon_Gaus.x < epsilon.in > epsilon.out
```

The full calculation should take around 40 secs. The output files are (i) **epsi_ZG-0.00K-si.dat** containing \(\varepsilon_2(\omega)\), (ii) **epsr_ZG-0.00K-si.dat** containing the real part of the dielectric function \(\varepsilon_1(\omega)\), and (iii) **eels_ZG-0.00K-si.dat** containing the electron energy loss spectrum; all calculated for the equilibrium structure. In each file the first column is the energy grid, and the rest three represent the observable along the three Cartesian directions.

Take the isotropic average of the dielectric function and plot the output file in the **gnuplot** directory using:

```
awk '{print $1,($2+$3+$4)/3}' epsi_ZG-0.00K-si.dat > epsi_si_333_ZG_30_av.dat
cp epsi_si_333_ZG_30_av.dat ../../gnuplot/.
cd ../../gnuplot/; gnuplot gp_coms.p; evince ZG_spectra.eps
gnuplot gp_coms_logscale.p; evince ZG_spectra_logscale.eps
```

Below, we also plot \(\varepsilon_2(\omega)\) as calculated with more **K**-points to obtain convergence.

What are the differences between the ZG and equilibrium spectra in Figure 16 **b** ? Are phonon-assisted transitions captured as expected? Phys. Rev. B 94, 075125, (2016) reports the convergence of the spectra with respect to the ZG supercell size and the agreement with experimental data for the absorption coefficient.

## Exercise 5¶

In this exercise we will discuss the physical meaning of the ZG configuration and learn how to calculate phonon-induced diffuse (or inelastic) scattering patterns using graphene. This is an example of an observable which allows us to computationally reach the thermodynamic limit and assess multi-phonon processes. More details about the theory of phonon-induced inelastic scattering and its connection to the special displacement method can be found in Phys. Rev. Lett. 127, 207401 (2021) and Phys. Rev. B 104, 205109 (2021).

We will evaluate the exact expression for the all-phonon inelastic scattering intensity and compare it with the ZG scattering intensity, i.e. when the scatterers (nuclei) are defined by the ZG displacements. Go to the directory **exercise5**, create your working directories, and copy the input files:

```
cd ../../exercise5/; mkdir workdir; cd workdir; mkdir exact; mkdir ZG
cd exact
cp ../../inputs/exact/graphene.881.fc .
cp ../../inputs/exact/disca.in .
```

The file `graphene.881.fc`

contains the IFCs of graphene calculated using an \(8\times8\times1\) **q**-grid.

Run a diffuse scattering calculation using `disca.x`

and the input:

```
-- disca.in
&input
asr = 'crystal', amass(1) = 12.011, flfrc = 'graphene.881.fc', atm_zg(1) = 'C',
dim1 = 40, dim2 = 40, dim3 = 1,
T = 300,
atmsf_a(1,1) = 0.1361,
atmsf_a(1,2) = 0.5482,
atmsf_a(1,3) = 1.2266,
atmsf_a(1,4) = 0.5971,
atmsf_a(1,5) = 0.0,
atmsf_b(1,1) = 0.3731,
atmsf_b(1,2) = 3.2814,
atmsf_b(1,3) = 13.0456,
atmsf_b(1,4) = 41.0202,
atmsf_b(1,5) = 0.0,
zero_one_phonon = .true., full_phonon = .true.
nks1 = 0, nks2 = 0, nks3 = 0, nksf1 = 5, nksf2 = 5, nksf3 = 1,
plane_val = 0.d0, plane_dir = 3
qstart = 1, qfinal = 40000,
/
&pp_disca
nrots = 6,
kres1 = 250,
kres2 = 250,
kmin = -10,
kmax = 10,
col1 = 1,
col2 = 2,
Np = 1600
/
```

and the command:

```
mpirun -np 56 $QE/bin/disca.x -nk 56 < disca.in > disca.out
```

The input format of `disca.in`

is the same as `ZG_333.in`

. Line 3 contains input variables similar to those in matdyn.in; in line 4 we provide the dimensions `dimX`

of the `q`

-grid used to sample the Brillouin zone; in line 5 we provide the temperature `T`

in Kelvin; lines 6-15 contain the parameters of the atomic scattering factor `atmsf_a`

and `atmsf_b`

from Micron 30, 625–648, (1999) in line 16 we specify whether we want to compute the one-phonon (`zero_one_phonon`

) and/or all-phonon (`full_phonon`

) structure factor; in line 17 `nksX`

and `nksfX`

define the range of reciprocal lattice vectors (in crystal coordinates) and are used to determine the extent of the calculated scattering vectors **Q**-points); in line 18 we provide `plane_val`

that defines the plane along which the structure factor is calculated (in units of 2pi/alat) and `plane_dir`

defines the Cartesian direction perpendicular to the plane (where 1 –> x, 2 –> y, and 3 –> z); and in line 19 we specify the range of **Q**-points (from `qstart`

to `qfinal`

) for which the structure factor is calculated.

In the namelist `&pp_disca`

we have: (i) `nrots`

which is the number of rotations to be applied on the all-phonon scattering intensity to
obtain the complete map. `nrots`

is based on the space group of each system. This step is important to avoid the extra effort of calculating the scattering intensity for wavevectors connected by rotation symmetry (see also plots at the end of the exercise). (ii) `kres1`

and `kres2`

define the resolution of the scattering intensity along the chosen Cartesian axes, (iii) `col1`

and `col2`

are integers defining two of the Cartesian directions of the scattering vectors (where 1 –> x, 2 –> y, and 3 –> z), (iv) (`kmin`

- `kmax`

) define the 2D momentum window in reciprocal space, and (v) `Np`

is the number of reduced wavevectors (**q**-points) used to sample each Brillouin zone.

More details about all input flags are provided in this link.

The standard outputs from a `disca.x`

run are: **strf_one-ph_broad.dat** and **strf_all-ph_broad.dat** containing the zero+one-phonon, and all-phonon scattering intensities, respectively, after convolution is applied. These files have three columns: the first two represent the Cartesian coordinates of the scattering vectors and the third column is the scattering intensity.

One can also print for each case the mode resolved and atom resolved scattering intensities by setting `mode_resolved = .true.`

and/or `atom_resolved = .true.`

.

For printing the raw data, i.e. before convolution is applied, one can use `print_raw_data = .true.`

to obtain the following output data: **Bragg_scattering.dat**, **strf_one-phonon.dat**, **strf_q_nu_one-phonon.dat**, and **strf_all-phonon.dat**, containing the Bragg, mode (or branch) resolved, one-phonon, and all-phonon scattering intensities, respectively. These files have four columns: the first three represent the Cartesian coordinates of the scattering vectors (Q_x, Q_y, and Q_z) and the fourth column is the scattering intensity. `disca.x`

also prints the **Q** vectors in crystal coordinates in **qpts_strf.dat** (to be used for computing the ZG scattering intensity).

Now copy the output file **strf_all-ph_broad.dat** in the `gnuplot`

directory:

```
cp strf_all-ph_broad.dat ../../gnuplot/exact/
```

Run a `ZG.x`

calculation that allows to compute the structure factor. Go to the `ZG`

directory and copy the input files:
::https://doi.org/10.1103/PhysRevLett.127.207401

cd ../ZG; cp ../../inputs/ZG/graphene.881.fc .; cp ../../inputs/ZG/ZG_strf.in . cp ../exact/qpts_strf.dat .

Note that we also copied **qpts_strf.dat** containing the **Q**-points generated by `disca.x`

.

The `ZG_strf.in`

should look like this:

```
-- ZG_strf.in
&input
asr='crystal', amass(1)= 12.011, flfrc= 'graphene.881.fc', atm_zg(1) = 'C'
dim1 = 40, dim2 = 40, dim3 = 1,
T = 300,
synch = .true., compute_error = .false.
incl_qA = .true.,
ZG_strf = .true.
/
&strf_ZG
qpts_strf = 40000
atmsf_a(1,1) = 0.1361, atmsf_a(1,2) = 0.5482, atmsf_a(1,3) = 1.2266, atmsf_a(1,4) = 0.5971,
atmsf_a(1,5) = 0.0,
atmsf_b(1,1) = 0.3731, atmsf_b(1,2) = 3.2814, atmsf_b(1,3) = 13.0456, atmsf_b(1,4) = 41.0202,
atmsf_b(1,5) = 0.0,
nrots = 6,
kres1 = 250,
kres2 = 250,
kmin = -10,
kmax = 10,
col1 = 1,
col2 = 2,
Np = 1600
/
```

Now run:

```
mpirun -np 56 $QE/bin/ZG.x -nk 56 < ZG_strf.in > ZG_strf.out
```

The input variables in the first 3 lines of `ZG_strf.in`

are the same as in `disca.in`

(but strictly speaking `dimX`

, here, is for the supercell size dimensions); in line 4 we ask the code to synchronize the modes `synch = .true.`

, but not to compute the minimization function, as we have seen in Exercise 1, i.e. `compute_error = .false.`

. We made this choice to demonstrate that the order of choosing the unique set of signs is not important as we approach the thermodynamic limit; in line 5 we make the choice to include **q**-points in set *A* by setting `incl_qA = .true`

(at the thermodynamic limit should make no difference); in line 6 we ask the code to compute the structure factor `ZG_strf = .true.`

.

The namelist `&strf_ZG`

is as `&pp_disca`

seen before, but here we include the parameters of the atomic scattering factor `atmsf_a`

and `atmsf_b`

. We also ask the code to compute the map for 40000 **Q**-points `qpts_strf = 40000`

from file `qpts_strf.dat`

.

More details about all input flags are provided in this link.

The calculation should take around 10 secs. The important outputs from this `ZG.x`

run are: **ZG-configuration_0.050.dat** and **strf_ZG_broad.dat**, containing the collection of ZG atomic coordinates (as before) and the associated scattering intensity, respectively. File **structure_factor_ZG_raw.dat** contains the raw data and has four columns: the first three represent the Cartesian coordinates of the scattering vectors and the fourth column is the scattering intensity.

Now copy the output file **strf_ZG_broad.dat** in the `gnuplot`

directory and plot your ZG and exact data:

```
cp strf_ZG_broad.dat ../../gnuplot/ZG
cd ../../gnuplot/ZG; gnuplot gp_pm3d.p; evince strf_ZG_40_40.eps
cd ../exact; gnuplot gp_pm3d.p; evince strf_exact_40_40.eps
```

Your results should look like the plots in Figure 17 (left column). We also add plots of the exact and ZG all-phonon structure maps using a \(100\times100\) and \(200\times200\) **q**-grid for Brillouin sampling. What happens to the ZG structure factor map as we approach the thermodynamic limit? Indeed, the choice of the order of signs is no longer important and the error in the ZG calculation is minimized at the thermodynamic limit. Our calculations also demonstrate the physical significance of the ZG configuration as the best collection of scatterers that reproduce the phonon-induced inelastic scattering.

## Exercise 6¶

In this exercise we will learn how to perform phonon unfolding and compute the phonon spectral function. This technique is useful to identify the effect on the phonons in the fundamental Brillouin zone when supercells are employed to accommodate disorder, defects, or any other perturbation that breaks the translational symmetry of the lattice. First, as a proof of concept, we will perform the perfect unfolding of the phonons calculated for a \(2\times2\times2\) silicon supercell; then we will show how to calculate the effect on the phonons when a \(2\times2\times2\) silicon supercell is doped by a phosphorous (P) atom.

Go to `exercise6`

, create **workdir**, and calculate the phonon dispersion of silicon in the unit cell using a \(2\times2\times2\) homogeneous **q**-grid:

```
cd workdir; cp -r ../inputs/unit_cell_222_grid/ .; cd unit_cell_222_grid/
mpirun -np 4 $QE/pw.x -nk 4 < si.scf.in > si.scf.out
mpirun -np 4 $QE/ph.x -nk 4 < si.ph.in > si.ph.out
mpirun -np 1 $QE/q2r.x < q2r.in > q2r.out
mpirun -np 1 $QE/matdyn.x < matdyn.in > matdyn.out
```

Input files are provided here:

```
-- si.scf.in
&control
calculation = 'scf'
restart_mode = 'from_scratch',
prefix = 'si',
pseudo_dir = './',
outdir='./'
/
&system
ibrav = 2,
celldm(1) = 10.20,
nat = 2,
ntyp = 1,
ecutwfc = 30.0
/
&electrons
diagonalization='david'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-7
/
ATOMIC_SPECIES
Si 28.086 Si.pz-vbc.UPF
K_POINTS automatic
6 6 6 0 0 0
ATOMIC_POSITIONS {crystal}
Si 0.00 0.00 0.00
Si 0.25 0.25 0.25
/
```

**si.ph.in**:

```
-- si.ph.in
Phonons of Silicon
&inputph
amass(1)= 28.0855,
prefix = 'si'
outdir = './'
ldisp = .true.
fildyn = 'si.dyn'
tr2_ph = 1.0d-12
nq1 = 2
nq2 = 2
nq3 = 2
/
```

**q2r.in**

```
-- q2r.in
&input
fildyn='si.dyn', flfrc = 'si.222.fc'
/
```

**matdyn.in**

```
-- matdyn.in
&input
asr='simple', amass(1)=28.0855,
flfrc='si.222.fc', fldyn='si.dyn.mat',
flfrq='si.freq', fleig='si.dyn.eig', q_in_cryst_coord = .false.
q_in_band_form = .true.
/
9
0 0 0 100
0.75 0.75 0.0 1
0.2500 1.0000 0.2500 100
0.000 1 0.000 100
0 0 0 100
0.5000 0.5000 0.5000 100
0 1 0 100
0.5000 1.0000 0 100
0.5000 0.5000 0.5000 100
/
```

After the calculation is completed use `plotband.x`

of Quantum Espresso to obtain the phonon dispersion data as an a xmgr file.

```
$QE/plotband.x
Input file > si.freq
Reading 6 bands at 702 k-points
Range: 0.0000 511.6990eV Emin, Emax, [firstk, lastk] > 0 600
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 0.0000
high-symmetry point: 0.7500 0.7500 0.0000 x coordinate 1.0607
...
high-symmetry point: 0.5000 0.5000 0.5000 x coordinate 5.3534
output file (gnuplot/xmgr) > si_ph_dispersion.xmgr
bands in gnuplot/xmgr format written to file si_ph_dispersion.xmgr
output file (ps) >
stopping ...
```

Now calculate the interatomic force constants of silicon using a \(2\times2\times2\) supercell and the \(\Gamma\)-point for the **q**-grid:

```
cd ../; cp -r ../inputs/222_super/ .; cd 222_super
mpirun -np 56 $QE/pw.x -nk 4 < si.scf_222.in > si.scf_222.out
mpirun -np 56 $QE/ph.x -nk 4 < si.ph_222.in > si.ph_222.out
mpirun -np 1 $QE/q2r.x < q2r.in > q2r.out
```

Input files are provided here:

```
-- si.scf_222.in
&control
calculation = 'scf'
restart_mode = 'from_scratch',
prefix = 'si',
pseudo_dir = './',
outdir='./'
/
&system
ibrav = 2,
celldm(1) = 20.40,
nat = 16,
ntyp = 1,
ecutwfc = 20.0
/
&electrons
diagonalization='david'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-7
/
ATOMIC_SPECIES
Si 28.086 Si.pz-vbc.UPF
K_POINTS automatic
3 3 3 0 0 0
ATOMIC_POSITIONS (crystal)
Si 0.00000000 0.00000000 0.00000000
Si -0.12500000 0.37500000 -0.12500000
Si 0.50000000 0.00000000 0.00000000
Si 0.37500000 0.37500000 -0.12500000
Si 0.00000000 0.50000000 0.00000000
Si -0.12500000 0.87500000 -0.12500000
Si 0.50000000 0.50000000 0.00000000
Si 0.37500000 0.87500000 -0.12500000
Si 0.00000000 0.00000000 0.50000000
Si -0.12500000 0.37500000 0.37500000
Si 0.50000000 0.00000000 0.50000000
Si 0.37500000 0.37500000 0.37500000
Si 0.00000000 0.50000000 0.50000000
Si -0.12500000 0.87500000 0.37500000
Si 0.50000000 0.50000000 0.50000000
Si 0.37500000 0.87500000 0.37500000
```

**si.ph.in**:

```
-- si.ph.in
Phonons of Silicon
&inputph
amass(1)= 28.0855,
prefix = 'si'
outdir = './'
ldisp = .true.
fildyn = 'si.dyn'
tr2_ph = 1.0d-12
nq1 = 1
nq2 = 1
nq3 = 1
/
```

**q2r.in**:

```
-- q2r.in
&input
fildyn='si.dyn', flfrc = 'si.222.fc'
/
```

Now perform phonon unfolding of the phonons of silicon calculated for a \(2\times2\times2\) supercell using `ZG.x`

:

```
mpirun -np 56 $QE/ZG.x -nk 56 < ZG_ph_unfold.in > ZG_ph_unfold.out
```

where the input file is:

```
-- ZG_ph_unfold.in
&input
flfrc='si.222.fc',
asr='simple', amass(1)=28.0855
q_in_cryst_coord = .false., q_in_band_form = .true.
ZG_conf = .false., ph_unfold = .true.
/
&phonon_unfold
dim1 = 2, dim2 = 2, dim3 = 2
ng1 = 14, ng2 = 14, ng3 = 14
flfrq='frequencies.dat', flweights='unfold_weights.dat'
/
9
0.000000 0.000000 0.000000 74
1.500000 1.500000 0.000000 1
0.500000 2.000000 0.500000 24
0.000000 2.000000 0.000000 70
0.000000 0.000000 0.000000 60
1.000000 1.000000 1.000000 60
0.000000 2.000000 0.000000 34
1.000000 2.000000 0.000000 50
1.000000 1.000000 1.000000 1
/
```

Here, lines 3-5 contain standard inputs for a phonon dispersion calculation (as for `matdyn.x`

). In line 4, we specify that we would like to use `ZG.x`

for phonon unfolding (`ph_unfold = .true.`

) and not to generate special displacements (`ZG_conf = .false.`

). In the namelist `&phonon_unfold`

we have: (i) **dimX** which are essentially the supercell dimensions, (ii) **ngX** which define the range of the reciprocal lattice vectors (and thus plane waves) that need to be employed for calculating the spectral weights. Check for convergence by increasing this number; already 14 is a tight value, (iii) **flfrq** is the name of the output file containing the frequencies for each band and **Q**-point and (iv) **flweights** is the name of the output file containing the spectral weights for each band and **Q**-point. At the end of the file we provide the high-symmetry path for which we would like to perform phonon unfolding. The input structure of path is in the same spirit with a `matdyn.x`

calculation. As for the electron band structure unfolding, the coordinates of the **Q**-points are obtained by multiplying the **q** wavevectors (defining the fundamental Brillouin zone) with the dimensions of the supercell, i.e. if **q** \(= [ x \,\,\, y \,\,\, z]\) then **Q** \(= [ {\rm dim1} \times x \,\, {\rm dim2} \times y \,\, {\rm dim3} \times z]\).

Now we have the phonon frequencies and their spectral weights, we can calculate the phonon spectral function using `pp_spctrlfn.x`

. Convert first the raw data in **frequencies.dat** and **unfold_weights.dat** into a different format using `plotband.x`

of QE (commands provided in `bands.sh`

):

```
$QE/plotband.x unfold_weights.dat < spectral_weights.in > spectral_weights.out
$QE/plotband.x frequencies.dat < energies.in > energies.out
sed -i '/^\s*$/d' spectral_weights.dat # remove empty lines
sed -i '/^\s*$/d' energies.dat
paste energies.dat spectral_weights.dat > tmp
awk '{print $1,$2,$4}' tmp > energies_weights.dat; rm tmp
```

You have just generated **energies_weights.dat** file in a similar format to a “.gnu” file where the first column has the momentum-axis values, the second the phonon frequencies, and the third the spectral weights. Now proceed with the calculation of the spectral function:

```
ibrun -n 28 $QE/pp_spctrlfn.x -nk 28 < pp_spctrlfn.in > pp_spctrlfn.out
cp spectral_function.dat ../../gnuplot/spectral_function_222.dat
cd ../../gnuplot/; gnuplot gp_pm3d_222.p; evince ph_unfolding_perfect_222.eps
```

where the input file should look like:

```
-- pp_spctrlfn.in
&input
flin = 'energies_weights.dat'
steps = 17952,
ksteps = 300,
esteps = 300,
kmin = 0,
kmax = 5.3534,
emin = -10
emax = 550
flspfn = 'spectral_function.dat'
/
```

Note

The description of the input flags is available in Exercise 2.

Figure 18 shows the silicon phonon spectral function (color map) which overlays with the phonon dispersion calculated in the unit cell (white lines), since the supercell is an exact periodic replica of the unit cell. One can also check the raw data in **unfold_weights.dat**, and observe that the weights of the phonon frequencies are either 1 or 0, making an one-to-one correspondence with the frequencies calculated in the unit cell.

Now go to `workdir`

, calculate the phonons of P-doped silicon, and perform phonon unfolding:

```
cd ../workdir/;cp -r ../inputs/222_super_P_doped .; cd 222_super_P_doped/
mpirun -np 56 $QE/pw.x -nk 4 < si.scf_222.in > si.scf_222.out
mpirun -np 56 $QE/ph.x -nk 4 < si.ph_222.in > si.ph_222.out
mpirun -np 1 $QE/q2r.x < q2r.in > q2r.out
mpirun -np 56 $QE/ZG.x -nk 56 < ZG_ph_unfold.in > ZG_ph_unfold.out
```

Note

Here we replace one Si atom with a P atom and provide the relaxed structure. Due to the excess electron the system becomes metalic.

The input files are provided:

```
-- si.scf_222.in
&control
calculation = 'scf'
restart_mode = 'from_scratch',
prefix = 'si',
pseudo_dir = './',
outdir='./'
/
&system
ibrav = 2,
celldm(1) = 20.40,
nat = 16,
ntyp = 2,
ecutwfc = 30.0
occupations='smearing'
degauss = 0.01
nosym =.true.
/
&electrons
diagonalization='david'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-7
/
ATOMIC_SPECIES
Si 28.086 Si.pz-vbc.UPF
P 30.974 P.pz-hgh.UPF
K_POINTS automatic
3 3 3 0 0 0
ATOMIC_POSITIONS (crystal)
P 0.0000000002 -0.0000000243 0.0000000412
Si -0.1261235427 0.3783706126 -0.1261235364
Si 0.4983582216 0.0016417759 0.0016417825
Si 0.3751699490 0.3751699540 -0.1255098692
Si 0.0016417706 0.4983582252 0.0016417800
Si -0.1261235463 0.8738764598 -0.1261235454
Si 0.4983582134 0.4983582133 0.0016417907
Si 0.3783706477 0.8738764517 -0.1261235465
Si 0.0016417985 0.0016418052 0.4983581945
Si -0.1255098925 0.3751699883 0.3751699378
Si 0.4983582056 0.0016417944 0.4983582063
Si 0.3751699659 0.3751699770 0.3751699562
Si 0.0016417902 0.4983582115 0.4983582106
Si -0.1261235587 0.8738764426 0.3783706794
Si 0.4999999969 0.5000000190 0.4999999746
Si 0.3751699807 0.8744900937 0.3751699437
```

**si.ph_222.in**:

```
-- si.ph_222.in
&inputph
amass(1)= 28.0855,
amass(2)= 30.974,
prefix = 'si'
outdir = './'
ldisp = .true.
fildyn = 'si.dyn'
tr2_ph = 1.0d-12
epsil = .false.
nq1 = 1
nq2 = 1
nq3 = 1
/
```

**q2r.in**:

```
-- q2r.in
&input
fildyn='si.dyn', flfrc = 'si.222.fc'
/
```

**ZG_ph_unfold.in**:

```
-- ZG_ph_unfold.in
&input
flfrc='si.222.fc',
asr='simple', amass(1)=28.0855, amass(2)= 30.974
q_in_cryst_coord = .false., q_in_band_form = .true.
ZG_conf = .false.
ph_unfold = .true.
/
&phonon_unfold
dim1 = 2, dim2 = 2, dim3 = 2
ng1 = 14, ng2 = 14, ng3 = 14
flfrq='frequencies.dat', flweights='unfold_weights.dat'
/
9
0.000000 0.000000 0.000000 74
1.500000 1.500000 0.000000 1
0.500000 2.000000 0.500000 24
0.000000 2.000000 0.000000 70
0.000000 0.000000 0.000000 60
1.000000 1.000000 1.000000 60
0.000000 2.000000 0.000000 34
1.000000 2.000000 0.000000 50
1.000000 1.000000 1.000000 1
```

Now we have the phonon frequencies and their spectral weights, we can calculate the phonon spectral function using `pp_spctrlfn.x`

as before. Convert first the raw data in **frequencies.dat** and **unfold_weights.dat** into a different format using `plotband.x`

of QE (commands provided in `bands.sh`

):

```
$QE/plotband.x unfold_weights.dat < spectral_weights.in > spectral_weights.out
$QE/plotband.x frequencies.dat < energies.in > energies.out
sed -i '/^\s*$/d' spectral_weights.dat # remove empty lines
sed -i '/^\s*$/d' energies.dat
paste energies.dat spectral_weights.dat > tmp
awk '{print $1,$2,$4}' tmp > energies_weights.dat; rm tmp
```

Proceed with the calculation and plot of the spectral function:

```
ibrun -n 28 $QE/pp_spctrlfn.x -nk 28 < pp_spctrlfn.in > pp_spctrlfn.out
cp spectral_function.dat ../../gnuplot/spectral_function_222_P_doped.dat
cd ../../gnuplot/; gnuplot gp_pm3d_222_P_doped.p; evince ph_unfolding_P_doped_222.eps
```

where the input file **pp_spctrlfn.in** is the same as before:

What are the changes in the phonon dispersion? Below we also provide the phonon spectral function calculated after replacing one Si atom with a P atom in a \(3\times3\times3\) supercell. Why are the changes less pronounced now?