# Phonons with DFPT
A small tutorial with cubic NaCl
## Download the material
Inputs and pseudo are available at this [link](https://sissa-my.sharepoint.com/:u:/g/personal/pdelugas_sissa_it/Ed00pvPgv7BFsxAIb7nSdlwBBeJz7bzdxK1g44rRI0sVdA?e=2miMV9)
* dowload with the browser the `NaCl_phons.tar.gz` file and expand it
```shell=
mkdir work-day2
cp NaCl_phons.tar.gz work-day2
cd work-dy2
tar xzf NaCl_phons.tar.gz
```
this should have create a directory `NaCl` with inside two sub-directories `inputs` and `pseudos`
* go inside `NaCl`
```shell=
cd NaCl
```
## important on KENET
* set the number to threads to 2 to avoid using hyperthreading
```shell=
export OMP_NUM_THREADS=2
```
* don't forget to activate conda qe environment:
```shell
conda activate /home/$USER/miniconda3/envs/qe
```
## scf calculation
The first input we use is the input for the scf calculation
```shell
cp inputs/nacl.scf.in ./
```
I have made few changes to the input of yesterday:
```
&system
ibrav=2
celldm(1)=10.58
nat=2
ntyp=2
ecutwfc = 60
ecutrho = 360
/
```
* use an optimized lattice constant **10.58**
* changed the wave function cutoff to **60** and set **ecutrho** to **360**
* this speeds up calcuation a bit
### lets run the calculation
```shell=
pw.x -i nacl.scf.in | tee out_scf
```
* then we run the phonon calculation for the $\Gamma$ point.
```shell=
ph.x -i phG.in | tee out_ph
```
* then we run the `dynmat` analysis to compute frequencies, dielectric tensor and other analysis
```shell
dynmat.x < input_dynmat | tee out_dyn
```
From this run we obtain the so call TO $\Gamma$ frequencies.
We obtain these values for the frequencies and the dielectric tensor:
```
# mode [cm-1] [THz] IR
1 -0.00 -0.0000 0.0000
2 -0.00 -0.0000 0.0000
3 -0.00 -0.0000 0.0000
4 166.89 5.0033 1.9591
5 166.89 5.0033 1.9591
6 166.89 5.0033 1.9591
Electronic dielectric permittivity tensor (F/m units)
2.589220 0.000000 0.000000
0.000000 2.589220 -0.000000
-0.000000 -0.000000 2.589220
... with zone-center polar mode contributions
6.008546 0.000000 0.000000
0.000000 6.008546 0.000000
0.000000 0.000000 6.008546
```
* now we run `dynmat.x` with a different input telling the program to add the non-analytic correction, simulating a long-wavelengh vector in the $x$ direction
```shell=
cp inputs/input_dynmat_LO ./
dynmat.x -i input_dynmat_LO | tee out_dynmat_LO
```
* in this case the frequencies of the polar modes' triplet are split with on mode mode that has a higher energy
```
IR activities are in (D/A)^2/amu units
# mode [cm-1] [THz] IR
1 -0.00 -0.0000 0.0000
2 -0.00 -0.0000 0.0000
3 0.00 0.0000 0.0000
4 166.89 5.0033 1.9591
5 166.89 5.0033 1.9591
6 254.24 7.6218 1.9591
```
* we should obtain a magnitude of the splitting of $90~\mathrm{cm}^{-1}$
## phonon dispersion
* we first copy the input for `ph.x` to compute the phonon frequencies in $4 \times 4 \times 4$ grid.
```shell=
cp ..inputs/phDisp.in ./
ph.x -i phDisp.in > out_phDisp &
```
* ... in the meantime that the calculation completes let's have a look at the input:
```
phonons on a q point mesh
&inputph
prefix='NaCl_B1'
tr2_ph = 1.0d-16
epsil = .true.
alpha_mix(10) = 0.3
outdir = 'temp_project'
fildyn = 'NaCl.dyn'
nk1 = 2, nk2 = 2, nk3 = 2
k1 = 0, k2 = 0, k3 = 0
ldisp = .true.
nq1 = 4
nq2 = 4
nq3 = 4
/
0.0 0.0 0.0
```
* `ldisp=.true.` enables the calculation on a uniform mesh of phonons'wave-vectors.
* `nq1,nq2, nq3` describes the mesh of phonons wave-vectors (q-points) that we want to compute. Only the points in the irreducible wedge will be actually computed, the other will be obtatined using the symmetries.
* `epsilon = .true.` enables the calculation of the response to macroscopic electric field.
* I also asked the program to use a reduced k-point mesh for the electonic calculations: `nk1,nk2,nk3`. This makes the calculation faster, we can finish it within the tutorial, but slightly inaccurate.
* `fildyn` value sets the names of the files where the program will save the dynamical matrices" one for each $\mathbf{q}$.
## Elaboration of the results with `q2r.x` and `mathdyn.x`
* we now use the results saved in the dynamical matrices files to compute the interatomic force constants. As we have used a $4 \times 4 \times 4$ $\mathbf{q}$-point mesh, we will be able to compute the interatomic force constants for all the ions contained in a $4\times4\times4$ supercell. This is done performing a inverse Fourier transform
$$\large
L_{\alpha,\{\tau,\mathbf{R}\}}^{\beta,\{\tau^\prime,\mathbf{R}^\prime \}} = \underset{\mathbf{q}}\sum{e^{i{(\mathbf{R-R^\prime})\cdot\mathbf{q}}}
D_{\alpha,\tau}^{\beta,\tau^\prime}(\mathbf{q})}=L_{\alpha,\tau}^{\beta,\tau^\prime}(\mathbf{R-R^\prime})
$$
1. because of the translational invariance the interatomic force constant depends only on the reciprocal lattice vector separating the unit cells containing the 2 ions.
2. The resolution of the interatomic force constants is restricted to separation vectors within the supercell. For longer vectors ( they can alway be seen as the sum of a superlattice vector plus a residual vector contained in the supercell) we just obtain the interatomic force constant corresponding the residual vector.
* To see this consider that $\mathbf{q}$-kpoints in the mesh are all of the form $\sum_i\frac{n_i}{4}\mathbf{B_i}$ where the $\mathbf{B_i}$s are the generators of the reciprocal lattice. By definition of reciprocal lattice we have that $$\large \frac{n_i}{4}\mathbf{B_i}\cdot 4\mathbf{A_j} = 2\pi \,\delta_{ij}\,n_i$$
* For polar systems, as our NaCl, the interaction contains also the long range electrostatic term, which `q2r.x` subtracts from the IFCs to shorten the ion-ion interaction range and improve the convergence of the Fourier expantion without using too much dense $\mathbf{q}$-points meshes.
* we copy the input inside the working directory and run the `q2r.x`
```shell=
cp inputs/input_q2r ./
q2r.x -i input_q2r | tee out_q2r
```
* the input of `q2r.x` is very simple we just need to tell the program the file names for the dynamical matrices `fildyn`, the name of the file where to save the interatomic force constants `flfrc`, and how to deal with the inaccuracies of the Born Charges, whether to impose the acustic sum rule and how `zasr`.
```
&input
fildyn = 'NaCl.dyn'
flfrc = 'NaCl.frc'
zasr = 'crystal'
/
```
* we will now use the interatomic force constants to compute the dynamical matrices at arbitrary $\mathbf{q}$-points with `matdyn.x`. This program reverses the operation done in `q2r.x`. $$\large D_{\alpha,\tau}^{\beta,\tau^\prime}(\mathbf{q}) = \underset{R}\sum{e^{-i\mathbf{q}\cdot\mathbf{R}}L_{\alpha,\tau}^{\beta,\tau^\prime}(\mathbf{R})}$$ now $\mathbf{q}$ is arbitrary, and even if the wavector $\mathbf{q}$ is not compatible with the supercell lenght, we are just assuming that outside the supercell the interatomic force constants vanish. In the case of polar systems this assumption makes sense only because we have removed the long-range electrostatic interaction, that now `matdyn.x` adds back to the dynamical matrix.
* copy the first input for `matdyn.x` and run it.
```shell=
cp inputs/input_matdyn_dos ./
matdyn.x -i input_matdyn_dos | tee out_matdyn_dos
```
* we ask matdyn to compute the vibration density of states so we need to specify a dense enough uniform grid of wave-vectors, in this case $8\times 8 \times 8$:
```
&input
flfrc='NaCl.frc'
asr ='crystal'
nk1 = 8
nk2 = 8
nk3 = 8
dos=.true.
fldos='NaCl.dos'
/
```
* The VDOS is saved in the file `NaCl.dos` that contains the total density of states and the projections on each atom species. Can be easily plotted for example with xmgrace:
```shell=
xmgrace -nxy NaCl.dos
```
Sign in with Wallet
Connect another wallet