Phonopy API for Python#

This is under development. Configurations may alter. Requests or suggestions are very welcome.

Import modules#

After setting the phonopy python path, the phonopy module is imported by:

from phonopy import Phonopy

Crystal structure is defined by the PhonopyAtoms class. The PhonopyAtoms module is imported by:

from phonopy.structure.atoms import PhonopyAtoms

The instance of PhonopyAtoms can be made by reading a crystal structure in a variety of calculator formats found at Interfaces to calculators.

from phonopy.interface.calculator import read_crystal_structure
unitcell, _ = read_crystal_structure("POSCAR-unitcell", interface_mode='vasp')

For VASP format, the keyword argument of interface_mode can be omitted. For QE,

unitcell, optional_structure_info = read_crystal_structure("NaCl.in", interface_mode='qe')

Note that read_crystal_structure returns a tuple and the first element is the PhonopyAtoms instance.

Work flow#

The work flow is schematically shown in Work flow.

Pre-process#

The first step is to create a Phonopy object with at least two arguments, a unit cell (PhonopyAtoms object, see PhonopyAtoms class) and a supercell matrix (3x3 array, see Supercell matrix). In the following example, a \(2\times 2\times 2\) supercell is created. The displacements to be introduced to the supercell are internally generated by the generate_displacements() method with the distance keyword argument. The supercells with displacements are obtained by get_supercells_with_displacements() method as a list of PhonopyAtoms objects.

import numpy as np
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms

a = 5.404
unitcell = PhonopyAtoms(symbols=['Si'] * 8,
                        cell=(np.eye(3) * a),
                        scaled_positions=[[0, 0, 0],
                                          [0, 0.5, 0.5],
                                          [0.5, 0, 0.5],
                                          [0.5, 0.5, 0],
                                          [0.25, 0.25, 0.25],
                                          [0.25, 0.75, 0.75],
                                          [0.75, 0.25, 0.75],
                                          [0.75, 0.75, 0.25]])
phonon = Phonopy(unitcell,
                 supercell_matrix=[[2, 0, 0], [0, 2, 0], [0, 0, 2]],
                 primitive_matrix=[[0, 0.5, 0.5],
                                   [0.5, 0, 0.5],
                                   [0.5, 0.5, 0]])
phonon.generate_displacements(distance=0.03)
supercells = phonon.supercells_with_displacements

In this example, the displacement distance is set to 0.03 (in Angstrom if the crystal structure uses the Angstrom unit and the default value is 0.01.) The supercells with displacements are given as a list of PhonopyAtoms. See Read and write crystal structures to write those into files in a crystal structure format.

The frequency unit conversion factor to THz has to be set by using the factor keyword in Phonopy class. The factors are VaspToTHz for VASP, Wien2kToTHz for Wien2k, AbinitToTHz for Abinit, PwscfToTHz for Pwscf, ElkToTHz for Elk, SiestaToTHz for Siesta, CrystalToTHz for CRYSTAL, FleurToTHz for Fleur, VaspToTHz, and DftbpToTHz for DFTB+ is the default value. For example:

from phonopy.units import AbinitToTHz

phonon = Phonopy(unitcell,
                 supercell_matrix=[[2, 0, 0], [0, 2, 0], [0, 0, 2]],
                 primitive_matrix=[[0, 0.5, 0.5],
                                   [0.5, 0, 0.5],
                                   [0.5, 0.5, 0]],
                 factor=AbinitToTHz)

Some more information on physical unit conversion is found at FREQUENCY_CONVERSION_FACTOR, Physical unit conversion, and Interfaces to calculators.

Post process#

Forces on atoms are supposed to be obtained by running force calculator (e.g. VASP) with each supercell with a displacement. Then the forces in the calculation outputs have to be collected by users. However output parsers for selected calculators are found under phonopy.interface, which may be useful. The forces have to be stored in a specific structure: a numpy array (or nested list) as follows:

[ [ [ f_1x, f_1y, f_1z ], [ f_2x, f_2y, f_2z ], ... ], # first supercell
  [ [ f_1x, f_1y, f_1z ], [ f_2x, f_2y, f_2z ], ... ], # second supercell
  ...                                                   ]

