Thermal properties in the harmonic approximation
This section of the tutorial demonstrates how an existing FCP can be employed in conjunction with phonopy to analyze the thermal properties of a material in the harmonic approximation.
Note that this analysis by definition invokes only the second-order force constants and primarily relies on phonopy.
Preparations
First a number of parameters are defined that influence the subsequent analysis. They have been moved to the top for clarity and convenience.
"""
Calculate thermal properties within the harmonic approximation using a
hiPhive force constant potential in combination with phonopy.
Runs in approximately 30 seconds on an Intel Core i5-4670K CPU.
"""
import numpy as np
import matplotlib.pyplot as plt
from ase import Atoms
from ase.build import bulk
from hiphive import ForceConstantPotential
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
# parameters
dim = 5 # dimension in phonopy calculation
Nq = 51 # number of q-points along each segment of the path through the BZ
mesh = [32, 32, 32] # q-point mesh for MSD calculation
temperatures = [300, 600, 900, 1200] # temperatures for evaluating MSD
Next we define the primitive structure (prim) along with a
corresponding PhonopyAtoms object (atoms_phonopy), and
initialize a Phonopy object (phonopy) that represents the
chosen supercell size.
prim = bulk('Ni')
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)
Now we can load the FCP generated previously from a file using the
ForceConstantPotential.read function and retrieve
the supercell (supercell) from the phonopy object, for which we
want to set up the second-order force constant matrix. The latter is
then generated using the
ForceConstantPotential.get_force_constants method.
fcp = ForceConstantPotential.read('fcc-nickel.fcp')
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)
print(fcs)
# access specific parts of the force constant matrix
fcs.print_force_constant((0, 1))
fcs.print_force_constant((10, 12))
At the end of the block, we use the print() function to get a compact view of
the force constant matrix, which yields an output akin to the following block:
=================== ForceConstants ===================
Orders: [2, 3, 4]
Atoms: Atoms(symbols='Ni125', pbc=True, cell=[[0.0, 8.8, 8.8], [8.8, 0.0, 8.8], [8.8, 8.8, 0.0]])
Cluster: (0, 0)
Atom('Ni', [0.0, 0.0, 0.0], index=0)
Atom('Ni', [0.0, 0.0, 0.0], index=0)
force constant:
[[11.55459 0. -0. ]
[ 0. 11.55459 -0. ]
[-0. -0. 11.55459]]
Cluster: (0, 1)
Atom('Ni', [0.0, 0.0, 0.0], index=0)
Atom('Ni', [0.0, 1.7600000000000002, 1.7600000000000002], index=1)
force constant:
[[ 0.04545 0. -0. ]
[ 0. -1.55551 -1.6022 ]
[-0. -1.6022 -1.55551]]
...
Using the ForceConstants.print_cluster function it is furthermore
possible to access the force constants for specific cluster. Here, a
cluster is referred to by the indices of its sites:
Cluster: (10, 12)
Atom('Ni', [3.5200000000000005, 0.0, 3.5200000000000005], index=10)
Atom('Ni', [3.5200000000000005, 3.5200000000000005, 7.040000000000001], index=12)
force constant:
[[ 0.00008 -0. 0. ]
[-0. 0.002 0.00463]
[ 0. 0.00463 0.002 ]]
This concludes the part that involves hiPhive.
Harmonic analysis
Mean-square displacements
Now we can conduct an analysis of the mean-square displacement (and potentially other thermal properties) at several temperatures.
phonopy.set_force_constants(fcs.get_fc_array(order=2))
phonopy.set_mesh(mesh, is_eigenvectors=True, is_mesh_symmetry=False)
phonopy.set_thermal_displacements(temperatures=temperatures)
_, msds = phonopy.get_thermal_displacements()
msds = np.sum(msds, axis=1) # sum up the MSD over x,y,z
for temperature, msd in zip(temperatures, msds):
print('T = {:4d} K MSD = {:.5f} A**2'.format(temperature, msd))
This yields:
T = 300 K MSD = 0.01326 A**2
T = 600 K MSD = 0.02550 A**2
T = 900 K MSD = 0.03797 A**2
T = 1200 K MSD = 0.05049 A**2
In the molecular dynamics (MD) tutorial section these data are compared with MD results.
Phonon dispersion
It is equally straightforward to extract the phonon dispersion. After one has defined a path through the Brillouin zone (BZ) this only requires two lines of code thanks to phonopy.
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)])
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]
# get phonon dispersion
phonopy.set_band_structure(bands)
qvecs, qnorms, freqs, _ = phonopy.get_band_structure()
Finally, we can plot the phonon dispersion using matplotlib.
fig = plt.figure()
kpts = [0.0, qnorms[0][-1], qnorms[1][-1], qnorms[2][-1]]
kpts_labels = ['$\\Gamma$', 'X', '$\\Gamma$', 'L']
plt.axvline(x=kpts[1], color='k', linewidth=0.9)
plt.axvline(x=kpts[2], color='k', linewidth=0.9)
for q, freq, in zip(qnorms, freqs):
plt.plot(q, freq, color='b', linewidth=2.0)
plt.xlabel('Wave vector $\\vec{q}$', fontsize=14.0)
plt.ylabel('Frequency $\\omega$ (THz)', fontsize=14.0)
plt.xticks(kpts, kpts_labels, fontsize=14.0)
plt.xlim([0.0, qnorms[-1][-1]])
plt.ylim([0.0, 12.0])
plt.tight_layout()
plt.savefig('phonon_dispersion.pdf')
This yields the following dispersion, which is virtually indistinguishable from the dispersion one obtains when using the EMT calculator directly.
Phonon dispersion of FCC nickel according to a force constant potential constructed using input data generated using an EMT potential.
Source code
The complete source code is available in
examples/tutorial/3_compute_harmonic_thermal_properties.py
"""
Calculate thermal properties within the harmonic approximation using a
hiPhive force constant potential in combination with phonopy.
Runs in approximately 30 seconds on an Intel Core i5-4670K CPU.
"""
import numpy as np
import matplotlib.pyplot as plt
from ase import Atoms
from ase.build import bulk
from hiphive import ForceConstantPotential
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
# parameters
dim = 5 # dimension in phonopy calculation
Nq = 51 # number of q-points along each segment of the path through the BZ
mesh = [32, 32, 32] # q-point mesh for MSD calculation
temperatures = [300, 600, 900, 1200] # temperatures for evaluating MSD
# set up phonopy
prim = bulk('Ni')
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)
# get force constants
fcp = ForceConstantPotential.read('fcc-nickel.fcp')
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)
print(fcs)
# access specific parts of the force constant matrix
fcs.print_force_constant((0, 1))
fcs.print_force_constant((10, 12))
# get mean square displacements
phonopy.set_force_constants(fcs.get_fc_array(order=2))
phonopy.set_mesh(mesh, is_eigenvectors=True, is_mesh_symmetry=False)
phonopy.set_thermal_displacements(temperatures=temperatures)
_, msds = phonopy.get_thermal_displacements()
msds = np.sum(msds, axis=1) # sum up the MSD over x,y,z
for temperature, msd in zip(temperatures, msds):
print('T = {:4d} K MSD = {:.5f} A**2'.format(temperature, msd))
# FCC q-point paths
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)])
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]
# get phonon dispersion
phonopy.set_band_structure(bands)
qvecs, qnorms, freqs, _ = phonopy.get_band_structure()
# plot dispersion
fig = plt.figure()
kpts = [0.0, qnorms[0][-1], qnorms[1][-1], qnorms[2][-1]]
kpts_labels = ['$\\Gamma$', 'X', '$\\Gamma$', 'L']
plt.axvline(x=kpts[1], color='k', linewidth=0.9)
plt.axvline(x=kpts[2], color='k', linewidth=0.9)
for q, freq, in zip(qnorms, freqs):
plt.plot(q, freq, color='b', linewidth=2.0)
plt.xlabel('Wave vector $\\vec{q}$', fontsize=14.0)
plt.ylabel('Frequency $\\omega$ (THz)', fontsize=14.0)
plt.xticks(kpts, kpts_labels, fontsize=14.0)
plt.xlim([0.0, qnorms[-1][-1]])
plt.ylim([0.0, 12.0])
plt.tight_layout()
plt.savefig('phonon_dispersion.pdf')