# 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. BCC tantalum.

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 500, 1500 and 2500 K. The phonon dispersion is constructed for each temperature and a clear 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¶

Run self-consistent phonons
examples/advanced_topics/temperature_dependent_phonons/2_run_scph.py

import os
import numpy as np
from hiphive import ClusterSpace, ForceConstantPotential
from hiphive.self_consistent_phonons import self_consistent_harmonic_model
from hiphive.calculators import ForceConstantCalculator

# parameters
cutoffs = [6.0]
temperatures = [2000, 1000, 300]
size = 6

# scp parameters
n_structures = 30
n_iterations = 50
alpha = 0.2

prim = fcp.primitive_structure

# setup scph
cs = ClusterSpace(prim, cutoffs)
supercell = prim.repeat(size)
fcs = fcp.get_force_constants(supercell)
calc = ForceConstantCalculator(fcs)

# run scph
os.makedirs('scph_trajs/', exist_ok=True)
for T in temperatures:
parameters_traj = self_consistent_harmonic_model(
supercell, calc, cs, T, alpha, n_iterations, n_structures)
fcp_scph = ForceConstantPotential(cs, parameters_traj[-1])
fcp_scph.write('fcps/scph_T{}.fcp'.format(T))
np.savetxt('scph_trajs/scph_parameters_T{}'.format(T), np.array(parameters_traj))


Plot SCPH phonon dispersions
examples/advanced_topics/temperature_dependent_phonons/4a_plot_scph_phonons.py

import numpy as np
import matplotlib.pyplot as plt
from hiphive import ForceConstantPotential
from helpers import get_dispersion

# parameters
temperatures = [2000, 1000, 300]
dim = 6

fcps = dict()
for T in temperatures:

# collect dispersion
dispersions = dict()
for T, fcp in fcps.items():
dispersions[T] = get_dispersion(fcp, dim)

# read parameter trajs
delta_parameters = [np.linalg.norm(p-p2) for p, p2 in zip(parameter_trajs, parameter_trajs[1:])]

# setup plot
fig = plt.figure(figsize=(8, 3.5))

# plot scph convergence
ax1.plot(parameter_trajs)
ax1.set_xlabel('SCPH iteration')
ax1.set_ylabel('Parameters')
ax1.set_xlim([0, len(parameter_trajs)])

ax2.plot(delta_parameters)
ax2.set_xlabel('SCPH iteration')
ax2.set_ylabel('$\\Delta$ Parameters')
ax2.set_xlim([0, len(parameter_trajs)])
ax2.set_ylim(bottom=0.0)

fig.tight_layout()
fig.savefig('scph_parameter_convergence.svg')

# setup plot dispersions
fig = plt.figure()

cmap = plt.get_cmap('viridis')
colors = {T: cmap(i/(len(dispersions)-0.9)) for i, T in enumerate(temperatures)}

# plotting dispersion
for T, (qnorms, freqs) in dispersions.items():
color = colors[T]
for q, freq, in zip(qnorms, freqs):
ax1.plot(q, freq, color=color)
ax1.plot(np.nan, np.nan, color=color, label='{}K'.format(T))

ax1.set_xlabel('Wave vector $\\vec{q}$')
ax1.set_ylabel('Frequency (meV)')
ax1.legend()

# set qpts as x-ticks
qpts = [0] + [q[-1] for q in qnorms]
qpts_names = ['N', r'$\Gamma$', 'H', r'$\Gamma$']
for q in qpts:
ax1.axvline(q, color='k', lw=0.6)
ax1.axhline(y=0.0, color='k', ls='-', lw=1.0)
ax1.set_xticks(qpts)
ax1.set_xticklabels(qpts_names)
ax1.set_xlim([qpts[0], qpts[-1]])
ax1.set_ylim([0, 20])

fig.tight_layout()
fig.savefig('scph_phonon_dispersions.svg')