This array (sets_of_forces) is set to the Phonopy object by:

phonon.forces = sets_of_forces

This is the case when the set of atomic displacements is generated internally. The information of displacements is already stored in the Phonopy object. But if you want to input the forces together with the corresponding custom set of displacements, displacement_dataset has to be prepared as a python dictionary as follows:

displacement_dataset =
   {'natom': number_of_atoms_in_supercell,
    'first_atoms': [
      {'number': atom index of displaced atom (starting with 0),
       'displacement': displacement in Cartesian coordinates,
       'forces': forces on atoms in supercell},
      {...}, ...]}

This is set to the Phonopy object by:

phonopy.dataset = displacement_dataset

From the set of displacements and forces, force constants internally with calculated supercell sets of forces by

phonon.produce_force_constants()

If you have force constants and don’t need to create force constants from forces and displacements, simply set your force constants by

phonon.force_constants = force_constants

The force constants matrix is given in 4 dimensional array (better to be a numpy array of dtype='double', order='C'). The shape of force constants matrix is (N, N, 3, 3) where N is the number of atoms in the supercell and 3 gives Cartesian axes. The compact force constants matrix with (Np, N, 3, 3) where Np is the number of atoms in the primitive cell is also supported. See the details at FORCE_CONSTANTS and force_constants.hdf5.

Phonon calculation#

Save parameters (phonopy.save)#

Basic information and parameters needed for phonon calculation are saved into a file by phonopy.save.

phonon.save()

The default file name is phonopy_params.yaml. Force sets, displacements, Born effective charges, and dielectric constant are written in the default behaviour. If force constants are needed to be written in the yaml file, the argument settings is set as follows:

phonon.save(settings={'force_constants': True})

Band structure#

Set band paths (run_band_structure()) and get the results (get_band_structure_dict()).

A dictionary with qpoints, distances, frequencies, eigenvectors, group_velocities is returned by get_band_structure_dict(). Eigenvectors can be obtained when with_eigenvectors=True at run_band_structure(). See the details at docstring of Phonopy.get_band_structure_dict. Phonon frequency is sqrt(eigenvalue). A negative eigenvalue has to correspond to the imaginary frequency, but for the plotting, it is set as the negative value in the above example. In addition, you need to multiply by your unit conversion factor. In the case of VASP to transform to THz, the factor is 15.633302.

In example/NaCl, the phonopy is executed from python script, e.g.,

import phonopy
from phonopy.phonon.band_structure import get_band_qpoints_and_path_connections

path = [[[0, 0, 0], [0.5, 0, 0.5], [0.625, 0.25, 0.625]],
        [[0.375, 0.375, 0.75], [0, 0, 0], [0.5, 0.5, 0.5], [0.5, 0.25, 0.75]]]
labels = ["$\\Gamma$", "X", "U", "K", "$\\Gamma$", "L", "W"]
qpoints, connections = get_band_qpoints_and_path_connections(path, npoints=51)
phonon = phonopy.load("phonopy_disp.yaml")
phonon.run_band_structure(qpoints, path_connections=connections, labels=labels)
phonon.plot_band_structure().show()

# To plot DOS next to band structure
phonon.run_mesh([20, 20, 20])
phonon.run_total_dos()
phonon.plot_band_structure_and_dos().show()

# To plot PDOS next to band structure
phonon.run_mesh([20, 20, 20], with_eigenvectors=True, is_mesh_symmetry=False)
phonon.run_projected_dos()
phonon.plot_band_structure_and_dos(pdos_indices=[[0], [1]]).show()

path_connections and labels are unnecessary to set unless nice looking plotting is needed. To obtain eigenvectors, it is necessary to inform to store eigenvectors by:

phonon.run_band_structure(bands, with_eigenvectors=True)

To obtain group velocities:

phonon.run_band_structure(bands, with_group_velocities=True)

Automatic selection of band paths using SeeK-path is invoked by

phonon.auto_band_structure()

and to plot

phonon.auto_band_structure(plot=True).show()

To use this method, seekpath python module is needed.

Mesh sampling#

