(phonopy_module)=
# Phonopy API for Python

## 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_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 (`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.,

```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

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.

```python
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, 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. To get the results
`get_total_dos_dict()` and `get_projected_dos_dict()` can be used.

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 by `get_thermal_properties_dict`,
where the results are given as a dictionary of temperatures, Helmholtz free
energy, entropy, and heat capacity with keys `temperatures`, `free_energy`,
`entropy`, and `heat_capacity`, respectively.

```python
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.

```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}
```

## 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,
```

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)=
#### `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.

#### `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
volume
```

where `volume` is the getter only.

### Methods

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

## Definitions of variables

(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 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{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_}
```