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.

Convergence of self-consistent phonon (SCPH) algorithm for 2500 K. The individual harmonic parameters shown on the left and the norm of the change in parameters on the right.

Phonon dispersions for BCC Ta from a series of SCPH models.

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

Source code

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


# read FCP
fcp = ForceConstantPotential.read('fcps/fcp_sixth_order.fcp')
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))

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


# read fcps
fcps = dict()
for T in temperatures:
    fcps[T] = ForceConstantPotential.read('fcps/scph_T{}.fcp'.format(T))

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

# read parameter trajs
parameter_trajs = np.loadtxt('scph_trajs/scph_parameters_T2000')
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))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

# 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()
ax1 = fig.add_subplot(111)

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')