Phonopy API for Python

Phonopy API for Python#

Note

The Python API is being restructured toward phonopy v5 and v6; see Development for the architecture principles, migration plan, and deprecation policy. Two points are directly relevant to API users:

Mutating a calculation input (force_constants, nac_params, masses, the displacement dataset, forces, …) invalidates all results derived from it. Run the corresponding run_* method again before accessing the result; accessing a cleared result raises RuntimeError.
APIs planned for removal emit DeprecationWarning. New code should not rely on them.

Three unit cells#

In the Phonopy class, mainly three different unit cells are used, unitcell, supercell, and primitive, which can be accessed with these attributes in the instance. unitcell is the initial input crystal structure. From unitcell, supercell is created by supercell_matrix (Supercell matrix). Then primitive is created from supercell with supercell_matrix and primitive_matrix (Primitive matrix).

Having the supercell_matrix (\(\mathrm{M}_\mathrm{s}\)) and the primitive matrix (\(\mathrm{M}_\mathrm{p}\)), for example,

\[\begin{split}\mathrm{M}_\mathrm{s} = \begin{pmatrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{pmatrix} \;\text{and}\;\; \mathrm{M}_\mathrm{p} = \begin{pmatrix} 0 & 1/2 & 1/2 \\ 1/2 & 0 & 1/2 \\ 1/2 & 1/2 & 0 \end{pmatrix},\end{split}\]

respectively, the basis vectors of supercell are built from those of unitcell 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} ) \mathrm{M}_\mathrm{s}\]

and the basis vectors of primitive are built from those of supercell by

\[( \mathbf{a}_\mathrm{p} \; \mathbf{b}_\mathrm{p} \; \mathbf{c}_\mathrm{p} ) = ( \mathbf{a}_\mathrm{s} \; \mathbf{b}_\mathrm{s} \; \mathbf{c}_\mathrm{s} ) \mathrm{M}_\mathrm{s}^{-1} \mathrm{M}_\mathrm{p}.\]

Once supercell and primitive are made, unitcell will not be used basically.

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 (0.01 by default). The units are determined by the calculator used (See Physical unit system for calculator for the complete list). This example uses the default calculator for the instance of Phonopy—VASP—which uses Angstroms for distance. 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.

Calculators#

If not using the default calculator ("vasp"), the calculator keyword argument must also be set in your instance of Phonopy (e.g. Phonopy(..., calculator="qe")).

The range of supported calculators use different units internally. If not using VASP, set set_factor_by_calculator to True for correct unit conversion. Some more information on physical unit conversion is found at FREQUENCY_CONVERSION_FACTOR, Physical unit conversion, and Interfaces to calculators.

Pre-processing Example#

import numpy as np
from phonopy import Phonopy
from phonopy.interface.calculator import (
    read_crystal_structure,
    write_crystal_structure,
)

calc = "qe"  # Quantum Espresso
unitcell, optional_structure_info = read_crystal_structure("pw.in", interface_mode=calc)

phonon = Phonopy(
    unitcell, supercell_matrix=np.eye(3), calculator=calc, set_factor_by_calculator=True
)

phonon.generate_displacements(distance=0.03)
supercells = phonon.supercells_with_displacements

for i, supercell in enumerate(supercells):
    write_crystal_structure(
        f"supercell-{i}.in",
        supercell,
        interface_mode=calc,
        optional_structure_info=optional_structure_info,
    )

phonon.save("phonopy_disp.yaml")

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:

phonon.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#

Result objects#

The Phonopy class holds the definition of the physical model: the three unit cells (Three unit cells), symmetry, force constants, and parameters of non-analytical term correction (Non-analytical term correction). Each phonon analysis is started by a run_* method, which computes the analysis and returns a self-contained result object. The same object is also accessible through the property of the same name (e.g. phonon.band_structure), which returns None before the corresponding run_* call. Mutating a calculation input (force constants, NAC parameters, masses, …) invalidates the stored results; run the analysis again.

Note

This part of the API is being restructured toward phonopy v5 / v6 (see Development). The run_* return value is the recommended access path for new code; the same-named properties on Phonopy are transitional and are planned for deprecation in a future major version. The data attributes of the result objects listed below are the stable surface.

bs = phonon.run_band_structure(paths)  # compute and return
bs.frequencies  # read data attributes
bs.write_yaml()  # file output
bs.plot(ax)  # draw into a matplotlib Axes

Internal machinery such as DynamicalMatrix and GroupVelocity is not a result object; those are stateful calculators wrapped by the run_* methods and are not part of the supported data-access API.

The result objects and their main data attributes:

