(phonopy_module)= # Phonopy API for Python ```{note} The Python API is being restructured toward phonopy v5 and v6; see {ref}`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. ``` (phonopy_three_unit_cells)= ## 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` ({ref}`variable_supercell_matrix`). Then `primitive` is created from `supercell` with `supercell_matrix` and `primitive_matrix` ({ref}`variable_primitive_matrix`). Having the `supercell_matrix` ($\mathrm{M}_\mathrm{s}$) and the `primitive matrix` ($\mathrm{M}_\mathrm{p}$), for example, ```{math} \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}, ``` respectively, the basis vectors of `supercell` are built from those of `unitcell` 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} ) \mathrm{M}_\mathrm{s} ``` and the basis vectors of `primitive` are built from those of `supercell` by ```{math} ( \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: ```python from phonopy import Phonopy ``` Crystal structure is defined by the `PhonopyAtoms` class. The `PhonopyAtoms` module is imported by: ```python 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`. ```python 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, ```python 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 {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. ```python 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 {ref}`interfaces-physical-units` 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 {ref}`phonopy_read_write_structure` 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 {ref}`frequency_conversion_factor_tag`, {ref}`physical_unit_conversion`, and {ref}`calculator_interfaces`. #### Pre-processing Example ```python 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: ```python [ [ [ 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: ```python 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: ```python 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: ```python phonon.dataset = displacement_dataset ``` From the set of displacements and forces, force constants internally with calculated supercell sets of forces by ```python 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 ```python 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_result_objects)= ### Result objects The `Phonopy` class holds the definition of the physical model: the three unit cells ({ref}`phonopy_three_unit_cells`), symmetry, force constants, and parameters of non-analytical term correction ({ref}`phonopy_nac_params`). 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 {ref}`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. ``` ```python 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. (phonopy_save_parameters)= ### Save parameters (`phonopy.save`) Basic information and parameters needed for phonon calculation are saved into a file by `phonopy.save`. ```python 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: ```python 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. ```python bs = phonon.run_band_structure(qpoints, path_connections=connections) frequencies = bs.frequencies ``` In `example/NaCl`, the phonopy is executed from python script, e.g., ```python 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: ```python phonon.run_band_structure(bands, with_eigenvectors=True) ``` To obtain group velocities: ```python phonon.run_band_structure(bands, with_group_velocities=True) ``` Automatic selection of band paths using [SeeK-path](https://seekpath.readthedocs.io/en/latest/) is invoked by ```python phonon.auto_band_structure() ``` and to plot ```python 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. ```python 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: ```python phonon.run_mesh([20, 20, 20], with_eigenvectors=True) ``` and for group velocities: ```python 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}`mesh_tag`, for example: ```python 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, ```python phonon.run_mesh([20, 20, 20]) phonon.run_total_dos() phonon.plot_total_dos().show() ``` and projected DOS ```python 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: ```python phonon.auto_total_dos(plot=True).show() ``` and ```python 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). ```python 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() ``` (phonopy_nac_params)= ### 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 {ref}`nac_default_value_interfaces`). | | `'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. | ```python 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). ```python 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. ```python 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) ``` ```python 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`). ```python 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. ```python phonon.run_irreps(q=[0.0, 0.0, 0.0]) phonon.show_irreps(show_irreps=True) phonon.write_yaml_irreps() ``` ### Dynamic structure factor See {ref}`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 {ref}`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 {ref}`mlp-sscha` 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. ```python 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. ```python 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 ```python 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 ```python 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: ```python 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 {ref}`mixed_species_sites`. (phonopy_Atoms_variables)= ### Variables The following variables are implemented in the `PhonopyAtoms` class in `phonopy/structure/atoms.py`. (phonopy_Atoms_cell)= #### `cell` Basis vectors are given in the matrix form in Cartesian coordinates. ```python [ [ 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. ```python [ [ x1_a, x1_b, x1_c ], [ x2_a, x2_b, x2_c ], [ x3_a, x3_b, x3_c ], ... ] ``` #### `positions` Cartesian positions of atoms. ```python positions = np.dot(scaled_positions, cell) ``` where `np` means the numpy module (`import numpy as np`). #### `symbols` Chemical symbols, e.g., ```python ['Zn', 'Zn', 'O', 'O'] ``` for the ZnO unit cell. #### `numbers` Atomic numbers, e.g., ```python [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., ```python [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)= ### Mixed-species sites and the Virtual Crystal Approximation ```{warning} **Experimental.** The mixed-species / Virtual Crystal Approximation (VCA) support — including `species_table` / `species_ids` / `has_mixtures` on `PhonopyAtoms`, `build_mixture_cell`, `build_species_table_from_mixtures`, the `--site-mixture` CLI option, the FC-time mixture-force reduction, and the VASP `expand_mixtures` writer — is experimental. APIs, file layouts (including the expanded `FORCE_SETS` format), and CLI flags may change without notice in upcoming releases. Currently only the VASP calculator interface is wired for mixture-expanded I/O. ``` Each `PhonopyAtoms` instance owns a deduplicated species table plus a per-atom index list (`species_ids`). A species can either be an ordinary chemical element or a weighted mixture of elements; the latter represents a single crystallographic site shared by several species in fixed proportions (e.g. a Virtual Crystal Approximation site for a Ge/Sn alloy). Mixed-site masses are the weight-averaged sum of constituent atomic masses. #### Building a cell with mixed sites `build_species_table_from_mixtures` packs a per-atom list of `(symbol, weight)` tuples into the canonical `(species_table, species_ids)` pair: ```python from phonopy.structure.atoms import ( PhonopyAtoms, build_species_table_from_mixtures, ) species_table, species_ids = build_species_table_from_mixtures( [ [("Si", 1.0)], # ordinary Si [("Ge", 0.5), ("Sn", 0.5)], # GeSn mixed-species site ] ) cell = PhonopyAtoms( cell=[[a, 0, 0], [0, a, 0], [0, 0, a]], scaled_positions=[[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]], species_table=species_table, species_ids=species_ids, ) assert cell.has_mixtures ``` Constituent weights of each entry must sum to 1.0; a single-component entry of weight 1.0 is canonicalized to a normal (single-element) species. Mixed species carry composite labels formed by concatenating constituent symbols in input order ("GeSn"); two distinct mixtures that would share the same composite label still get distinct ids because the underlying `(symbol, weight)` tuples differ. #### `build_mixture_cell` utility `phonopy.structure.cells.build_mixture_cell` collapses overlapping atoms in an ordinary cell into mixed-species sites. The Virtual Crystal Approximation (VCA) is the typical use case. ```python from phonopy.structure.cells import build_mixture_cell mixed_cell = build_mixture_cell(cell, weights=[0.5, 0.5, 0.5, 0.5]) ``` `weights` must have one entry per atom in the input order. Atoms whose fractional positions agree (modulo lattice translations) within `symprec` are merged into a single site whose weights must sum to 1.0; isolated atoms must carry weight 1.0. When several distinct mixtures within the same cell would collide on the same composite label, all colliding sites get 1-based suffixes (`"GeSn1"`, `"GeSn2"`, ...). The returned cell can be fed to `Phonopy(...)` as a unit cell. The per-atom convention here is distinct from VASP's INCAR `VCA` tag, which lists one weight per element row in POSCAR; see the CLI section below for how `phonopy` generates a VASP-compatible POSCAR/INCAR pair. #### CLI: `--site-mixture` The same merge can be requested at the command line: ```bash phonopy-init --dim "2 2 2" --site-mixture "0.5 0.5 0.5 0.5" -d ``` `--site-mixture` takes a space-separated list of per-atom weights in the input order, with the same validation rules as `build_mixture_cell`. The merge is applied immediately after the unit cell is read, so the rest of the workflow (displacement generation, supercell construction, force-constant building) sees only the merged cell. When the calculator is VASP, the supercell and displaced supercells are written with each mixed-species site expanded into one POSCAR row per constituent at the same fractional coordinates, and a hint is printed showing the matching POTCAR concatenation order and the INCAR `VCA = ...` line. For example, a 50/50 GeSn cell with `--dim "2 2 2"` produces `SPOSCAR` with `Ge Sn / 16 16` and prints `INCAR: VCA = 0.5 0.5`. #### Expansion ordering The mapping from expanded row indices back to phonopy sites is fixed by the same routine that writes the VASP POSCAR: for each entry in `cell.species_table`, the atoms belonging to that entry are emitted once per constituent (in `mixture` order) at the original site coordinates. The shared helper `phonopy.structure.cells.get_mixture_expansion(cell)` returns the ordered list of `(site_index, weight)` pairs that this row order corresponds to. ## Definitions of variables (variable_supercell_matrix)= ### Supercell matrix Supercell matrix {math}`\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 ```{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} ) \mathrm{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 ```python supercell_lattice = (original_lattice.T @ supercell_matrix).T, ``` (variable_primitive_matrix)= ### Primitive matrix Primitive matrix {math}`\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 ```{math} ( \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} ``` where {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 ```python 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. (phonopy_load)= ## 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: ```python 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`: ```python 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: ```python 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: ```bash % ls BORN FORCE_SETS phonopy.yaml ``` then both the `BORN` and `FORCE_SETS` files will be read automatically by ```python phonon = phonopy.load("phonopy.yaml") ``` For more details, see the function's docstring: ```python In [1]: import phonopy In [2]: help(phonopy.load) ``` (phonopy_read_write_structure)= ## 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. ```ipython 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](https://github.com/phonopy/phonopy/blob/develop/phonopy/interface/calculator.py#L123). ## 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. ```python 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_} ```