Unofficial way to use Phonopy with ASE

If you've ever done any phonon calculations with VASP, then you are most certainly familiar with phonopy. In many scenarios you may want to use interatomic potentials for the forces instead of ab-initio¹ calculated forces. This allows you to more quickly converge and calculate ground-state bandstructures and thermal quantities. You can even look at finite temperature effects with lammps using fix phonon. Phonopy does have a good amount of support for several DFT² and LAMMPS.

I recently have been using atomic simulaton environment (ASE) to leverage some of its calculator interfaces. While ASE does have built in finite-difference method support for phonon calculations, phonopy is more comprehensive. So I wanted to explore how much effort it would take to get the two to "talk", turns out its extremly straightfoward.

The Workflow

First things first, you'll need to run pip install ase and pip install phonopy, its also ideal to install the Brillouin zone path library via pip install seekpath. Now, let me create a structure with ASE and define a calculator, for simplicity I'm going to use the effective medium potential (EMT) calculator thats built into ASE.

from ase.build import bulk
from ase.calculators.emt import EMT
calculator = EMT()
structure = bulk('Cu', 'fcc', a=3.6)

Now to import the phonopy python API tools:

from phonopy import Phonopy
from phonopy.interface.calculator import get_force_sets
from phonopy.structure.atoms import PhonopyAtoms

The import thing to notice is we need to use PhonopyAtoms to create the supercells and displacements. At first this may seem like a big hurdle, in that you'll need a lot of wrapper code, but this is not the case. Simply do:

phnpy_struct = PhonopyAtoms(
    symbols=structure.get_chemical_symbols(),
    positions=structure.get_positions(),
    cell=structure.get_cell(),
    )

Now we construct our phonopy object:

phonons = Phonopy(
    phnpy_struct,
    supercell_matrix=[[5, 0, 0],
                      [0, 5, 0],
                      [0, 0, 5]],
    )

Then we need to specify the displacement magnitude that is used when using the finite-difference method.

phonon.generate_displacements(distance=0.5)

Note

One thing I've noticed is that the distance values to converge the phonon calculations are significantly larger than I've used in the past. Usually displacements of 0.01-0.05 angstroms, but I'm noticing 0.5-0.8 angstroms are need to converge! Seems concerning since we need to remain within the harmonic regime.

Okay straightforward, now what about the forces? How do we get them into our phonopy object phonons. Turns out this is also relatively simple:

import numpy as np
sets_of_forces = []
supercells = phonons.get_supercells_with_displacements()
for i,d in enumerate(supercells):
    # Convert back to ASE
    d_ase = Atoms(symbols=d.get_chemical_symbols(),
        positions=d.get_positions(),
        cell=d.get_cell()
        )
    d_ase.set_calculator(calculator)
    forces = d_ase.get_forces()
    sets_of_forces.append(forces)

Now that we have the forces, fortunately, all we need to do to get the force constants to calculate the phonon properties is:

sets_of_forces = np.array(sets_of_forces)
phonons.forces = sets_of_forces
phonons.produce_force_constants()
phonons.save()

The last command saves the details for phonopy which I believe you can use the outputs to run the standalone phonopy command-line tool. Now we can view the phonon bandstructure:

phonons.phonons.auto_band_structure(plot=True)

which gives:

Phonon bandstructure of Cu FCC using ASE and Phonopy together.

Using the phonopy python api you should also be able to calculate the DOS and thermal properties.

Footnotes

---

1. If the forces on atoms are calculated using the Hellmann-Feynman theorem, then they are referred to as first-principles or ab-initio because they arise from the expectation value of the Hamiltonian about some continous parameter, ex., $\mathbf{F}_i = -\frac{dE_i}{d\mathbf{r}_i} = \langle \psi_{\mathbf{r_i}} \lvert \frac{d\hat{H}_{\mathbf{r}_i}}{{d\mathbf{r}_i}} \rvert \psi_{\mathbf{r_i}} \rangle$. ↩

2. Density functional theory is the primary compute for most high-quality phonon calculations. ↩

Reuse and Attribution