(lammps_interface)=

# LAMMPS & phonopy calculation

## How to handle LAMMPS input and output files in phonopy

- Phonopy assumes the LAMMPS calculation is performed in `units metal` and `atom_style atomic`.
- Format of version 15Sep2022 or later of LAMMPS is assumed.
- Force calculation has to performed by a specific setting as presented at
  {ref}`lammps_input_script_format`.

(lammps_structure_input_format)=

## Supported LAMMPS structure input format

Two important limitations are:

- Crystal structure has to be described in the similar format for
  [read_data](https://docs.lammps.org/read_data.html).
- Basis vectors are rotated to match the LAMMPS structure file format of
  [triclinic simulation box](https://docs.lammps.org/Howto_triclinic.html).

### Supported `read_data` keywords in the header

```
atoms
atom types
xlo xhi
ylo yhi
zlo zhi
xy xz yz
```

### Supported `read_data` keywords in the body

```
Atoms
Atom Type Labels
```

`Atom Type Labels` is the keyword new in version 15Sep2022 of LAMMPS. See
[Type labels](https://docs.lammps.org/Howto_type_labels.html).

`Masses` has not been supported yet.

### Example

The crystal structure is converted to the LAMMPS structure format of [triclinic
simulation box](https://docs.lammps.org/Howto_triclinic.html).

LAMMPS triclinic simulation box requires basis vectors to be described in the
following way

```
a = (a_x 0 0)
b = (b_x b_y 0)
c = (c_x c_y c_z)
```

An example of HCP crystal structure is given as follows:

```
#

2 atoms
1 atom types

0.0 2.923479689273095 xlo xhi   # xx
0.0 2.531807678358337 ylo yhi   # yy
0.0 4.624022835916574 zlo zhi   # zz

-1.461739844636547 0.000000000000000 0.000000000000000 xy xz yz

Atom Type Labels

1 Ti

Atoms

1 Ti 0.000000000000001 1.687871785572226 3.468017126937431
2 Ti 1.461739844636549 0.843935892786111 1.156005708979144
```

(lammps_input_script_format)=

## LAMMPS input script for force calculation

Phonopy can reads forces only from LAMMPS output file written in a specific text
format. An example of LAMMPS input script for phonopy force calculation is shown
below.

```
units metal

read_data supercell-001

pair_style  polymlp
pair_coeff * * mlp.lammps dummy

dump phonopy all custom 1 force.* id type x y z fx fy fz
dump_modify phonopy format line "%d %d %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f"
run 0
```

where `supercell-001` is that generated by phonopy. The `pair_style` of
`polymlp` is a LAMMPS module of the polynomial machine learning potentials
provided at <https://sekocha.github.io/lammps/index-e.html>. For the HCP Ti
calculation found in the [example
directory](https://github.com/phonopy/phonopy/tree/develop/example), `mlp.lammp`
of gtinv-294 was obtained from [Polynomial Machine Learning Potential Repository
at Kyoto
University](http://cms.mtl.kyoto-u.ac.jp/seko/mlp-repository/index.html).

## How to run Ti-lammps example

1. Read a lammps input structure file and create supercells with

   ```
   % phonopy --lammps -c lammps_structure_Ti -d --dim 4 4 3
           _
     _ __ | |__   ___  _ __   ___   _ __  _   _
    | '_ \| '_ \ / _ \| '_ \ / _ \ | '_ \| | | |
    | |_) | | | | (_) | | | | (_) || |_) | |_| |
    | .__/|_| |_|\___/|_| |_|\___(_) .__/ \__, |
    |_|                            |_|    |___/
                                         2.18.0

   Compiled with OpenMP support (max 10 threads).
   Python version 3.10.8
   Spglib version 2.0.2

   Calculator interface: lammps
   Crystal structure was read from "lammps_structure_Ti".
   Unit of length: angstrom
   Displacements creation mode
   Settings:
     Supercell: [4 4 3]
   Spacegroup: P6_3/mmc (194)
   Use -v option to watch primitive cell, unit cell, and supercell structures.

   "phonopy_disp.yaml" and supercells have been created.

   Summary of calculation was written in "phonopy_disp.yaml".
                    _
      ___ _ __   __| |
     / _ \ '_ \ / _` |
    |  __/ | | | (_| |
     \___|_| |_|\__,_|
   ```

2. Run LAMMPS

   ```bash
   % lmp_serial -in in.polymlp
   ```

   Suppose that the LAMMPS output file name is renamed to `lammps_forces_Ti.0`
   after the LAMMPS calculation.

3. Make `FORCE_SETS`

   ```
   % phonopy -f lammps_forces_Ti.0
           _
     _ __ | |__   ___  _ __   ___   _ __  _   _
    | '_ \| '_ \ / _ \| '_ \ / _ \ | '_ \| | | |
    | |_) | | | | (_) | | | | (_) || |_) | |_| |
    | .__/|_| |_|\___/|_| |_|\___(_) .__/ \__, |
    |_|                            |_|    |___/
                                         2.18.0

   Compiled with OpenMP support (max 10 threads).
   Python version 3.10.8
   Spglib version 2.0.2

   Calculator interface: lammps
   Displacements were read from "phonopy_disp.yaml".
   1. Drift force of "lammps_forces_Ti.0" to be subtracted
    -0.00000000  -0.00000000   0.00000000
   Forces parsed from LAMMPS output were rotated by F=R.F(lammps) with R:
     1.00000 0.00000 0.00000
     0.00000 0.00000 0.00000
     0.00000 1.00000 1.00000
   "FORCE_SETS" has been created.
                    _
      ___ _ __   __| |
     / _ \ '_ \ / _` |
    |  __/ | | | (_| |
     \___|_| |_|\__,_|
   ```

## Supercell generation from unit cell defined in yaml

LAMMPS input structure that phonopy can read must follow the convention of the
basis vector orientation in Cartesian coordinates. If it is unconformable for
the phonon calculation setting, we can generate supercells with displacements
from a structure format in yaml. A silicon primitive cell example is given as
follows:

```yaml
lattice:
- [0.000000000000000, 2.733099421887393, 2.733099421887393] # a
- [2.733099421887393, 0.000000000000000, 2.733099421887393] # b
- [2.733099421887393, 2.733099421887393, 0.000000000000000] # c
points:
- symbol: Si # 1
  coordinates: [0.875000000000000, 0.875000000000000, 0.875000000000000]
- symbol: Si # 2
  coordinates: [0.125000000000000, 0.125000000000000, 0.125000000000000]
```

With this saved in `phonopy_unitcell.yaml` file, we can generate 2x2x2 supercell
of this primitive cell orientation by using a python script:

```python
import phonopy
from phonopy.interface.phonopy_yaml import read_cell_yaml
from phonopy.interface.calculator import write_supercells_with_displacements

cell = read_cell_yaml("phonopy_unitcell.yaml")
ph = phonopy.load(unitcell=cell, primitive_matrix='auto', supercell_matrix=[2, 2, 2], calculator='lammps')
ph.generate_displacements()
ph.save("phonopy_disp.yaml")
write_supercells_with_displacements("lammps", ph.supercell, ph.supercells_with_displacements)
```

The primitive and supercell structures are stored in `phonopy_disp.yaml` in the
original orientation. But `supercell-001` follows the convention of the LAMMPS
input structure file format.

## Appendix: Structure optimization using LAMMPS

It is necessary to relax crystal structure before starting phonon calculation.
At least vanishing residual forces on atoms are expected. Using LAMMPS, crystal
structure can be optimized under different constraints. The following is the
simplest optimization where only internal atomic positions are relaxed.

```
units metal

read_data unitcell

pair_style  polymlp
pair_coeff * * mlp.lammps dummy

variable etol equal 0.0
variable ftol equal 1e-8
variable maxiter equal 1000
variable maxeval equal 100000

minimize ${etol} ${ftol} ${maxiter} ${maxeval}

write_data dump.unitcell
```

More instruction is found at
<https://gist.github.com/lan496/e9dff8449cd7489f6722b276282e66a0>.