`run_*` method	Result class	Main data attributes
`run_band_structure()`	`BandStructure`	`qpoints`, `distances`, `frequencies`, `eigenvectors`, `group_velocities`, `labels`, `path_connections`
`run_mesh()`	`Mesh`	`qpoints`, `weights`, `frequencies`, `eigenvectors`, `group_velocities`, `mesh_numbers`
`run_qpoints()`	`QpointsPhonon`	`qpoints`, `frequencies`, `eigenvectors`, `group_velocities`, `dynamical_matrices`
`run_total_dos()`	`TotalDos`	`frequency_points`, `dos`
`run_projected_dos()`	`ProjectedDos`	`frequency_points`, `projected_dos`
`run_thermal_properties()`	`ThermalProperties`	`temperatures`, `free_energy`, `entropy`, `heat_capacity`, `zero_point_energy`
`run_thermal_displacements()`	`ThermalDisplacements`	`temperatures`, `thermal_displacements`
`run_thermal_displacement_matrices()`	`ThermalDisplacementMatrices`	`temperatures`, `thermal_displacement_matrices`, `thermal_displacement_matrices_cif`
`run_moment()`	`PhononMoment`	`moment`
`run_modulations()`	`Modulation`	`modulations`, `supercell`, `modulated_supercells`, `modulated_supercell`, `frequencies`, `eigenvectors`
`run_irreps()`	`IrReps`	`band_indices`, `characters`, `frequencies`, `eigenvectors`, `qpoint`
`run_dynamic_structure_factor()`	`DynamicStructureFactor`	`qpoints`, `frequencies`, `dynamic_structure_factors`

File output lives on the result objects (write_yaml, write_hdf5, write, write_cif, depending on the class), and ax-level plotting on those that have a graphical representation (plot(ax)). Figure-level convenience plot functions are provided in phonopy.phonon.plot and as Phonopy.plot_* methods.

Save parameters (`phonopy.save`)#

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

phonon.save()

Force sets, displacements, Born effective charges, and dielectric constant are written in the default behaviour.

The default file name is phonopy_params.yaml, but this can be changed with the filename keword argument, which may be necessary if using certain CUI commands that expect a particular filename

The force constants can be written as follows:

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

Band structure#

Set band paths with run_band_structure(), which returns a BandStructure result object with attributes qpoints, distances, frequencies, eigenvectors, and group_velocities (the same object is also accessible via the band_structure property). Eigenvectors are included when with_eigenvectors=True is passed to run_band_structure(). Frequencies are returned in the unit set by Phonopy.unit_conversion_factor (THz by default for the VASP calculator with displacements in Angstrom and forces in eV/Angstrom). Imaginary frequencies are encoded as negative real numbers.

bs = phonon.run_band_structure(qpoints, path_connections=connections)
frequencies = bs.frequencies

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 optional unless nice looking plotting is needed. To obtain eigenvectors, the corresponding keyword argument must be set:

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#

Run the sampling-mesh phonon calculation with run_mesh() in reciprocal space. It returns a Mesh result object with attributes qpoints, weights, frequencies, eigenvectors, and group_velocities for the irreducible q-points (the same object is also accessible via the mesh property). mesh gives the sampling mesh in the Monkhorst-Pack scheme. The keyword shift gives the fractional mesh shift with respect to the neighboring grid points.

m = phonon.run_mesh([20, 20, 20])
qpoints = m.qpoints
weights = m.weights
frequencies = m.frequencies
eigenvectors = m.eigenvectors
group_velocities = m.group_velocities

To obtain eigenvectors, the corresponding keyword argument must be set:

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. The results are accessible via the total_dos and projected_dos properties.

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 from the returned ThermalProperties object with attributes temperatures, free_energy, entropy, and heat_capacity (the same object is also accessible via the thermal_properties property).

phonon.run_mesh([20, 20, 20])
tp = phonon.run_thermal_properties(t_step=10, t_max=1000, t_min=0)
temperatures = tp.temperatures
free_energy = tp.free_energy
entropy = tp.entropy
heat_capacity = tp.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.

The parameters are set to the nac_params attribute as a dictionary (typed as NacParams in phonopy.harmonic.dynamical_matrix) with the following keys:

Key	Type	Description
`'born'`	array_like, shape=(primitive cell atoms, 3, 3)	Born effective charge tensors in Cartesian coordinates, in the order of atoms of the primitive cell.
`'dielectric'`	array_like, shape=(3, 3)	High-frequency dielectric constant tensor in Cartesian coordinates.
`'factor'`	float, optional	Unit conversion factor of the non-analytical term. When omitted, the value for the calculator interface is used (see Default unit conversion factor for non-analytical term correction).
`'method'`	str, optional	NAC method, either `'gonze'` (Gonze-Lee, default) or `'wang'`.
`'G_cutoff'`	float, optional	Cutoff distance of reciprocal-space sampling of the Gonze-Lee method. When omitted, determined automatically.
`'Lambda'`	float, optional	Smearing parameter of the Ewald-like sum of the Gonze-Lee method. When omitted, determined automatically.

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}

