# Self-consistent phonons¶

In this section the usage of the basic self-consistent phonons module in
**hiphive** is demonstrated. We use the same example as in the
effective harmonic models tutorial
i.e. FCC nickel with an EMT calculator as reference potential.

The basic algorithm for self-consistent phonons in **hiphive** resembles
the algorithm in [ErrCalMau14].

First an harmonic model \(\Phi(\mathbf{x}_0)\) is generated from the initial parameter vector \(\mathbf{x}_0\) obtained from fitting with structures with small random displacements.

This initial harmonic model, \(\Phi(\mathbf{x}_0)\), can be used to generate phonon rattled structures at the desired temperature. These structures are used to fit new parameters \(\mathbf{x}_\text{new}\). In order for the algorithm to converge smoothly we compute \(\mathbf{x}_1\) as

Next, structures are generated using \(\Phi(\mathbf{x}_1)\) and the procedure is repeated. The \(n\)-th step in the algorithm would formally be

where \(H\) denotes the reference potential.

In this example, self-consistent phonons are obtained for 200, 500, 800, 1100 and 1400 K. The phonon dispersion is constructed for each temperature and a clear yet small temperature dependency can be observed. This can be compared to the effective harmonic models.

Please note that in practice you usually would use a larger number of structures leading to smoother convergence.

## Source code¶

`tutorial/advanced/self_consistent_phonons/run_self_consistent_phonons.py`

```
from hiphive.structure_generation import generate_rattled_structures
from ase.build import fcc110
from ase.io import write
from ase.calculators.emt import EMT
from ase.optimize import BFGS
# parameters
number_of_structures = 5
rattle_std = 0.05
surface_size = (3, 4, 8)
structures_fname = 'rattled_structures.extxyz'
# setup atoms and calculator
atoms_ideal = fcc110('Cu', size=surface_size)
calc = EMT()
atoms_ideal.center(vacuum=20, axis=2)
atoms_ideal.pbc = True
# relax structure
atoms_ideal.set_calculator(calc)
dyn = BFGS(atoms_ideal)
converged = dyn.run(fmax=0.0001, steps=1000)
# generate rattled structures
structures = generate_rattled_structures(
atoms_ideal, number_of_structures, rattle_std)
for structure in structures:
structure.set_calculator(calc)
forces = structure.get_forces()
displacements = structure.positions - atoms_ideal.get_positions()
structure.new_array('displacements', displacements)
structure.new_array('forces', forces)
structure.positions = atoms_ideal.get_positions()
structure.calc = None
# save structures
write(structures_fname, structures)
```