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.
We use cutoffs of 6.5 and 3.0 Å for second and third order, respectively. Note that the third order term is only included to reduce the noise in the fit. These cutoffs yields the following 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
7 | 3 | Ni Ni Ni | 1.1062 | (0, 0, 1) | 12 | 5
8 | 3 | Ni Ni Ni | 1.4370 | (0, 1, 5) | 8 | 7
==========================================================================================
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)
200
11.3
800
94.9
1400
221.1
The phonon dispersions obtained in this fashion are shown in the following figure.
Phonon dispersions for FCC Ni 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/effective_harmonic_models/1_run_md_simulations.py
"""
This script carries a series of molecular dynamics simulations at different
temperatures and stores snapshots to file.
This scripts runs in approximately 6 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 = 20
# 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')
Effective harmonic models are constructed and their phonon dispersions plotted in
examples/advanced_topics/effective_harmonic_models/2_generate_ehm_dispersions.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, 3.0]
temperatures = [200, 800, 1400]
dim = 5
prim = bulk('Ni')
# set up phonon band paths
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)
# construct effective harmonic dispersions
band_structures = []
for temperature in temperatures:
# Construct ForceConstantPotential
structure_indices = [i for i, fs in enumerate(sc)
if fs.user_tag == user_tag.format(temperature)]
opt = Optimizer(sc.get_fit_data(structure_indices), train_size=1.0)
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])
# plot dispersions
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_ehm.svg')