"""POSCAR <-> Bilbao-style CIF (crystod-group --poscar2cif / --cif2poscar).
--poscar2cif converts a POSCAR into a CIF laid out like the files produced by
the Bilbao Crystallographic Server (http://www.cryst.ehu.es): conventional
cell with 4-decimal lattice parameters, the space-group number and
Hermann-Mauguin symbol, the full list of conventional-cell symmetry
operations as compact x,y,z strings, and one representative site per Wyckoff
orbit with 5-decimal fractional coordinates -- unlike the pymatgen CifWriter
layout (which quotes the operators, adds formula/volume/Z lines, and lists
sites with a multiplicity column). The output file is written next to the
input as <POSCAR>.cif.
--cif2poscar is the inverse: it reads any CIF (Bilbao or pymatgen flavour),
expands the symmetry operations, and writes the spglib-standardized
primitive cell as a POSCAR (the working format of the other crystod
commands; --conventional writes the conventional cell instead). The output
is the input path without the .cif suffix.
"""
from __future__ import annotations
import argparse
import os
from datetime import datetime
def build_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(
prog="crystod-group --poscar2cif",
description="Convert a POSCAR into a Bilbao-style CIF (<POSCAR>.cif).",
)
parser.add_argument(
"--cell",
"-c",
"--poscar",
dest="cell",
required=True,
metavar="POSCAR",
help="Input POSCAR file.",
)
parser.add_argument(
"--tolerance",
type=float,
default=0.01,
help="Symmetry-detection tolerance (symprec) in Angstrom (default: 0.01).",
)
parser.add_argument(
"--output",
default=None,
metavar="FILE",
help="Output CIF path (default: <POSCAR>.cif next to the input).",
)
return parser
def _operation_string(rotation, translation) -> str:
"""Compact Bilbao operator string 'x+1/2,-y,z' from an integer rotation
matrix and a fractional translation (both in the conventional basis)."""
from fractions import Fraction
labels = ("x", "y", "z")
parts = []
for i in range(3):
term = ""
for j in range(3):
entry = int(round(rotation[i][j]))
if entry == 1:
term += ("+" if term else "") + labels[j]
elif entry == -1:
term += "-" + labels[j]
elif entry != 0:
raise SystemExit(
"ERROR: non-crystallographic rotation entry in the "
"standardized setting; please report this case."
)
twelfths = translation[i] * 12
if abs(twelfths - round(twelfths)) > 1e-4:
raise SystemExit(
"ERROR: translation is not a multiple of 1/12 in the "
"standardized setting; please report this case."
)
fraction = Fraction(int(round(twelfths)), 12) % 1
if fraction:
term += f"+{fraction}"
parts.append(term)
return ",".join(parts)
def chemical_formula_parts(atomic_numbers, equivalent_atoms):
"""Reduced chemical formula as ordered (symbol, count) parts.
Ordering follows the conventional chemical notation:
- cations first, anions last (MgO, Li2O, MoS2);
- among the cations, the element on the most special Wyckoff site
first -- smallest orbit multiplicity, i.e. the letter closest to
``a`` (La3Ni2O7: La on 2b before Ni on 4e) -- then increasing
valence (SrTiO3, KNbO3, PbZrO3);
- anions by increasing valence (LaOF: O before F).
Valences come from the pymatgen oxidation-state guesser, with
electronegativity as the fallback; ``equivalent_atoms`` is the
spglib orbit assignment of the same cell as ``atomic_numbers``.
"""
import numpy as np
from pymatgen.core.periodic_table import Element
symbols = [Element.from_Z(int(z)).symbol for z in atomic_numbers]
counts: dict[str, int] = {}
for symbol in symbols:
counts[symbol] = counts.get(symbol, 0) + 1
# most special site of each element = smallest orbit multiplicity
min_multiplicity = {symbol: len(symbols) + 1 for symbol in counts}
equivalent = np.asarray(equivalent_atoms)
for orbit_id in set(int(v) for v in equivalent):
members = np.nonzero(equivalent == orbit_id)[0]
symbol = symbols[int(members[0])]
min_multiplicity[symbol] = min(min_multiplicity[symbol], len(members))
divisor = 0
for value in counts.values():
divisor = int(np.gcd(divisor, value))
reduced = {symbol: count // divisor for symbol, count in counts.items()}
def electronegativity(symbol: str) -> float:
try:
value = float(Element(symbol).X)
return value if value == value else 0.0 # NaN for noble gases
except Exception:
return 0.0
valence: dict[str, float] = {}
try:
from pymatgen.core import Composition
guesses = Composition(reduced).oxi_state_guesses(max_sites=-1)
if guesses:
valence = {symbol: float(v) for symbol, v in guesses[0].items()}
except Exception:
valence = {}
if valence:
anions = {symbol for symbol, v in valence.items() if v < 0}
elif len(reduced) >= 2:
anions = {max(reduced, key=electronegativity)}
else:
anions = set()
cations = sorted(
(symbol for symbol in reduced if symbol not in anions),
key=lambda symbol: (
min_multiplicity[symbol],
valence.get(symbol, electronegativity(symbol)),
electronegativity(symbol),
symbol,
),
)
anion_list = sorted(
anions,
key=lambda symbol: (
valence.get(symbol, -electronegativity(symbol)),
min_multiplicity[symbol],
symbol,
),
)
return [(symbol, reduced[symbol]) for symbol in cations + anion_list]
def format_chemical_formula(parts) -> str:
"""(('K', 2), ('Se', 1), ('O', 4)) -> 'K2SeO4'."""
return "".join(
f"{symbol}{count}" if count != 1 else symbol for symbol, count in parts
)
def format_chemical_formula_sum(parts) -> str:
"""(('K', 2), ('Se', 1), ('O', 4)) -> 'K2 Se O4' (IUCr
_chemical_formula_sum style)."""
return " ".join(
f"{symbol}{count}" if count != 1 else symbol for symbol, count in parts
)
def _wrap_coordinate(value: float) -> float:
wrapped = round(value % 1.0, 5) % 1.0
return abs(wrapped) # -0.0 -> 0.0
[docs]
def bilbao_cif_lines(structure, tolerance: float, title: str):
"""Bilbao-style CIF of a structure, as lines.
The conversion behind ``crystod-group --poscar2cif -c POSCAR``: the
structure is brought to the spglib-standardized (idealized)
conventional cell, which follows the International Tables (ITA)
setting and origin -- the convention of the Bilbao Crystallographic
Server -- and written with 4-decimal lattice parameters, the
space-group number and symbol, the full list of conventional-cell
operations as compact ``x,y,z`` strings, and one representative site
per Wyckoff orbit.
Args:
structure: A pymatgen ``Structure`` (e.g. ``Structure.from_file``).
tolerance: Symmetry-detection tolerance (symprec, Angstrom); the
command's default is ``0.01``. A distortion below it is
symmetrized away.
title: The ``data_`` block name (the command uses the POSCAR file
name).
Returns:
``(lines, info)``: the CIF lines without line terminators, and a
dict with ``number``, ``symbol``, ``n_operations`` and ``n_sites``.
Raises:
SystemExit: spglib could not determine the symmetry at this
tolerance (``ValueError`` when called through
``crystod.group``).
Example:
>>> from pymatgen.core import Structure
>>> from crystod import group
>>> from crystod.examples import example_path
>>> structure = Structure.from_file(example_path("221_PPOSCAR_ScF3"))
>>> lines, info = group.bilbao_cif_lines(structure, 0.01, "ScF3")
>>> info
{'number': 221, 'symbol': 'Pm-3m', 'n_operations': 48, 'n_sites': 2}
>>> lines[10]
'_symmetry_Int_Tables_number 221'
"""
import numpy as np
import spglib
from pymatgen.core import Lattice
from pymatgen.core.periodic_table import Element
cell = (
np.asarray(structure.lattice.matrix),
np.asarray(structure.frac_coords),
[site.specie.Z for site in structure],
)
standardized = spglib.standardize_cell(
cell, to_primitive=False, no_idealize=False, symprec=tolerance
)
if standardized is None:
raise SystemExit(
"ERROR: spglib could not determine the symmetry of the structure "
f"(tolerance {tolerance}); try another --tolerance."
)
lattice_matrix, positions, atomic_numbers = standardized
dataset = spglib.get_symmetry_dataset(standardized, symprec=1e-5)
def field(name): # spglib < 2.4 returns a dict, >= 2.4 an object
if isinstance(dataset, dict):
return dataset.get(name)
return getattr(dataset, name, None)
number = int(field("number"))
symbol = str(field("international")).replace("_", "")
rotations = np.asarray(field("rotations"))
translations = np.asarray(field("translations"))
# Bilbao presentation: proper operations first, then the improper block
order_key = sorted(
range(len(rotations)),
key=lambda i: 0 if np.linalg.det(rotations[i]) > 0 else 1,
)
operations = [
_operation_string(rotations[i], translations[i]) for i in order_key
]
equivalent = np.asarray(field("equivalent_atoms"))
orbit_representatives = []
seen = set()
for index, orbit_id in enumerate(equivalent):
if orbit_id not in seen:
seen.add(orbit_id)
orbit_representatives.append(index)
formula_sum = format_chemical_formula_sum(
chemical_formula_parts(atomic_numbers, equivalent)
)
now = datetime.now()
lattice = Lattice(lattice_matrix)
lines = [
"# Created by CrystOD (Bilbao Crystallographic Server style)",
"# https://www.cryst.ehu.es",
f"# Date: {now.strftime('%d/%m/%Y %H:%M:%S')}",
"",
f"# {title} -- non-magnetic",
"",
f"data_{title}",
f"{'_audit_creation_date':<35}{now.strftime('%Y-%m-%d')}",
f"{'_audit_creation_method':<35}\"CrystOD (Bilbao style)\"",
f"{'_chemical_formula_sum':<35}\"{formula_sum}\"",
f"{'_symmetry_Int_Tables_number':<35}{number}",
f"{'_symmetry_space_group_name_H-M':<35}\"{symbol}\"",
f"{'_cell_length_a':<35}{lattice.a:.4f}",
f"{'_cell_length_b':<35}{lattice.b:.4f}",
f"{'_cell_length_c':<35}{lattice.c:.4f}",
f"{'_cell_angle_alpha':<35}{lattice.alpha:.4f}",
f"{'_cell_angle_beta':<35}{lattice.beta:.4f}",
f"{'_cell_angle_gamma':<35}{lattice.gamma:.4f}",
"",
"loop_",
"_symmetry_equiv_pos_site_id",
"_symmetry_equiv_pos_as_xyz",
]
for index, operation in enumerate(operations, start=1):
lines.append(f"{index:>4} {operation}")
lines.extend(
[
"",
"loop_",
"_atom_site_label",
"_atom_site_type_symbol",
"_atom_site_fract_x",
"_atom_site_fract_y",
"_atom_site_fract_z",
"_atom_site_occupancy",
]
)
representatives = sorted(
orbit_representatives,
key=lambda index: (
Element.from_Z(int(atomic_numbers[index])).symbol,
tuple(np.round(positions[index], 5)),
),
)
counters: dict[str, int] = {}
for index in representatives:
element = Element.from_Z(int(atomic_numbers[index])).symbol
counters[element] = counters.get(element, 0) + 1
x, y, z = (_wrap_coordinate(value) for value in positions[index])
lines.append(
f"{element}{counters[element]} {element} "
f"{x:.5f} {y:.5f} {z:.5f} 1.0000"
)
lines.append("")
info = {
"number": number,
"symbol": symbol,
"n_operations": len(operations),
"n_sites": len(representatives),
}
return lines, info
def build_cif2poscar_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(
prog="crystod-group --cif2poscar",
description="Convert a CIF into a POSCAR (spglib-standardized cell).",
)
parser.add_argument(
"--cell",
"-c",
"--cif",
dest="cell",
required=True,
metavar="FILE.cif",
help="Input CIF file.",
)
parser.add_argument(
"--tolerance",
type=float,
default=0.01,
help="Symmetry-detection tolerance (symprec) in Angstrom (default: 0.01).",
)
parser.add_argument(
"--conventional",
action="store_true",
help="Write the conventional cell instead of the primitive cell.",
)
parser.add_argument(
"--output",
default=None,
metavar="FILE",
help="Output POSCAR path (default: the input path without .cif).",
)
return parser
[docs]
def poscar_lines(lattice_matrix, positions, atomic_numbers) -> list[str]:
"""POSCAR content of a cell, as lines.
The writer behind ``crystod-group --cif2poscar``: the crystod test-file
style (six decimals, ``direct`` coordinates, element tag on every
coordinate line, species grouped by first appearance).
Args:
lattice_matrix: Lattice vectors as rows, Angstrom, shape
``(3, 3)``.
positions: Fractional coordinates, shape ``(n_atoms, 3)``.
atomic_numbers: Atomic number of every atom.
Returns:
The POSCAR lines without line terminators, the last one empty (so
that joining with newlines ends the file with a newline).
Example:
>>> from crystod import group
>>> lines = group.poscar_lines([[4, 0, 0], [0, 4, 0], [0, 0, 4]],
... [[0, 0, 0], [0.5, 0.5, 0.5]], [55, 17])
>>> lines[:2] + lines[5:8]
['Cs1 Cl1', '1.0', 'Cs Cl', '1 1', 'direct']
>>> lines[8:]
['0.000000 0.000000 0.000000 Cs', '0.500000 0.500000 0.500000 Cl', '']
"""
from pymatgen.core.periodic_table import Element
symbols = [Element.from_Z(int(z)).symbol for z in atomic_numbers]
ordered_elements = []
for symbol in symbols:
if symbol not in ordered_elements:
ordered_elements.append(symbol)
counts = {element: symbols.count(element) for element in ordered_elements}
lines = [
" ".join(f"{element}{counts[element]}" for element in ordered_elements),
"1.0",
]
for row in lattice_matrix:
lines.append(" ".join(f"{value:.6f}" for value in row))
lines.append(" ".join(ordered_elements))
lines.append(" ".join(str(counts[element]) for element in ordered_elements))
lines.append("direct")
for element in ordered_elements:
for symbol, coords in zip(symbols, positions):
if symbol == element:
x, y, z = (_wrap_coordinate(value) for value in coords)
lines.append(f"{x:.6f} {y:.6f} {z:.6f} {element}")
lines.append("")
return lines
def cif2poscar_main(argv: list[str] | None = None) -> None:
args = build_cif2poscar_parser().parse_args(argv)
if not os.path.isfile(args.cell):
raise SystemExit(f"ERROR: CIF file not found: {args.cell}")
try:
import warnings
with warnings.catch_warnings():
warnings.simplefilter("ignore")
from pymatgen.core import Structure
structure = Structure.from_file(args.cell)
except ImportError as exc:
raise SystemExit(
"ERROR: --cif2poscar requires pymatgen (pip install pymatgen)."
) from exc
except Exception as exc:
raise SystemExit(f"ERROR: could not read {args.cell} as a CIF: {exc}")
import numpy as np
import spglib
cell = (
np.asarray(structure.lattice.matrix),
np.asarray(structure.frac_coords),
[site.specie.Z for site in structure],
)
standardized = spglib.standardize_cell(
cell,
to_primitive=not args.conventional,
no_idealize=False,
symprec=args.tolerance,
)
if standardized is None:
raise SystemExit(
"ERROR: spglib could not determine the symmetry of the structure "
f"(tolerance {args.tolerance}); try another --tolerance."
)
lattice_matrix, positions, atomic_numbers = standardized
dataset = spglib.get_symmetry_dataset(standardized, symprec=1e-5)
def field(name): # spglib < 2.4 returns a dict, >= 2.4 an object
if isinstance(dataset, dict):
return dataset.get(name)
return getattr(dataset, name, None)
if args.output:
output = args.output
elif args.cell.endswith(".cif"):
output = args.cell[: -len(".cif")]
else:
output = f"{args.cell}_POSCAR"
if os.path.exists(output):
print(f"NOTE: overwriting existing {output}")
with open(output, "w") as handle:
handle.write("\n".join(poscar_lines(lattice_matrix, positions, atomic_numbers)))
cell_kind = "conventional" if args.conventional else "primitive"
symbol = str(field("international")).replace("_", "")
print("\n* CIF -> POSCAR *")
print(f"input : {args.cell}")
print(f"output : {output}")
print(
f"space group: {symbol} (No. {int(field('number'))}), "
f"tolerance {args.tolerance}"
)
print(f"{cell_kind} cell, {len(atomic_numbers)} atoms\n")
def main(argv: list[str] | None = None) -> None:
args = build_parser().parse_args(argv)
if not os.path.isfile(args.cell):
raise SystemExit(f"ERROR: POSCAR file not found: {args.cell}")
try:
from pymatgen.core import Structure
except ImportError as exc:
raise SystemExit(
"ERROR: --poscar2cif requires pymatgen (pip install pymatgen)."
) from exc
try:
structure = Structure.from_file(args.cell)
except Exception as exc:
raise SystemExit(f"ERROR: could not read {args.cell} as a structure: {exc}")
title = os.path.basename(args.cell)
lines, info = bilbao_cif_lines(structure, args.tolerance, title)
output = args.output or f"{args.cell}.cif"
if os.path.exists(output):
print(f"NOTE: overwriting existing {output}")
with open(output, "w") as handle:
handle.write("\n".join(lines))
print("\n* POSCAR -> Bilbao-style CIF *")
print(f"input : {args.cell}")
print(f"output : {output}")
print(
f"space group: {info['symbol']} (No. {info['number']}), "
f"tolerance {args.tolerance}"
)
print(
f"{info['n_operations']} symmetry operations, "
f"{info['n_sites']} independent sites\n"
)
_warn_if_symmetrized(structure, args.tolerance, info)
def _warn_if_symmetrized(structure, tolerance: float, info: dict) -> None:
"""Warn when the CIF came out MORE symmetric than the input POSCAR.
The CIF carries the spglib-standardized structure, so a displacement
smaller than the tolerance is averaged away -- silently turning a
distorted structure into its own parent. A symmetry-mode analysis then
reports no distortion at all, with nothing in the output to explain why
(GaN F-43m -> I-42d, whose whole distortion is 0.0002 A against the
0.01 A default)."""
import numpy as np
import spglib
cell = (
np.asarray(structure.lattice.matrix),
np.asarray(structure.frac_coords),
[site.specie.Z for site in structure],
)
try:
tight = spglib.get_symmetry_dataset(cell, symprec=1e-6)
except Exception:
return
if tight is None:
return
number = tight["number"] if isinstance(tight, dict) else tight.number
if int(number) >= int(info["number"]):
return
symbol = (
tight["international"] if isinstance(tight, dict) else tight.international
)
print(
f"WARNING: at tolerance {tolerance} this structure is "
f"{info['symbol']} (No. {info['number']}), but its\n"
f" own coordinates carry only "
f"{str(symbol).replace('_', '')} (No. {int(number)}). The CIF holds "
"the SYMMETRIZED\n structure, so any distortion below "
f"{tolerance} A has been averaged away. Pass a\n"
" smaller --tolerance to keep it.\n"
)
if __name__ == "__main__":
main()
