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.

../_images/self_consistent_phonons.svg

Convergence of the harmonic parameters as a function of iterations (left) in the self-consistent phonon algorithm for 1400K. Phonon dispersions for different temperatures (right).

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

Source code

Run self-consistent phonons
tutorial/advanced/self_consistent_phonons/run_self_consistent_phonons.py
import numpy as np
import matplotlib.pyplot as plt

from ase import Atoms
from ase.build import bulk
from ase.calculators.emt import EMT

from hiphive import ClusterSpace, ForceConstantPotential
from hiphive.self_consistent_phonons import self_consistent_harmonic_model
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms


THz2meV = 4.13567


def get_band(q_start, q_stop, N):
    """ Return path between q_start and q_stop """
    return np.array([q_start + (q_stop-q_start)*i/(N-1) for i in range(N)])


def get_dispersion(prim, fcp):
    # setup phonopy and get FCs
    atoms_phonopy = PhonopyAtoms(symbols=prim.get_chemical_symbols(),
                                 scaled_positions=prim.get_scaled_positions(),
                                 cell=prim.cell)
    phonopy = Phonopy(atoms_phonopy, supercell_matrix=size*np.eye(3),
                      primitive_matrix=None)
    supercell = phonopy.get_supercell()
    supercell = Atoms(cell=supercell.cell, numbers=supercell.numbers, pbc=True,
                      scaled_positions=supercell.get_scaled_positions())

    fcs = fcp.get_force_constants(supercell)
    phonopy.set_force_constants(fcs.get_fc_array(order=2))

    # get phonon dispersion
    phonopy.set_band_structure(bands)
    qvecs, qnorms, freqs, _ = phonopy.get_band_structure()
    freqs = [f * THz2meV for f in freqs]
    return qnorms, freqs


# parameters
cutoffs = [6.0]
temperatures = [200, 500, 800, 1100, 1400]
size = 6

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

# phonopy parameters
Nq = 51
G2X = get_band(np.array([0, 0, 0]), np.array([0.5, 0.5, 0]), Nq)
X2K2G = get_band(np.array([0.5, 0.5, 1.0]), np.array([0, 0, 0]), Nq)
G2L = get_band(np.array([0, 0, 0]), np.array([0.5, 0.5, 0.5]), Nq)
bands = [G2X, X2K2G, G2L]

# atoms and calculator
prim = bulk('Ni')
atoms = prim.repeat(size)
calc = EMT()

# run scp
cs = ClusterSpace(atoms, cutoffs)
dispersions = dict()
parameters = dict()
for T in temperatures:
    parameters_traj = self_consistent_harmonic_model(
        atoms, calc, cs, T, alpha, n_iterations, n_structures)
    fcp = ForceConstantPotential(cs, parameters_traj[-1])
    dispersions[T] = get_dispersion(prim, fcp)
    parameters[T] = parameters_traj

# plotting
lw = 2.0
fs = 14.0
cmap = plt.get_cmap('viridis')
colors = {T: cmap(i/(len(dispersions)-1)) for i, T in enumerate(temperatures)}

fig = plt.figure(figsize=(14, 4.8))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)

# plot parameter trajectory
T = 1400
dp = [np.linalg.norm(p-p2) for p, p2 in zip(parameters[T], parameters[T][1:])]
ax1.plot(dp, lw=lw)
ax1.set_xlabel('Iteration', fontsize=fs)
ax1.set_ylabel('$\\Delta$ Parameters', fontsize=fs)
ax1.set_title('T={}K'.format(T), fontsize=fs)
ax1.set_xlim([0, n_iterations])
ax1.set_ylim(bottom=0.0)
ax1.tick_params(labelsize=fs)

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

qpts = [0.0, qnorms[0][-1], qnorms[1][-1], qnorms[2][-1]]
qpts_labels = ['$\\Gamma$', 'X', '$\\Gamma$', 'L']
ax2.axvline(x=qpts[1], color='k', linewidth=0.9)
ax2.axvline(x=qpts[2], color='k', linewidth=0.9)

ax2.set_xlabel('Wave vector $\\vec{q}$', fontsize=fs)
ax2.set_ylabel('Frequency (meV)', fontsize=fs)
ax2.set_xticks(qpts, qpts_labels)
ax2.set_xlim([0.0, qnorms[-1][-1]])
ax2.set_ylim([0.0, 50.0])
ax2.legend(loc=1, bbox_to_anchor=(0.45, 0.45), fontsize=fs)
ax2.tick_params(labelsize=fs)

plt.tight_layout()
plt.savefig('self_consistent_phonons.svg')