Examples of calculations with ZG.x ================================== Author: Marios Zacharias .. highlight:: none .. 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. .. _backtotop: :ref:`Exercise 1 `: Generation of the ZG configuration and DFT-ZG calculation. :ref:`Exercise 2 `: Temperature-dependent band structures. :ref:`Exercise 3 `: Anharmonic lattice dynamics. :ref:`Exercise 4 `: Temperature-dependent optical spectra. :ref:`Exercise 5 `: Multiphonon diffuse scattering intensity. :ref:`Exercise 6 `: Phonon unfolding. :doc:`Input flags ` A complete list of tutorials and files are available in this `link `_. .. _ref1: 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 :math:`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 :ref:`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 :math:`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 .. figure:: /doc/images/TutorialZG/ZG_si_ph_dispersion.png :width: 750px **Figure 1**: Phonon dispersion relation of silicon. 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 :math:`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-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 (:ref:`Exercise 3 `). For the next exercises we will only need the charge-density files from **equil-si.save** and **ZG-si.save**. Thus remove: :: rm -r *wfc* si.save _ph0/ equil-si.save/*wfc* ZG-si.save/*wfc* :ref:`Back to top ` .. _ref2: 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:`3\times3\times3` 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 .. figure:: /doc/images/TutorialZG/ZG_si_band_structure.png :width: 350px **Figure 2**: Band structure of silicon along Γ-X using 50 **k**-points. 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-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** *=* [ :math:`m \times x n \times y p \times z` ], where integers :math:`m, n, p` define an :math:`m \times n \times p` supercell. If you did not complete :ref:`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-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 :math:`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 .. figure:: /doc/images/TutorialZG/ZG_TE_spctrl_fn.png :width: 750px **Figure 3**: Left: Thermal ellipsoids of silicon obtained by folding the nuclei of a :math:`50\times50\times50` ZG configuration back to the unit-cell. Arrows represent the ZG displacements. Right: Convergence of the zero-point renormalization of silicon as a function of supercell size. Green and black discs represent data obtained using ZG configurations generated with and without modes in set *A*, respectively. In the inset we show the corresponding spectral functions calculated using the band structure unfolding technique. Now prepare a band structure unfolding calculation using the :math:`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-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-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-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 :math:`8\times8\times8` ZG supercell, more **K**-points, and empty bands. .. figure:: /doc/images/TutorialZG/ZG_spectral_fn_BSU.png :width: 750px **Figure 4**: Left: Spectral function of silicon along Γ-Χ using a :math:`3\times3\times3` ZG configuration generated for T = 0 K. Right: Converged result as reported in Phys. Rev. Research 2, 013357, (2020). :ref:`Back to top ` .. _ref3: Exercise 3a ----------- In this exercise we will learn how to calculate temperature-dependent anharmonic phonons of Zr at *T=1188* K using :math:`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) <10.1080/14786440408520575>`_. 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: .. math:: :nowrap: \begin{eqnarray} \tilde{C}_{p \kappa \alpha, p' \kappa' \alpha'}(T) = \bigg< \frac{\partial U} {\partial \tau_{p \kappa \alpha} \partial \tau_{p' \kappa' \alpha'}} \bigg>_T. \end{eqnarray} 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 :math:`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 .. figure:: /doc/images/TutorialZG/ZG_Zr_harm_ph_dispersion.png :width: 650px **Figure 5**: Harmonic phonon dispersion relation of bcc Zr. Is the matrix of IFCs positive definite (i.e. all phonon frequencies are positive)? Unlike Si phonon dispersion in :ref:`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 :ref:`Exercise 1 `. We remark that here we are using a :math:`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 :ref:`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 :math:`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 :math:`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 :ref:`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 .. figure:: /doc/images/TutorialZG/Zr_poly_ph_dispersion.png :width: 650px **Figure 6**: Comparison of the harmonic phonon dispersion of polymorphous and monomorphous Zr. 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 .. figure:: /doc/images/TutorialZG/Zr_01_ph_dispersion.png :width: 650px **Figure 7**: Comparison of the phonon dispersion of Zr at *T=1188* K obtained at the first iteration (blue) with the phonon dispersion of the polymorphous network (red) and the harmonic phonon dispersion (green). 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 .. figure:: /doc/images/TutorialZG/Zr_02_ph_dispersion.png :width: 650px **Figure 8**: Comparison of the phonon dispersion of Zr at *T=1188* K obtained at the second iteration (brown) with the one at the first iteration (blue), the phonon dispersion of the polymorphous network (red) and the harmonic phonon dispersion (green). 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: .. figure:: /doc/images/TutorialZG/Zr_03_04_ph_dispersion.png :width: 650px **Figure 9**: Converged anharmonic phonon dispersion of Zr at *T=1188* K as calculated with A-SDM for a :math:`2\times2\times2` supercell. 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 :math:`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 :math:`2\times2\times2` supercells (i.e. using **1188.00K\_iter\_03.fc**). The result is provided in the plot below. .. figure:: /doc/images/TutorialZG/ZG_Zr_333_ph_dispersion.png :width: 650px **Figure 10**: Converged anharmonic phonon dispersion of Zr at *T=1188* K as calculated with A-SDM for a :math:`3\times3\times3` supercell. 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 (:math:`||.||`) of the leading component of the IFCs matrix (:math:`\tilde{C}_{1 1 \alpha, 1 1 \alpha} (T)`). This information is found from the first :math:`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 :math:`||\tilde{C}_{1 1 \alpha, 1 1 \alpha}||` and the difference is taken with respect to :math:`|| \tilde{C}_{1 1 \alpha, 1 1 \alpha} (T) ||` of the last iteration. Note also that IFCs unit has been converted from Ry/Bohr :math:`^2` to eV/ Å :math:`^2`. .. figure:: /doc/images/TutorialZG/ZG_Zr_norm.png :width: 650px **Figure 11**: Frobenius norm of the leading component of the IFCs matrix calculated with A-SDM for a :math:`2\times2\times2` supercell as a function of iterations. + 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 :math:`\big< \partial U / \partial \tau_{\kappa \alpha} \big>_T = 0`. This requires to solve self-consistently the equation: .. math:: :nowrap: \begin{eqnarray} \tau^{\rm new}_{p \kappa \alpha} (T) = \tau_{p \kappa \alpha} (T) + \tilde{C}^{-1}_{p \kappa \alpha, p \kappa \alpha} (T) \big< F_{p \kappa \alpha} \big>_T, \end{eqnarray} using at each iterative step the calculated :math:`\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 :math:`m\bar{3}m`), then the code will give :math:`\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 :ref:`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. :ref:`Back to top ` .. _ref3b: Exercise 3b ----------- In this exercise we will learn how to calculate temperature-dependent anharmonic phonons of PbTe at *T=300* K using :math:`2\times2\times2` ZG supercells and the A-SDM. Unlike Zr, seen in :ref:`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. Chem. Theory Comput. 15, 5845 (2019) `_, 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: :: &control scf.in 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**. .. figure:: /doc/images/TutorialZG/ZG_PbTe_harm_ph_dispersion.png :width: 650px **Figure 12**: Harmonic phonon dispersion of PbTe. 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 :math:`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 = 700 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**. .. figure:: /doc/images/TutorialZG/ZG_PbTe_01_ph_dispersion.png :width: 650px **Figure 13**: Comparison of the phonon dispersion of PbTe at *T=300* K obtained at the first iteration (red) with the harmonic phonon dispersion (green). 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**: .. figure:: /doc/images/TutorialZG/ZG_PbTe_iters_ph_dispersion.png :width: 650px **Figure 14**: A-SDM phonon dispersions of PbTe at *T=300* K obtained at the first (red), second (blue), and third (orange) iterations. The harmonic phonon dispersion (green) is shown for comparison. 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**: .. figure:: /doc/images/TutorialZG/ZG_PbTe_T_ph_dispersion.png :width: 650px **Figure 15**: A-SDM phonon dispersions of PbTe at *T=300* K (red), *T=500* K (blue), and *T=700* K (orange) obtained for the first iteration. The harmonic phonon dispersion (green) is shown for comparison. :ref:`Back to top ` .. _ref4: Exercise 4 ---------- In this exercise we will learn how to calculate the phonon-assisted optical spectra using the example of silicon and the :math:`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 :ref:`Exercise 1 ` and :ref:`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 :math:`3\times3\times3` equilibrium supercell and 30 random **K**-points. The aim is to calculate the imaginary part of the dielectric function [:math:`\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 :math:`\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 :math:`\varepsilon_2(\omega)`, (ii) **epsr_equil-si.dat** containing the real part of the dielectric function :math:`\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 :math:`3\times3\times3` ZG supercell and 30 random **K**-points. We will essentially repeat all steps performed for calculating the optical spectra of the :math:`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-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 :ref:`Exercise 1 ` and :ref:`Exercise 2 ` you can copy **ZG-si.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-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-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-si'``. Now run the following to calculate :math:`\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-si.dat** containing :math:`\varepsilon_2(\omega)`, (ii) **epsr_ZG-si.dat** containing the real part of the dielectric function :math:`\varepsilon_1(\omega)`, and (iii) **eels_ZG-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-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 :math:`\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. .. figure:: /doc/images/TutorialZG/ZG_epsilon2.png :width: 750px **Figure 16**: **a** Imaginary part of the dielectric function of silicon for photon energies up to 4.5 eV. Red and blue lines represent the spectra calculated with the nuclei at their equilibrium positions (direct transitions only) and at their ZG coordinates (indirect + direct transitions) for 0 K. **b** Same as in **a** but in logarithmic scale. 30 random **K**-points are used to sample the Brillouin zone. **c** Variation of the imaginary part of the dielectric function of silicon in logarithmic scale with the number of random **K**-points used to sample the Brillouin zone. **d** Same as in **b** but for 200 random **K**-points. :ref:`Back to top ` .. _ref5: 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 :math:`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 :ref:`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 :math:`100\times100` and :math:`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. .. figure:: /doc/images/TutorialZG/ZG_strf.png :width: 750px **Figure 17**: All-phonon diffuse scattering intensity of graphene at *T* = 300 K calculated using the ZG configuration (upper part) and the exact LBJ expression (lower part) for different **q**-grid samplings. :ref:`Back to top ` .. _ref6: 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 :math:`2\times2\times2` silicon supercell; then we will show how to calculate the effect on the phonons when a :math:`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 :math:`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 :math:`2\times2\times2` supercell and the :math:`\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 :math:`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** :math:`= [ x \,\,\, y \,\,\, z]` then **Q** :math:`= [ {\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 :ref:`Exercise 2 `. .. figure:: /doc/images/TutorialZG/ZG_ph_unfolding_perfect_222.png :width: 750px **Figure 18**: Phonon spectral function of silicon obtained by unfolding the phonons calculated with the atoms clamped at their high-symmetry position in a :math:`2\times2\times2` supercell. 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: .. figure:: /doc/images/TutorialZG/ZG_ph_unfolding_P_doped_222.png :width: 750px **Figure 19**: Phonon spectral function of P-doped silicon calculated in a :math:`2\times2\times2` supercell. 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 :math:`3\times3\times3` supercell. Why are the changes less pronounced now? .. figure:: /doc/images/TutorialZG/ZG_ph_unfolding_P_doped_333.png :width: 750px **Figure 20**: Phonon spectral function of P-doped silicon calculated in a :math:`3\times3\times3` supercell. :ref:`Back to top `