Set sampling mesh (set_mesh) in reciprocal space. The irreducible q-points and corresponding q-point weights, eigenvalues, and eigenvectors are obtained by get_mesh_dict(). mesh gives the sampling mesh with Monkhorst-Pack scheme. The keyword shift gives the fractional mesh shift with respect to the neighboring grid points.

mesh = [20, 20, 20]
phonon.run_mesh(mesh)
mesh_dict = phonon.get_mesh_dict()
qpoints = mesh_dict['qpoints']
weights = mesh_dict['weights']
frequencies = mesh_dict['frequencies']
eigenvectors = mesh_dict['eigenvectors']
group_velocities = mesh_dict['group_velocities']

To obtain eigenvectors, it is necessary to inform to store eigenvectors by:

phonon.run_mesh([20, 20, 20], with_eigenvectors=True)

and for group velocities:

phonon.run_mesh([20, 20, 20], with_group_velocities=True)

The first argument of run_mesh() can be a float value, which is a length measure as explained at MESH, MP, or MESH_NUMBERS, for example:

phonon.run_mesh(100.0)

DOS and PDOS#

Before starting mesh sampling has to be finished. Then set parameters (run_total_dos() or run_projected_dos()) and write the results into files (write_total_dos() and write_projected_dos()). In the case of PDOS, the eigenvectors have to be calculated in the mesh sampling. To get the results get_total_dos_dict() and get_projected_dos_dict() can be used.

To plot total DOS,

phonon.run_mesh([20, 20, 20]) phonon.run_total_dos()
phonon.plot_total_dos().show()

and projected DOS

phonon.run_mesh([20, 20, 20], with_eigenvectors=True, is_mesh_symmetry=False)
phonon.run_projected_dos()
phonon.plot_projected_dos().show()

Convenient shortcuts exist as follows:

phonon.auto_total_dos(plot=True).show()

and

phonon.auto_projected_dos(plot=True).show()

Thermal properties#

Before starting the thermal property calculation, the mesh sampling calculation has to be done in the THz unit. The unit conversion factor for phonon frequency is set in the pre-process of Phonopy with the factor keyword. Calculation range of temperature is set by the parameters run_thermal_properties. Helmholtz free energy, entropy, heat capacity at constant volume at temperatures are obtained by get_thermal_properties_dict, where the results are given as a dictionary of temperatures, Helmholtz free energy, entropy, and heat capacity with keys temperatures, free_energy, entropy, and heat_capacity, respectively.

phonon.run_mesh([20, 20, 20])
phonon.run_thermal_properties(t_step=10,
                              t_max=1000,
                              t_min=0)
tp_dict = phonon.get_thermal_properties_dict()
temperatures = tp_dict['temperatures']
free_energy = tp_dict['free_energy']
entropy = tp_dict['entropy']
heat_capacity = tp_dict['heat_capacity']

for t, F, S, cv in zip(temperatures, free_energy, entropy, heat_capacity):
    print(("%12.3f " + "%15.7f" * 3) % ( t, F, S, cv ))

phonon.plot_thermal_properties().show()

Non-analytical term correction#

To apply non-analytical term correction, Born effective charge tensors for all atoms in primitive cell, dielectric constant tensor, and the unit conversion factor have to be correctly set. The tensors are given in Cartesian coordinates.

born = [[[1.08878299, 0, 0],
         [0, 1.08878299, 0],
         [0, 0, 1.08878299]],
        [[-1.08878299, 0, 0],
         [0, -1.08878299, 0],
         [0, 0, -1.08878299]]]
epsilon = [[2.56544559, 0, 0],
           [0, 2.56544559, 0],
           [0, 0, 2.56544559]]
factors = 14.400
phonon.nac_params = {'born': born,
                     'factor': factors,
                     'dielectric': epsilon}

Data structure#

Eigenvectors#

Eigenvectors are given as the column vectors. Internally phonopy uses numpy.linalg.eigh and eigh is a wrapper of LAPACK. So eigenvectors follow the convention of LAPACK, which can be shown at http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigh.html

Eigenvectors corresponding to phonopy yaml output are obtained as follows.

Band structure#

if eigvecs is not None:
    for eigvecs_on_path in eigvecs:
        for eigvecs_at_q in eigvecs_on_path:
            for vec in eigvecs_at_q.T:
                print(vec)