Phonon at arbitrary q-points#

To evaluate phonons at an explicit list of q-points (without imposing a band path or a mesh), use run_qpoints(). It returns a QpointsPhonon result object with attributes qpoints, frequencies, eigenvectors, group_velocities, and dynamical_matrices (the same object is also accessible via the qpoints property).

qpts = phonon.run_qpoints(
    [[0.0, 0.0, 0.0], [0.5, 0.0, 0.0], [0.5, 0.5, 0.5]],
    with_eigenvectors=True,
    with_group_velocities=True,
)
frequencies = qpts.frequencies
eigenvectors = qpts.eigenvectors
group_velocities = qpts.group_velocities

run_qpoints() also covers one-off evaluations at a single q-point, e.g. phonon.run_qpoints([q]).frequencies[0]. The legacy single-point helpers (get_frequencies, get_frequencies_with_eigenvectors, get_group_velocity_at_q, get_dynamical_matrix_at_q) are deprecated.

Thermal displacements#

Mean-square atomic displacements and the corresponding 3x3 matrices (useful for cif export) are calculated from a converged mesh sampling.

phonon.run_mesh([20, 20, 20], with_eigenvectors=True, is_mesh_symmetry=False)
td = phonon.run_thermal_displacements(t_min=0, t_max=1000, t_step=10)
temperatures = td.temperatures
u2 = td.thermal_displacements  # shape=(temperatures, atoms*3)

tdm = phonon.run_thermal_displacement_matrices(t_min=0, t_max=1000, t_step=10)
temperatures = tdm.temperatures
phonon.write_thermal_displacement_matrix_to_cif(temperature_index=10)

Modulations#

run_modulations() builds atomic-displacement patterns of selected phonon modes on a chosen supercell. Each requested mode is a list [q-point, band_index, amplitude, phase]. The returned Modulation object gives the per-mode modulated cells (modulated_supercells), the cell modulated by all modes summed (modulated_supercell), the complex displacement fields (modulations), and the perfect supercell (supercell).

mod = phonon.run_modulations(
    dimension=[2, 2, 2],
    phonon_modes=[
        [[0.0, 0.0, 0.0], 0, 1.0, 0.0],
        [[0.5, 0.5, 0.5], 3, 0.5, 0.0],
    ],
)
modulated_cells = mod.modulated_supercells
phonon.write_modulations()  # MPOSCAR-001, MPOSCAR-002, ...

Irreducible representations#

For mode analysis at a chosen q-point, irreducible representations of the little co-group can be computed.

phonon.run_irreps(q=[0.0, 0.0, 0.0])
phonon.show_irreps(show_irreps=True)
phonon.write_yaml_irreps()

Dynamic structure factor#

See Dynamic structure factor for the formulation and the run_dynamic_structure_factor() / init_dynamic_structure_factor() workflow. A mesh calculation with eigenvectors and without mesh symmetry is required as a prerequisite.

Random displacements at finite temperature#

Phonopy can sample displacement snapshots from the canonical ensemble of harmonic oscillators (init_random_displacements() and get_random_displacements_at_temperature()). The same functionality is also exposed through generate_displacements(temperature=...). See Random displacements for details and recommended use.

Machine learning potentials#

Force constants can be obtained via a machine-learning potential trained on a type-2 displacement-force dataset. The relevant methods are develop_mlp(), save_mlp(), load_mlp(), and evaluate_mlp(), which work in concert with the mlp_dataset and mlp attributes. See Temperature dependent force constants using pypolymlp and symfc for a worked example.

Force-constants symmetrization#

symmetrize_force_constants() refines force constants produced from displacements and forces. Two schemes are available.

Traditional scheme (default). Translational and permutation symmetries are applied successively. The level keyword sets how many times this successive (translation -> permutation) application is repeated. Because the two symmetries are applied one after the other rather than simultaneously, the resulting force constants can break space-group symmetry slightly.

phonon.symmetrize_force_constants(level=1)

symfc projector. Set use_symfc_projector=True to apply the symfc projector instead. Space-group, translational, and permutation symmetries are imposed simultaneously in a single shot, so no trade-off is made between them. The level keyword is not used in this mode.

