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