Effective harmonic models

Using hiPhive it is straightforward to generate effective harmonic models (EHMs) from molecular dynamics (MD) trajectories. These models enable one to describe one for example the temperature dependence of the phonon dispersion as well as derived quantities such as the harmonic free energy.

This example comprises two scripts: The first one carries out a series of MD simulations using the EMT calculator from ASE, while the second scripts generates EHMs for each temperature and analyzes the resulting phonon dispersions.

The EHMs constructed here include interactions up to the second neighbor shell as can be seen from the list of orbits:

===================================== List of Orbits =====================================
index | order |      elements      |  radius  |     prototype      | clusters | parameters
------------------------------------------------------------------------------------------
  0   |   2   |       Ni Ni        |  0.0000  |       (0, 0)       |    1     |     1
  1   |   2   |       Ni Ni        |  1.2445  |       (0, 1)       |    6     |     3
  2   |   2   |       Ni Ni        |  2.4890  |       (0, 2)       |    6     |     3
  3   |   2   |       Ni Ni        |  2.1556  |       (0, 3)       |    12    |     4
  4   |   2   |       Ni Ni        |  3.0484  |      (0, 10)       |    4     |     2
  5   |   2   |       Ni Ni        |  2.7828  |      (0, 11)       |    12    |     4
  6   |   2   |       Ni Ni        |  1.7600  |      (0, 16)       |    3     |     2
==========================================================================================

With increasing temperature the ability of an EHM to map the forces present in the structure decreases as evident from the root-mean-square errors over the training and test sets:

Temperature (K) Training RMSE (meV/A) Test RMSE (meV/A)
200 74.6 75.3
800 289.8 272.8
1400 525.8 497.2

The phonon dispersions obtained in this fashion are shown in the following figure.

../_images/phonon_dispersions_from_effective_harmonic_models.svg

Phonon dispersions for FCC Ni from a series of effective harmonic models generated from MD trajectories recorded at different temperatures.

Source code

Rattled structures are generated in
tutorial/advanced/effective_harmonic_models/prepare_data.py
'''
This script carries a series of molecular dynamics simulations at different
temperatures and stores snapshots to file.

This scripts runs in approximately 4 minutes on an Intel Core i7-6700K CPU.
'''

import os
from ase.io import write, read
from ase.build import bulk
from ase.calculators.emt import EMT
from ase.io.trajectory import Trajectory

from ase import units
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.langevin import Langevin
from ase.md import MDLogger

if not os.path.exists('md_runs'):
    os.makedirs('md_runs')


# parameters
temperatures = [200, 800, 1400]
number_of_equilibration_steps = 500
number_of_production_steps = 500
time_step = 5
dump_interval = 50

# setup system
atoms_ideal = bulk('Ni').repeat(5)
calc = EMT()

for temperature in temperatures:

    print('Temperature: {}'.format(temperature))

    # set up molecular dynamics simulation
    atoms = atoms_ideal.copy()
    atoms.set_calculator(calc)
    dyn = Langevin(atoms, time_step * units.fs, temperature * units.kB, 0.02)

    # equilibration run
    MaxwellBoltzmannDistribution(atoms, temperature * units.kB)
    dyn.run(number_of_equilibration_steps)

    # production run
    log_file = 'md_runs/T{:}.log'.format(temperature)
    traj_file = 'md_runs/T{:}.traj'.format(temperature)
    logger = MDLogger(dyn, atoms, log_file,
                      header=True, stress=False,
                      peratom=True, mode='w')
    traj_writer = Trajectory(traj_file, 'w', atoms)
    dyn.attach(logger, interval=dump_interval)
    dyn.attach(traj_writer.write, interval=dump_interval)
    dyn.run(number_of_production_steps)

    # prepare snapshots for later use
    frames = []
    for atoms in read(traj_file, ':'):
        forces = atoms.get_forces()
        displacements = atoms.positions - atoms_ideal.get_positions()
        atoms.positions = atoms_ideal.get_positions()
        atoms.new_array('displacements', displacements)
        atoms.new_array('forces', forces)
        frames.append(atoms.copy())
    print(' Number of snapshots: {}'.format(len(frames)))
    write('md_runs/snapshots_T{:}.xyz'.format(temperature),
          frames, format='extxyz')
The force distributions are analyzed in
tutorial/advanced/effective_harmonic_models/generate_harmonic_models.py
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.cm as cmx

from ase import Atoms
from ase.build import bulk
from ase.io import read

from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms

from hiphive import ClusterSpace, StructureContainer, ForceConstantPotential
from hiphive.fitting import Optimizer


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


# parameters
THz2meV = 4.13567
cutoffs = [6.5]
temperatures = range(200, 1401, 600)
dim = 5
prim = bulk('Ni')
# ------------------------------

# set up phonon dispersion
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]

# set up ClusterSpace and StructureContainer
cs = ClusterSpace(prim, cutoffs)
print(cs)
cs.print_orbits()

sc = StructureContainer(cs)
user_tag = 'md-T{}'
for temperature in temperatures:
    structures = read('md_runs/snapshots_T{:}.xyz@:'.format(temperature))
    for structure in structures:
        sc.add_structure(structure, user_tag=user_tag.format(temperature))
print(sc)

band_structures = []
for temperature in temperatures:

    # Construct ForceConstantPotential
    structure_indices = sc.get_structure_indices(user_tag.format(temperature))
    opt = Optimizer(sc.get_fit_data(structure_indices))
    opt.train()
    print(opt)
    fcp = ForceConstantPotential(cs, opt.parameters)

    # 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=dim*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()
    band_structures.append([temperature, qnorms, freqs])

# -----------------------------------------

fig = plt.figure()
lw = 2.0
fs = 14.0

kpts = [0.0, qnorms[0][-1], qnorms[1][-1], qnorms[2][-1]]
kpts_labels = ['$\\Gamma$', 'X', '$\\Gamma$', 'L']

cm = plt.get_cmap('viridis')
cNorm = mcolors.Normalize(vmin=0, vmax=1)
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)
colors = [scalarMap.to_rgba(float(k) / (len(band_structures) - 1))
          for k in range(len(band_structures))]

plt.axvline(x=kpts[1], color='k', linewidth=0.9)
plt.axvline(x=kpts[2], color='k', linewidth=0.9)

for color, (T, qnorms, freqs) in zip(colors, band_structures):

    for q, freq, in zip(qnorms, freqs):
        plt.plot(q, freq * THz2meV, color=color, lw=lw)
    plt.plot(np.nan, np.nan, color=color, lw=lw, label='T = {}'.format(T))

plt.xlabel('Wave vector $\\vec{q}$', fontsize=fs)
plt.ylabel('$\omega$ (meV)', fontsize=fs)
plt.xticks(kpts, kpts_labels, fontsize=fs)
plt.xlim([0.0, qnorms[-1][-1]])
plt.ylim([0.0, 50.0])

plt.legend()
plt.tight_layout()
plt.savefig('phonon_dispersions_from_effective_harmonic_models.pdf')