Mesh sampling#

if eigvecs is not None:
    for eigvecs_at_q in eigvecs:
        for vec in eigvecs_at_q.T:
            print(vec)

PhonopyAtoms class#

Initialization#

The usable keywords in the initialization are:

cell=None,
scaled_positions=None,
positions=None,
numbers=None,
symbols=None,
masses=None,
magnetic_moments=None,

At least three arguments have to be given at the initialization, which are

  • cell

  • positions or scaled_positions

  • symbols or numbers

Variables#

The following variables are implemented in the PhonopyAtoms class in phonopy/structure/atoms.py.

cell#

Basis vectors are given in the matrix form in Cartesian coordinates.

[ [ a_x, a_y, a_z ], [ b_x, b_y, b_z ], [ c_x, c_y, c_z ] ]

scaled_positions#

Atomic positions in fractional coordinates.

[ [ x1_a, x1_b, x1_c ], [ x2_a, x2_b, x2_c ], [ x3_a, x3_b, x3_c ], ... ]

positions#

Cartesian positions of atoms.

positions = np.dot(scaled_positions, cell)

where np means the numpy module (import numpy as np).

symbols#

Chemical symbols, e.g.,

['Zn', 'Zn', 'O', 'O']

for the ZnO unit cell.

numbers#

Atomic numbers, e.g.,

[30, 30, 8, 8]

for the ZnO unit cell.

masses#

Atomic masses, e.g.,

[65.38, 65.38, 15.9994, 15.9994]

for the ZnO unit cell.

Attributes#

cell
positions
scaled_positions
masses
magnetic_moments
symbols
numbers
volume

where volume is the getter only.

Methods#

