Phonon lifetimes

This section of the tutorial demonstrates how an existing FCP can be employed in conjunction with phono3py to analyze properties that are related to the third-order force constants such as phonon lifetimes (considering only phonon-phonon scattering). Note that this analysis invokes also the second-order force constants and requires a working installation of phono3py.

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 third order properties using a ForceConstantPotential and feeding
the resulting force constants to phono3py.

Runs in approximately 180 seconds on an Intel Core i5-4670K CPU.
"""
import h5py
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from ase import Atoms
from ase.io import write
from ase.build import bulk
from hiphive import ForceConstantPotential
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms

# parameters
dim = 5
mesh = 14
temperatures = [1200]

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 and third-order force constant matrices. The latter are 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)

Since phono3py is controlled via a command line interface the force constant matrices have to be written to file.

fcs.write_to_phonopy('fc2.hdf5')
fcs.write_to_phono3py('fc3.hdf5')
write('POSCAR', prim)

This concludes the part that involves hiPhive.

Lifetimes

To obtain the phonon lifetimes phono3py is called after which the results are read from file.

phono3py_cmd = 'phono3py --dim="{0} {0} {0}" --fc2 --fc3 --br --mesh="'\
               '{1} {1} {1}" --ts="{2}"'.format(
                dim, mesh, ' '.join(str(T) for T in temperatures))
subprocess.call(phono3py_cmd, shell=True)

# collect phono3py data
with h5py.File('kappa-m{0}{0}{0}.hdf5'.format(mesh), 'r') as hf:
    temperatures = hf['temperature'][:]
    frequency = hf['frequency'][:]
    gamma = hf['gamma'][:]

Finally, we can plot the phonon lifetimes versus the respective frequencies using matplotlib.

ms = 4
fs = 14
plt.figure()
plt.plot(frequency.flatten(), gamma[0].flatten(), 'o', ms=ms)

plt.xlabel('Frequency (THz)', fontsize=fs)
plt.ylabel('Gamma (THz)', fontsize=fs)
plt.xlim(left=0)
plt.ylim(bottom=0)

plt.title('T={:d}K'.format(int(temperatures[0])))
plt.gca().tick_params(labelsize=fs)
plt.tight_layout()

plt.savefig('frequency_gamma.pdf')

Source code

The complete source code is available in examples/tutorial/4_compute_third_order_properties.py

"""
Calculate third order properties using a ForceConstantPotential and feeding
the resulting force constants to phono3py.

Runs in approximately 180 seconds on an Intel Core i5-4670K CPU.
"""
import h5py
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from ase import Atoms
from ase.io import write
from ase.build import bulk
from hiphive import ForceConstantPotential
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms

# parameters
dim = 5
mesh = 14
temperatures = [1200]

# get phono3py supercell
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)

# write force constants and POSCAR
fcs.write_to_phonopy('fc2.hdf5')
fcs.write_to_phono3py('fc3.hdf5')
write('POSCAR', prim)

# call phono3py
phono3py_cmd = 'phono3py --dim="{0} {0} {0}" --fc2 --fc3 --br --mesh="'\
               '{1} {1} {1}" --ts="{2}"'.format(
                dim, mesh, ' '.join(str(T) for T in temperatures))
subprocess.call(phono3py_cmd, shell=True)

# collect phono3py data
with h5py.File('kappa-m{0}{0}{0}.hdf5'.format(mesh), 'r') as hf:
    temperatures = hf['temperature'][:]
    frequency = hf['frequency'][:]
    gamma = hf['gamma'][:]

# generate plot of lifetimes
ms = 4
fs = 14
plt.figure()
plt.plot(frequency.flatten(), gamma[0].flatten(), 'o', ms=ms)

plt.xlabel('Frequency (THz)', fontsize=fs)
plt.ylabel('Gamma (THz)', fontsize=fs)
plt.xlim(left=0)
plt.ylim(bottom=0)

plt.title('T={:d}K'.format(int(temperatures[0])))
plt.gca().tick_params(labelsize=fs)
plt.tight_layout()

plt.savefig('frequency_gamma.pdf')