# 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

$\mathbf{x}_1 = \alpha\mathbf{x}_\text{new} + (1-\alpha)\mathbf{x}_0$

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

$\mathbf{x}_{n+1} = \alpha \, \underset{\mathbf{x}}{\text{argmin}} (\left < (H - \Phi(\mathbf{x}))^2 \right >_{\Phi(\mathbf{x}_{n})}) + (1-\alpha)\mathbf{x}_{n},$

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¶

Prepare data
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)