LAMMPS & phonopy calculation#
How to handle LAMMPS input and output files in phonopy#
Phonopy assumes the LAMMPS calculation is performed in
units metalandatom_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#
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". _ ___ _ __ __| | / _ \ '_ \ / _` | | __/ | | | (_| | \___|_| |_|\__,_|
Run LAMMPS
% lmp_serial -in in.polymlp
Suppose that the LAMMPS output file name is renamed to
lammps_forces_Ti.0after the LAMMPS calculation.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.