"""VASP POSCAR I/O — read into :class:`CrystalData` and write supercell geometry.
Layer: **VASP-specific** (POSCAR format; also the geometry source for the WIEN2k path).
"""
from __future__ import annotations
import logging
from typing import List
import numpy as np
from mag4.constants import CrystalData
from mag4.lattice import lattice_matrix
log = logging.getLogger(__name__)
[docs]
def read_poscar(poscar: str, species_index: int) -> CrystalData:
"""Parse a VASP POSCAR / CONTCAR file.
Parameters
----------
poscar : str
Path to the POSCAR.
species_index : int
1-based index of the magnetic species in the POSCAR header.
"""
with open(poscar) as f:
lines = f.readlines()
scale = float(lines[1])
vec_a = np.array([float(x) for x in lines[2].split()[:3]]) * scale
vec_b = np.array([float(x) for x in lines[3].split()[:3]]) * scale
vec_c = np.array([float(x) for x in lines[4].split()[:3]]) * scale
a = np.linalg.norm(vec_a)
b = np.linalg.norm(vec_b)
c = np.linalg.norm(vec_c)
alpha = np.degrees(np.arccos(np.clip(np.dot(vec_b, vec_c)/(b*c), -1, 1)))
beta = np.degrees(np.arccos(np.clip(np.dot(vec_a, vec_c)/(a*c), -1, 1)))
gamma = np.degrees(np.arccos(np.clip(np.dot(vec_a, vec_b)/(a*b), -1, 1)))
cell = [a, b, c, alpha, beta, gamma]
lat_mat = lattice_matrix(a, b, c, alpha, beta, gamma)
lat_vecs = np.array([vec_a, vec_b, vec_c])
elements = lines[5].split()
nb_species = [int(x) for x in lines[6].split()]
line_idx = 7
hdr = lines[line_idx].strip()
if hdr.lower().startswith("selective"):
line_idx += 1
hdr = lines[line_idx].strip()
is_cartesian = hdr.lower().startswith(("c", "k"))
coord_start = line_idx + 1
lat_mat_direct = np.column_stack([vec_a, vec_b, vec_c])
inv_lat_mat = np.linalg.inv(lat_mat_direct)
all_species: List[list] = []
total = 1
line_idx = coord_start
for gi, count in enumerate(nb_species):
group: list = []
for li in range(1, count + 1):
raw = np.array([float(x) for x in lines[line_idx].split()[:3]])
frac = (inv_lat_mat @ raw) if is_cartesian else raw
label = f"{elements[gi]}{li}"
group.append([label, *frac.tolist(), 0, total])
line_idx += 1
total += 1
all_species.append(group)
log.info("Coordinate type: %s", "Cartesian" if is_cartesian else "Direct")
return CrystalData(cell, all_species, all_species[species_index - 1],
nb_species, total, elements, lat_mat, lat_vecs)
[docs]
def write_poscar(crystal: CrystalData, name: str) -> str:
"""Write a single POSCAR for the supercell geometry. Returns the file path."""
lat = crystal.lattice_matrix
all_atoms = [atom for group in crystal.all_species for atom in group]
poscar_path = f"POSCAR_{name}"
with open(poscar_path, "w") as f:
f.write(f"Supercell {name}\n")
f.write("1.0\n")
for k in range(3):
v = lat[:, k]
f.write(f" {v[0]:20.10f} {v[1]:20.10f} {v[2]:20.10f}\n")
f.write(" " + " ".join(crystal.elements) + "\n")
f.write(" " + " ".join(str(n) for n in crystal.nb_species) + "\n")
f.write("Direct\n")
for atom in all_atoms:
f.write(f" {atom[1]:16.10f} {atom[2]:16.10f} {atom[3]:16.10f}\n")
log.info("POSCAR written: %s", poscar_path)
return poscar_path