unitcell.get_number_of_atoms() is equivalent to len(unitcell). An instance can be deep-copied by unitcell.copy(). Human-readable crystal structure in Yaml format is shown by print(unitcell). unitcell.to_tuple converts to spglib crystal structure (https://spglib.github.io/spglib/python-spglib.html#crystal-structure-cell).

Definitions of variables#

Primitive matrix#

Primitive matrix \(M_\mathrm{p}\) is a tranformation matrix from lattice vectors to those of a primitive cell if there exists the primitive cell in the lattice vectors. Following a crystallography convention, the transformation is given by

\[( \mathbf{a}_\mathrm{p} \; \mathbf{b}_\mathrm{p} \; \mathbf{c}_\mathrm{p} ) = ( \mathbf{a}_\mathrm{u} \; \mathbf{b}_\mathrm{u} \; \mathbf{c}_\mathrm{u} ) M_\mathrm{p}\]

where \(\mathbf{a}_\mathrm{u}\), \(\mathbf{b}_\mathrm{u}\), and \(\mathbf{c}_\mathrm{u}\) are the column vectors of the original lattice vectors, and \(\mathbf{a}_\mathrm{p}\), \(\mathbf{b}_\mathrm{p}\), and \(\mathbf{c}_\mathrm{p}\) are the column vectors of the primitive lattice vectors. Be careful that the lattice vectors of the PhonopyAtoms class are the row vectors (cell). Therefore the phonopy code, which relies on the PhonopyAtoms class, is usually written such as

primitive_lattice = np.dot(original_lattice.T, primitive_matrix).T,

or equivalently,

primitive_lattice = np.dot(primitive_matrix.T, original_lattice)

Supercell matrix#

Supercell matrix \(M_\mathrm{s}\) is a transformation matrix from lattice vectors to those of a super cell. Following a crystallography convention, the transformation is given by

\[( \mathbf{a}_\mathrm{s} \; \mathbf{b}_\mathrm{s} \; \mathbf{c}_\mathrm{s} ) = ( \mathbf{a}_\mathrm{u} \; \mathbf{b}_\mathrm{u} \; \mathbf{c}_\mathrm{u} ) M_\mathrm{s}\]

where \(\mathbf{a}_\mathrm{u}\), \(\mathbf{b}_\mathrm{u}\), and \(\mathbf{c}_\mathrm{u}\) are the column vectors of the original lattice vectors, and \(\mathbf{a}_\mathrm{s}\), \(\mathbf{b}_\mathrm{s}\), and \(\mathbf{c}_\mathrm{s}\) are the column vectors of the supercell lattice vectors. Be careful that the lattice vectors of the PhonopyAtoms class are the row vectors (cell). Therefore the phonopy code, which relies on the PhonopyAtoms class, is usually written such as

supercell_lattice = np.dot(original_lattice.T, supercell_matrix).T,

or equivalently,

supercell_lattice = np.dot(supercell_matrix.T, original_lattice)

Symmetry search tolerance#

Symmetry search tolerance (often the name symprec is used in phonopy) is used to determine symmetry operations of the crystal structures. The physical unit follows that of input crystal structure.

Load phonopy settings phonopy.load#

phonopy.load is a convenient function for creating a Phonopy instance by loading data from a phonopy_xxx.yaml file, which may contain the necessary information to run phonopy. A typical usage example is:

phonon = phonopy.load("phonopy_params.yaml")

To avoid producing force constants,

phonon = phonopy.load("phonopy_params.yaml", produce_fc=False)

The details are found in the docstring:

In [1]: import phonopy
In [2]: help(phonopy.load)

More detailed control can be performed using phonopy.load.

phonon = phonopy.load(supercell_matrix=[2, 2, 2],
                      primitive_matrix='auto',
                      unitcell_filename="POSCAR",
                      force_constants_filename="force_constants.hdf5")

With is_nac=True (default), BORN file in the current directory is read if it exists. If supercell is passed and primitive matrix and supercell_matrix are not set, the primitive cell is automatically searched:

phonon = phonopy.load(supercell_filename="SPOSCAR",
                      force_constants_filename="force_constants.hdf5")

If FORCE_SETS exists in the current directory, this below works to be ready for post-process calculation with automatic choice of primitive matrix:

phonon = phonopy.load(supercell_filename="SPOSCAR")

For example, in the example/NaCl directory, phonon band structure of NaCl is easily plotted by

In [1]: import phonopy

In [2]: ph = phonopy.load(supercell_filename="SPOSCAR", log_level=1)
Supercell structure was read from "SPOSCAR".
NAC params were read from "BORN".
Force constants were read from "FORCE_CONSTANTS".

In [3]: print(ph.primitive)
lattice:
- [     0.000000000000000,     2.845150738087836,     2.845150738087836 ] # a
- [     2.845150738087836,     0.000000000000000,     2.845150738087836 ] # b
- [     2.845150738087836,     2.845150738087836,     0.000000000000000 ] # c
points:
- symbol: Na # 1
  coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
  mass: 22.989769
- symbol: Cl # 2
  coordinates: [  0.500000000000000,  0.500000000000000,  0.500000000000000 ]
  mass: 35.453000

In [4]: ph.nac_params
Out[4]:
{'born': array([[[ 1.08703000e+00, -5.17677526e-34, -1.06309751e-33],
         [-5.45419984e-34,  1.08703000e+00,  1.06309751e-33],
         [ 0.00000000e+00,  3.08148791e-33,  1.08703000e+00]],

        [[-1.08672000e+00, -2.93244455e-35,  5.15939995e-34],
         [ 5.45264441e-34, -1.08672000e+00, -5.15939995e-34],
         [ 0.00000000e+00,  0.00000000e+00, -1.08672000e+00]]]),
 'factor': 14.4,
 'dielectric': array([[2.43533967, 0.        , 0.        ],
        [0.        , 2.43533967, 0.        ],
        [0.        , 0.        , 2.43533967]])}

In [5]: ph.auto_band_structure(plot=True).show()

Read and write crystal structures#

There is a function to write the PhonopyAtoms instance into crystal structure formats of different force calculators, write_crystal_structure. This works as a partner of read_crystal_structure. Taking an example of QE interface, how to use these functions is shown below.

In [1]: from phonopy.interface.calculator import read_crystal_structure, write_crystal_structure

In [2]: !cat "NaCl.in"
 &control
    calculation = 'scf'
    tprnfor = .true.
    tstress = .true.
    pseudo_dir = '/home/togo/espresso/pseudo/'
 /
 &system
    ibrav = 0
    nat = 8
    ntyp = 2
    ecutwfc = 70.0
 /
 &electrons
    diagonalization = 'david'
    conv_thr = 1.0d-9
 /
ATOMIC_SPECIES
 Na  22.98976928 Na.pbe-spn-kjpaw_psl.0.2.UPF
 Cl  35.453      Cl.pbe-n-kjpaw_psl.0.1.UPF
ATOMIC_POSITIONS crystal
 Na   0.0000000000000000  0.0000000000000000  0.0000000000000000
 Na   0.0000000000000000  0.5000000000000000  0.5000000000000000
 Na   0.5000000000000000  0.0000000000000000  0.5000000000000000
 Na   0.5000000000000000  0.5000000000000000  0.0000000000000000
 Cl   0.5000000000000000  0.5000000000000000  0.5000000000000000
 Cl   0.5000000000000000  0.0000000000000000  0.0000000000000000
 Cl   0.0000000000000000  0.5000000000000000  0.0000000000000000
 Cl   0.0000000000000000  0.0000000000000000  0.5000000000000000
CELL_PARAMETERS angstrom
 5.6903014761756712 0 0
 0 5.6903014761756712 0
 0 0 5.6903014761756712
K_POINTS automatic
 8 8 8 1 1 1

In [3]: cell, optional_structure_info = read_crystal_structure("NaCl.in", interface_mode='qe')

In [4]: optional_structure_info
Out[4]:
('NaCl.in',
 {'Na': 'Na.pbe-spn-kjpaw_psl.0.2.UPF', 'Cl': 'Cl.pbe-n-kjpaw_psl.0.1.UPF'})

In [5]: write_crystal_structure("NaCl-out.in", cell, interface_mode='qe', optional_structure_info=optional_structure_info)

In [6]: !cat "NaCl-out.in"
!    ibrav = 0, nat = 8, ntyp = 2
CELL_PARAMETERS bohr
   10.7531114272216008    0.0000000000000000    0.0000000000000000
    0.0000000000000000   10.7531114272216008    0.0000000000000000
    0.0000000000000000    0.0000000000000000   10.7531114272216008
ATOMIC_SPECIES
 Na   22.98977   Na.pbe-spn-kjpaw_psl.0.2.UPF
 Cl   35.45300   Cl.pbe-n-kjpaw_psl.0.1.UPF
ATOMIC_POSITIONS crystal
 Na   0.0000000000000000  0.0000000000000000  0.0000000000000000
 Na   0.0000000000000000  0.5000000000000000  0.5000000000000000
 Na   0.5000000000000000  0.0000000000000000  0.5000000000000000
 Na   0.5000000000000000  0.5000000000000000  0.0000000000000000
 Cl   0.5000000000000000  0.5000000000000000  0.5000000000000000
 Cl   0.5000000000000000  0.0000000000000000  0.0000000000000000
 Cl   0.0000000000000000  0.5000000000000000  0.0000000000000000
 Cl   0.0000000000000000  0.0000000000000000  0.5000000000000000

Depending on calculator interfaces, all the information can not be recovered from the information obtained from read_crystal_structure. More details about how write_crystal_structure works may need to read directly the code.

Getting parameters for non-analytical term correction#

Parameters for non-analytical term correction may be made as follows. This example assumes that the user knows what are the unit cell and primitive cell and that the Born effective charge and dielectric constant were calculated using VASP code by the unit cell.

import io
import numpy as np
from phonopy.units import Hartree, Bohr
from phonopy.structure.symmetry import symmetrize_borns_and_epsilon
from phonopy.interface.vasp import VasprunxmlExpat

with io.open("vasprun.xml", "rb") as f:
    vasprun = VasprunxmlExpat(f)
    vasprun.parse():
    epsilon = vasprun.epsilon
    borns = vasprun.born
    unitcell = vasprun.cell

borns_, epsilon_ = symmetrize_borns_and_epsilon(
    borns,
    epsilon,
    unitcell,
    primitive_matrix=[[0, 0.5, 0.5],
                      [0.5, 0, 0.5],
                      [0.5, 0.5, 0]],
    supercell_matrix=np.diag([2, 2, 2]),
    symprec=1e-5)

nac_params = {'born': borns_,
              'factor': Hartree * Bohr,
              'dielectric': epsilon_}