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.

In this example we look construct EHMs from MD simulations for BCC tantalum. We use cutoffs of 6.0

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 set:

Temperature (K)	Train RMSE (meV/A)
300	41.3
1000	148.7
2000	315.2

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

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

Source code

MD structures are generated in
examples/advanced_topics/temperature_dependent_phonons/3a_run_fcp_md.py

"""
This script carries a series of molecular dynamics simulations at different
temperatures and stores snapshots to file.
"""

import os
import numpy as np
from ase.io import write, read
from ase.io.trajectory import Trajectory

from hiphive import ForceConstantPotential
from hiphive.calculators import ForceConstantCalculator

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
size = 6
temperatures = [2000, 1000, 300]
number_of_equilibration_steps = 3000
number_of_production_steps = 3000
time_step = 1.0
dump_interval = 100

# read FCP
fcp = ForceConstantPotential.read('fcps/fcp_sixth_order.fcp')
prim = fcp.primitive_structure

# setup calculator
atoms_ideal = prim.repeat(size)
fcs = fcp.get_force_constants(atoms_ideal)
calc = ForceConstantCalculator(fcs)

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
    rs = np.random.RandomState(2020)
    MaxwellBoltzmannDistribution(atoms, temperature * units.kB, rng=rs)
    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')

Effective harmonic models are constructed in
examples/advanced_topics/temperature_dependent_phonons/3b_train_ehms.py

from ase.io import read
from hiphive import ClusterSpace, StructureContainer, ForceConstantPotential
from trainstation import Optimizer


temperatures = [2000, 1000, 300]
cutoffs = [6.0]

for T in temperatures:
    structures = read('md_runs/snapshots_T{:}.xyz'.format(T), index=':')

    cs = ClusterSpace(structures[0], cutoffs)
    sc = StructureContainer(cs)
    for s in structures:
        sc.add_structure(s)

    opt = Optimizer(sc.get_fit_data(), train_size=1.0)
    opt.train()
    print(opt)
    fcp = ForceConstantPotential(cs, opt.parameters)
    fcp.write('fcps/ehm_T{}.fcp'.format(T))

Phonon dispersion are plotted using
examples/advanced_topics/temperature_dependent_phonons/4b_plot_ehm_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

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

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

# 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('ehm_phonon_dispersions.svg')