# Phonon calculation (finite displacement) {#sec:phonon}

RESCU can perform phonon computations using the finite-displacement
method [\[PLK\]](#PLK){.citation}. The workflow of a phonon computation
goes as follows. Firstly, RESCU extends the structure to a designated
supercell structure. Then it determines its symmetries, generates
non-equivalent displacements and carries out a self-consistent
calculation for every atomic configuration. Input/output files are
generated for each non-equivalent displacement such that the
self-consistent calculations can be executed in parallel. Upon
completion, the output files contain the Hellmann-Feynman forces
experienced by each atom. Finally, RESCU calculates the force constants,
dynamical matrices and related quantities such as the phonon band
structure. In this section, we present the example of a bulk fcc
silicon.

## Example: silicon {#sec:phonon_si}

Let's create the file `si_real_phonon.input` as follows.

``` matlab
info.savepath = './results/si_real_phonon'
info.calculationType = 'phonon'
atom.element = ones(2,1)
atom.xyz = 10.25*[0 0 0;0.25 0.25 0.25]
domain.latvec=10.25/2*[0 1 1;1 0 1;1 1 0]
domain.lowres = 0.4
element(1).species = 'Si'
element(1).path = './Si_TM_LDA.mat'
functional.list = {'XC_LDA_X','XC_LDA_C_PW'}
kpoint.gridn = [3,3,3]
ph.mode = 'phonon'
ph.dr = 0.2
ph.supercellDimension = [2 2 2]
ph.ASR = 1
ph.qpointSympoints = {{[0,0,0] 0.5*[1,0,1]} {0.5*[2,1,1] [0,0,0]} {[0,0,0] 0.5*[1,1,1]}}         
ph.labels = {'G','X','X1','G','G','L'}
```

<figure>
<img src="data/si-phonon.png" id="fig:si_ph_bs"
alt="data/si-phonon.png" />
<figcaption>Phonon band structure of fcc silicon</figcaption>
</figure>

The meaning of the phonon parameters (`ph.xxx`) is described below.

-   `info.calculationType` is set to phonon so that RESCU performs a
    phonon calculation
-   `ph.mode` is a string that determines which step of a phonon
    calculation is to be done. You will most often use the
    `displacement`, `phonon` and `post` modes. There are eight options:
    -   `bands`: RESCU reads the FORCE_CONSTANTS file, constructs the
        effective dynamical matrix and computes the phonon frequencies
        along the designated k-point line.
    -   `displacement`: RESCU determines the non-equivalent
        displacements and prints out input and XYZ files for each of
        them.
    -   `gamma`: is of boolean type. When it is true, the supercell size
        is automatically
        1$\times$`<!-- -->`{=html}1$\times$`<!-- -->`{=html}1. This is
        used to obtain exact phonon frequencies at the $\Gamma$-point
        only.
    -   `forcesets`: RESCU combines the atomic forces stored in the
        files FORCESN ($N=1,...,N_{disp}$) into a single file:
        FORCE_SETS.
    -   `forceconstants`: RESCU uses the FORCE_SETS file together with
        symmetry operations to calculate the force constants of the
        supercell structure, and then stores the results in a file
        FORCE_CONSTANTS.
    -   `phonon`: RESCU performs all calculations necessary to obtain a
        phonon band structure:
        `symmetry`,`displacement`,`self-consistent`,`forcesets`,`forceconstants`
        and `bands`.
    -   `post`: RESCU uses the self-consistent
        results`forces to perform the following calculations:`forcesets`,`forceconstants`, and`bands\`\`.
    -   `symmetry`: RESCU analyzes the structure and gives information
        about the number of symmetrical operations.
-   `ph.dr` determines the magnitude of the atomic displacements in
    Bohr. A typical value is around 0.1 Bohr. It is advisable to take dr
    to be the same order as domain.lowres to mitigate the eggbox effect.
-   `ph.supercellDimension` is a most important parameter which
    determines the size of the supercells used to calculate the
    interatomic force constants. Large supercells generally yield better
    results, but demand more computational resources.
-   `ph.ASR` switch for the Acoustic-Sum-Rule (ASR) correction scheme.
    If ph.ASR is true, the ASR is enforced. Theoretically, the
    $\Gamma$-point of a 3D bulk system should have 3 zero frequencies,
    which are part of the so-called acoustic phonon branches. The ASR
    correction scheme usually improves the accuracy at the
    $\Gamma$-point.
-   `ph.qpointSympoints` is a cell array of high symmetry point symbols
    or row-vectors in fractional coordinates that gives the path of
    k-points when calculating phonon band structures. One may write, for
    instance, ph.qpointSympoints={$k_0$ $k_1$ \... $k_n$}, where tuples
    are typically are high-symmetry points of the first Brillouin zone.
    Here we define three such lines which are all contained in a cell
    array.
-   `ph.labels` is a cell array of strings labeling the k-points {$k_0$
    $k_1$ \... $k_n$} given in `ph.qpointSympoints`. The labels will
    show on the phonon band structure plot.
-   `ph.qpointGridn` (not used in example) is an integer that determines
    the number of k-points along the line given by `ph.qpointSympoints`.
-   `ph.plot` (not used in example) is of boolean type and it tells
    RESCU whether to plot the phonon band structure and save it as a
    .png file.

You may perform the calculation typing

``` shell
rescu -i si_real_phonon.input
```

Upon completion, the following files have been produced: disp.dat,
FORCE_SETS and FORCE_CONSTANTS. Their units are Å, eV/Å and eV/Å$^2$,
respectively. The phonon band structure of silicon is shown in
`fig:si_ph_bs`{.interpreted-text role="numref"}.

## Example: black phosphorus

<figure>
<img src="data/bp_ph_bs.png" id="fig:bp_ph_bs"
alt="data/bp_ph_bs.png" />
<figcaption>Phonon band structure of a black phosphorus
monolayer.</figcaption>
</figure>

In section [1.1](#sec:phonon_si), we show how the `phonon` mode of
phonon calculations computes the phonon band structure of a material in
one shot. It is often impractical or inefficient to use the `phonon`
mode for research activities. In this section, we thus show how to
perform a phonon band structure calculation step-by-step using black
phosphorus as an example. First, create the following input file under
the name `bp_real_disp.input`.

``` shell
info.calculationType = 'phonon'
info.savepath = 'results/bp_real_disp'
atom.element = ones(4,1)
atom.xyz = [0.0 0.053831313 4.704675810
1.660810589 1.549175746 4.704675810
1.660810589 2.343841383 2.572460227
0.0 3.839185675 2.572460227]
domain.latvec = [3.3216214769 0.0 0.0
0.0 4.5800202998 0.0
0.0 0.0 10.0]
domain.lowres = 0.4
units.length = 'a'
element(1).species = 'P'
element(1).path = './P_TM_LDA.mat'
ph.supercellDimension = [4 4 1]
ph.dr = 0.4
ph.mode = 'displacement'
```

Then execute it as follows

``` shell
rescu -i bp_real_disp.input
```

The structure of black phosphorus is such that there are four
non-equivalent displacements. RESCU thus creates four input files,
`1displaced.input`, `2displaced.input`, `3displaced.input` and
`4displaced.input` and four corresponding XYZ files containing the
atomic coordinates and magnetic moments. For example, `1displaced.input`
looks as follows:

``` shell
domain.lowres = 0.4
units.length = 'a'
element(1).species = 'P'
element(1).path = './P_TM_LDA.mat'
kpoint.gridn = [1,1,1]
mixing.tol = [1e-6,1e-6]
info.calculationType = 'self-consistent'
info.savepath = 'results/1displaced'
domain.latvec = [+1.32864859076e+01 0.0 0.0
0.0 +1.83200811992e+01 0.0
0.0 0.0 10.0]
atom.xyz = '1displaced.xyz'
atom.magmom = '1displaced.xyz'
atom.nvalence = 320
units.latvec = 'a'
units.xyz = 'a'
```

RESCU reproduces the phonon calculation input file at the exception of a
few parameters. Firstly, `info.calculationType` is changed to
`self-consistent` to compute the atomic forces and `info.savepath` is
set to a displacement-dependent name to avoid overwriting other results.
Next, the lattice vectors `domain.latvec`, atomic coordinates
`atom.xyz`, magnetic moments `atom.magmom` and total valence charge
`atom.nvalence` are changed to their supercell counterparts. Finally,
`units.latvec` and `units.xyz` are set to angstroms since the phonon
calculator interface employs those units.

One may then modify the self-consistent input files as he sees fit and
proceed with the self-consistent calculations.

``` shell
rescu -i 1displaced.input
rescu -i 2displaced.input
rescu -i 3displaced.input
rescu -i 4displaced.input
```

Note that these calculations can be submitted independently on a
cluster. After completing all calculations, the forces are saved in
`results/xdisplaced.mat`.

Finally, we obtain the band structure of black phosphorus by performing
a phonon calculation in `post` mode (which stands for post-calculation).
As mentioned in section [1.1](#sec:phonon_si), the `post` mode reads the
supercell forces, computes the force sets, force constants, dynamical
matrix and band structure. In the present exercice, this is achieved
creating the file `bp_real_post.input` as follows.

``` shell
info.calculationType = 'phonon'
info.savepath = 'results/bp_real_post'
atom.element = ones(4,1)
atom.xyz = [0.0 0.053831313 4.704675810
1.660810589 1.549175746 4.704675810
1.660810589 2.343841383 2.572460227
0.0 3.839185675 2.572460227]
domain.latvec = [3.3216214769 0.0 0.0
0.0 4.5800202998 0.0
0.0 0.0 10.0]
domain.lowres = 0.4
units.length = 'a'
element(1).species = 'P'
element(1).path = './P_TM_LDA.mat'
ph.supercellDimension = [4 4 1]
ph.dr = 0.4
ph.mode = 'post'
ph.qpointSympoints = {[0 0 0] [0 1/2 0] [1/2 1/2 0] [1/2 0 0] [0 0 0]}
ph.labels = {'G','Y','M','X','G'}
ph.plot = true
```

The file looks just like `bp_real_disp.input`, but we have defined the
k-point line via the keywords `ph.qpointSympoints` and `ph.labels`. We
also set `ph.plot` to `true` to force RESCU to display the band
structure after computing it. The band structure is shown in
`fig:bp_ph_bs`{.interpreted-text role="numref"}.

::: {#citations}

[PLK]{#PLK .citation-label}

:   Parlinski, K., Z. Q. Li, and Y. Kawazoe.
    `First-principles determination of the soft mode in cubic ZrO 2. <10.1103/PhysRevLett.78.4063>`{.interpreted-text
    role="doi"} Physical Review Letters 78.21 (1997): 4063.
:::
