Ctrl+K

LAMMPS & phonopy calculation

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 LAMMPS input script for force calculation.

Supported LAMMPS structure input format#

Two important limitations are:

  • Crystal structure has to be described in the similar format for read_data.

  • Basis vectors are rotated to match the LAMMPS structure file format of triclinic simulation box.

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.

Masses has not been supported yet.

Example#

The crystal structure is converted to the LAMMPS structure format of triclinic simulation box.

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 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, mlp.lammp of gtinv-294 was obtained from Polynomial Machine Learning Potential Repository at Kyoto University.

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

    % 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:

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:

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.

Contents