# 15. Phonon calculation (finite displacement)

RESCU can perform phonon computations using the finite-displacement method [PLK]. 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.

## 15.1. Example: silicon

Let’s create the file `si_real_phonon.input`

as follows.

```
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'}
```

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\)1\(\times\)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

```
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. 15.1.1.

## 15.2. Example: black phosphorus

In section 1.1, 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`

.

```
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

```
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:

```
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.

```
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, 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.

```
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. 15.2.1.

Parlinski, K., Z. Q. Li, and Y. Kawazoe. First-principles determination of the soft mode in cubic ZrO 2. Physical Review Letters 78.21 (1997): 4063.