phonon.symmetrize_force_constants(use_symfc_projector=True)

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,)
species_table = (None,)
species_ids = (None,)

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

cell
positions or scaled_positions
symbols, or numbers, or (species_table and species_ids)

symbols, numbers, and species_table are mutually exclusive. The species_table / species_ids pair is the canonical internal representation and is required to construct cells with mixed-species sites such as those arising from the Virtual Crystal Approximation; see Mixed-species sites and the Virtual Crystal Approximation.

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.

For cells that contain mixed-species sites (e.g. VCA), cell.numbers raises RuntimeError because a mixture has no single atomic number; use cell.species_ids instead.

`species_ids`#

Per-atom indices into cell.species_table. Each id is an opaque non-negative integer that refers to a _Species entry; two atoms share an id iff they are the same chemical species (including any suffix index or mixture content). Available on every cell, including those without mixed sites.

`species_table`#

Deduplicated list of _Species entries indexed by species_ids. Each entry holds either a single-element species (atomic_number set) or a mixed-species site (mixture set to a tuple of (symbol, weight) pairs summing to 1.0). The list is a shallow copy on each access; entries are frozen and safe to share.

`has_mixtures`#

Boolean property. True when any species in the cell is a weighted mixture of constituents (e.g. a VCA virtual-crystal site).

