.. _phonopy_module: Phonopy API for Python ================================= .. contents:: :depth: 2 :local: **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 :ref:`calculator_interfaces`. :: 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 :ref:`workflow`. Pre-process ^^^^^^^^^^^^ The first step is to create a ``Phonopy`` object with at least two arguments, a unit cell (``PhonopyAtoms`` object, see :ref:`phonopy_Atoms`) and a supercell matrix (3x3 array, see :ref:`variable_supercell_matrix`). In the following example, a :math:`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 :ref:`frequency_conversion_factor_tag`, :ref:`physical_unit_conversion`, and :ref:`calculator_interfaces`. 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 :ref:`file_force_constants`. Phonon calculation ------------------- .. _phonopy_save_parameters: 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}) .. _phonopy_load: 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 :ref:`mp_tag`, 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 .. _phonopy_Atoms: ``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`` or ``scaled_positions`` - ``symbols`` or ``numbers`` .. _phonopy_Atoms_variables: Variables ^^^^^^^^^^ The following variables are implemented in the ``PhonopyAtoms`` class in ``phonopy/structure/atoms.py``. .. _phonopy_Atoms_cell: ``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``). ``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. 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 :ref:`phonopy_Atoms_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 ------------------------- .. _variable_primitive_matrix: Primitive matrix ^^^^^^^^^^^^^^^^^ Primitive matrix :math:`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 .. math:: ( \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 :math:`\mathbf{a}_\mathrm{u}`, :math:`\mathbf{b}_\mathrm{u}`, and :math:`\mathbf{c}_\mathrm{u}` are the column vectors of the original lattice vectors, and :math:`\mathbf{a}_\mathrm{p}`, :math:`\mathbf{b}_\mathrm{p}`, and :math:`\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 (:ref:`phonopy_Atoms_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) .. _variable_supercell_matrix: Supercell matrix ^^^^^^^^^^^^^^^^^ Supercell matrix :math:`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 .. math:: ( \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 :math:`\mathbf{a}_\mathrm{u}`, :math:`\mathbf{b}_\mathrm{u}`, and :math:`\mathbf{c}_\mathrm{u}` are the column vectors of the original lattice vectors, and :math:`\mathbf{a}_\mathrm{s}`, :math:`\mathbf{b}_\mathrm{s}`, and :math:`\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 (:ref:`phonopy_Atoms_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. 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()