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 keyward argument of interface_mode
can be
omitted. For QE,
unitcell, _ = 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,
[[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 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,
[[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.set_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 suuprcell 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})
Shortcut to load input files (phonopy.load
)¶
phonopy.load
is a convenient python method to create Phonopy
instance loading forces, displacements, and parameters for
non-analytical term correction. The details are found in the docstring
that can be seen by (e.g., in ipython)
In [1]: import phonopy
In [2]: help(phonopy.load)
Examples of how to use phonopy.load
are listed below.
phonopy_params.yaml
may contain all information needed to prepare
phonon calculation:
phonon = phonopy.load("phonopy_params.yaml")
More detailed configuration can be given as follows:
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()
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(unitcell_filename="POSCAR",
supercell_matrix=[2, 2, 2],
primitive_matrix='F')
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.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
calclation 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 contant volume at temperaturs are obtained by
get_thermal_properties_dict
, where the results are given as a
dictionary of temperaturs, 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:
symbols=None,
positions=None,
numbers=None,
masses=None,
scaled_positions=None,
cell=None
At least three arguments have to be given at the initialization, which are
cell
positions
orscaled_positions
symbols
ornumbers
Variables¶
The following variables are implemented in the PhonopyAtoms
class
in phonopy/structure/atoms.py
.
lattice_vectors
¶
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, lattice_vectors)
where np
means the numpy module (import numpy as np
).
Methods and attributes¶
set_cell(basis_vectors)
get_cell()
set_positions(positions)
get_positions()
set_scaled_positions(scaled_positions)
get_scaled_positions()
set_masses(masses)
get_masses()
set_chemical_symbols(symbols)
get_chemical_symbols()
get_magnetic_moments()
set_magnetic_moments(magmoms)
get_atomic_numbers()
set_atomic_numbers(numbers)
get_volume()
get_number_of_atoms()
The arguments have to be set in the structures shown in Variables. There are the respective attributes:
cell
positions
scaled_positions
masses
magnetic_moments
symbols
numbers
volume
unitcell.get_number_of_atoms()
is equivalent to
len(unitcell)
. The instance can be deep-copied by
unitcell.copy()
. print(unitcell)
shows human-readable crystal
structure in Yaml format. 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
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
(lattice_vectors). 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 tranformation matrix from lattice vectors to those of a super cell. Following a crystallography convention, the transformation is given by
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
(lattice_vectors). 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.
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_}
PhononDB at Kyoto university¶
The phonon calculation database at
http://phonondb.mtl.kyoto-u.ac.jp/ph20180417/index.html can be easily
used from phonopy-API. Downloading the raw data, e.g.,
mp-361-20180417.tar.lzma
and expand it. In the directory
mp-361-20180417
,
% ipython
or we can use jupyter notebook. The data is loaded by
In [1]: import phonopy
In [2]: ph = phonopy.load("phonon.yaml")
For example, the band structure is plotted by
In [3]: ph.auto_band_structure(plot=True).show()
and similarly for PDOS
In [4]: ph.auto_projected_dos(plot=True).show()