`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
species_ids
species_table
has_mixtures
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).

Mixed-species sites and the Virtual Crystal Approximation#

Warning

Experimental. The support of mixed-species sites, including the Virtual Crystal Approximation (VCA) of VASP, is experimental. The site_mixture_scheme parameter and the unmerged cells of Phonopy, the --site-mixture and --split-site-mixture options of phonopy-init, the site_mixture_scheme entry of phonopy_disp.yaml, and the VASP for_vca writer may change without notice in upcoming releases. Currently only the VASP calculator interface writes supercells for VCA.

A mixed-species site is a crystallographic site that is shared by several elements in fixed proportions, for example a site of a Ge/Sn alloy that is occupied by Ge with the weight 0.9 and by Sn with the weight 0.1. Phonopy reads such a site in the same way as VASP does for VCA. The input structure has one atom per element at the same fractional coordinates, and each of these co-located atoms has a weight. The weights of the atoms of one site sum to 1.0, and an atom that is not co-located with another atom has the weight 1.0.

For a zincblende GeSn 50/50 cell, the input POSCAR has two Ge atoms and two Sn atoms, and the Ge atom and the Sn atom of each site are at the same fractional coordinates:

GeSn 50/50
0
000000 2.894478 2.894478
894478 0.000000 2.894478
894478 2.894478 0.000000
Ge Sn
2
Direct
00 0.00 0.00
25 0.25 0.25
00 0.00 0.00
25 0.25 0.25

Weighted atoms and `Phonopy`#

phonopy.structure.cells.apply_site_mixture attaches the weights to the atoms of a cell. The weights are given per atom in the order of the atoms in the cell. The atoms are kept as they are, so the returned cell has the same atoms in the same order as the input structure.

from phonopy import Phonopy
from phonopy.interface.calculator import read_crystal_structure
from phonopy.structure.cells import apply_site_mixture

cell, _ = read_crystal_structure("POSCAR", interface_mode="vasp")
weighted = apply_site_mixture(cell, [0.5, 0.5, 0.5, 0.5])
phonon = Phonopy(weighted, supercell_matrix=[2, 2, 2], primitive_matrix="P")

The site_mixture_scheme parameter of Phonopy decides how the co-located atoms are calculated. Its default is "merge".

With site_mixture_scheme="merge", the co-located atoms of each site are treated as one site. phonon.unitcell, phonon.primitive and phonon.supercell have one atom per site, and the force constants are those between the sites. The mass of a site is the weighted average of the masses of its atoms. The cells with the atoms of the input structure are phonon.unmerged_unitcell, phonon.unmerged_primitive and phonon.unmerged_supercell. A displacement moves all the atoms of a site by the same amount, so the GeSn example has one displacement.
With site_mixture_scheme="split", every atom of the input structure is displaced and calculated by itself, and the force constants are those between the atoms. The weights scale the masses in the dynamical matrix. phonon.unmerged_unitcell, phonon.unmerged_primitive and phonon.unmerged_supercell are None. The GeSn example has two displacements, one of a Ge atom and one of a Sn atom.

With both schemes, the supercells in phonon.supercells_with_displacements have the atoms of the input structure in the same order, and the forces given to phonon.forces are those on these atoms. With the merge scheme, the forces on the atoms of each site are summed to the force on the site before the force constants are calculated. For VASP, the forces are summed as they are, because the VCA calculation already includes the weights in the forces.

A cell whose sites are already merged, which build_mixture_cell returns, is not accepted by Phonopy.

CLI: `--site-mixture`#

At the command line, the weights are given to phonopy-init with --site-mixture, per atom in the order of the atoms in the input structure:

phonopy-init --dim 2 2 2 --pa P -d --site-mixture 0.5 0.5 0.5 0.5

This uses the merge scheme. Adding --split-site-mixture uses the split scheme. The scheme is written in phonopy_disp.yaml,

phonopy:
  version: "..."
  site_mixture_scheme: merge  # merge or split

and the cells in phonopy_disp.yaml have the atoms of the input structure. The phonopy command reads the scheme from this file, and it does not accept --site-mixture and --split-site-mixture.

VASP POSCAR and INCAR for VCA#

For VASP, the supercells are written in POSCAR files in the order of their atoms. Consecutive atoms of one element with one weight are written in one species row of the POSCAR, and each species row has one weight of the INCAR VCA tag. The command prints the species rows, the order in which the POTCAR files are concatenated, and the VCA line. For the GeSn example, both schemes write SPOSCAR with the species rows Ge Sn and the counts 16 16, and the command prints

VASP VCA hint:
  POSCAR species rows: Ge Sn
  POSCAR counts:       16 16
  POTCAR order:        Ge Sn  (concat one POTCAR per row; duplicate as needed)
  INCAR:               VCA = 0.5 0.5

The same weights are returned by phonopy.interface.vasp.get_vasp_vca_weights(cell) for a cell with weighted atoms.

Definitions of variables#

Supercell matrix#

Supercell matrix \(\mathrm{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} ) \mathrm{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 = ((original_lattice.T @ supercell_matrix).T,)

Primitive matrix#

Primitive matrix \(\mathrm{M}_\mathrm{p}\) is a transformation 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

where \(\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 = (
    (supercell_lattice.T @ np.linalg.inv(supercell_matrix) @ primitive_matrix).T,
)

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 that creates a Phonopy instance by loading data from a phonopy_xxx.yaml file, which may include all the necessary information to run phonopy. A typical usage is:

import phonopy

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

If phonopy_params.yaml contains a displacement-force dataset and you want to avoid producing force constants, set produce_fc=False:

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

Alternatively, if you have either phonopy.yaml (or phonopy_disp.yaml) along with a FORCE_SETS file, you can create a Phonopy object like this:

phonon = phonopy.load("phonopy.yaml", force_sets_filename="FORCE_SETS")

In this case, the command reads the structure information from phonopy.yaml and the displacement-force data from FORCE_SETS to create the Phonopy instance.

If your current directory contains the following files:

% ls
BORN  FORCE_SETS  phonopy.yaml

then both the BORN and FORCE_SETS files will be read automatically by

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

For more details, see the function’s docstring:

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

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.physical_units import get_physical_units
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)

units = get_physical_units()
nac_params = {'born': borns_,
              'factor': units.Hartree * units.Bohr,
              'dielectric': epsilon_}

Phonopy API for Python

Contents

Phonopy API for Python#

Three unit cells#

Import modules#

Work flow#

Pre-process#

Calculators#

Pre-processing Example#

Post process#

Phonon calculation#

Result objects#

Save parameters (phonopy.save)#

Band structure#

Mesh sampling#

DOS and PDOS#

Thermal properties#

Non-analytical term correction#

Phonon at arbitrary q-points#

Thermal displacements#

Modulations#

Irreducible representations#

Dynamic structure factor#

Random displacements at finite temperature#

Machine learning potentials#

Force-constants symmetrization#

Data structure#

Eigenvectors#

Band structure#

Mesh sampling#

PhonopyAtoms class#

Initialization#

Variables#

cell#

scaled_positions#

positions#

symbols#

numbers#

species_ids#

species_table#

has_mixtures#

masses#

Attributes#

Methods#

Mixed-species sites and the Virtual Crystal Approximation#

Weighted atoms and Phonopy#

CLI: --site-mixture#

VASP POSCAR and INCAR for VCA#

Definitions of variables#

Supercell matrix#

Primitive matrix#

Symmetry search tolerance#

Load phonopy settings phonopy.load#

Read and write crystal structures#

Getting parameters for non-analytical term correction#

Save parameters (`phonopy.save`)#

`PhonopyAtoms` class#

`cell`#

`scaled_positions`#

`positions`#

`symbols`#

`numbers`#

`species_ids`#

`species_table`#

`has_mixtures`#

`masses`#

Weighted atoms and `Phonopy`#

CLI: `--site-mixture`#

Load phonopy settings `phonopy.load`#