#!/usr/bin/env python3
# Pure-matplotlib POSCAR visualizer (no numpy) + bonds (minimum-image)
# Usage: python plot_poscar_with_bonds.py POSCAR [--point_size N] [--bond_length FLOAT] [--scale FLOAT]
import sys
import math
import argparse
import numpy as np
from mpl_toolkits.mplot3d import Axes3D # noqa: F401 (needed for 3D)
import matplotlib.pyplot as plt
from EPWpy.default.default_dicts import *
# --- small covalent radii table (Å). Expand if you need more elements.
COVALENT_RADII = {
'H': 0.31, 'He': 0.28,
'Li': 1.28, 'Be': 0.96, 'B': 0.84, 'C': 0.76, 'N': 0.71, 'O': 0.66, 'F': 0.57, 'Ne': 0.58,
'Na': 1.66, 'Mg': 1.41, 'Al': 1.21, 'Si': 1.11, 'P': 1.07, 'S': 1.05, 'Cl': 1.02, 'Ar': 1.06,
'K': 2.03, 'Ca': 1.76, 'Sc': 1.70, 'Ti': 1.60, 'V': 1.53, 'Cr': 1.39, 'Mn': 1.39, 'Fe': 1.32,
'Co': 1.26, 'Ni': 1.24, 'Cu': 1.32, 'Zn': 1.22, 'Ga': 1.22, 'Ge': 1.20, 'As': 1.19, 'Se': 1.20,
'Br': 1.20, 'Kr': 1.16, 'Rb': 2.20, 'Sr': 1.95, 'Y': 1.90, 'Zr': 1.75, 'Nb': 1.64, 'Mo': 1.54,
'Tc': 1.47, 'Ru': 1.46, 'Rh': 1.42, 'Pd': 1.39, 'Ag': 1.45, 'Cd': 1.44, 'In': 1.42, 'Sn': 1.39,
'Sb': 1.39, 'Te': 1.38, 'I': 1.39, 'Xe': 1.40, 'Cs': 2.44, 'Ba': 2.15, 'La': 2.07, 'Ce': 2.04,
'W': 1.62, 'Pt': 1.36, 'Au': 1.36, 'Pb': 1.46,
}
FALLBACK_RADIUS = 0.75
# ----------------- utility math (no numpy) -----------------
[docs]
def dot(a, b): return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
[docs]
def norm(v): return math.sqrt(dot(v, v))
[docs]
def add(a, b): return [a[0]+b[0], a[1]+b[1], a[2]+b[2]]
[docs]
def sub(a, b): return [a[0]-b[0], a[1]-b[1], a[2]-b[2]]
[docs]
def scale_vec(v, s): return [v[0]*s, v[1]*s, v[2]*s]
[docs]
def mat_det_3(m):
(a,b,c),(d,e,f),(g,h,i) = m
return a*(e*i - f*h) - b*(d*i - f*g) + c*(d*h - e*g)
[docs]
def mat_scale(m, s):
return [[m[0][0]*s, m[0][1]*s, m[0][2]*s],
[m[1][0]*s, m[1][1]*s, m[1][2]*s],
[m[2][0]*s, m[2][1]*s, m[2][2]*s]]
[docs]
def mat_vec(m, v):
return [
m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2],
m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2],
m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2],
]
[docs]
def vec_add(a, b):
return [a[0]+b[0], a[1]+b[1], a[2]+b[2]]
[docs]
def vec_sub(a, b):
return [a[0]-b[0], a[1]-b[1], a[2]-b[2]]
[docs]
def vec_min(a, b):
return [min(a[0], b[0]), min(a[1], b[1]), min(a[2], b[2])]
[docs]
def vec_max(a, b):
return [max(a[0], b[0]), max(a[1], b[1]), max(a[2], b[2])]
# ----------------- POSCAR parsing (your original) -----------------
[docs]
def read_nonempty_lines(path):
with open(path, "r", encoding="utf-8") as f:
lines = [ln.strip() for ln in f.readlines()]
return [ln for ln in lines if ln != ""]
[docs]
def is_all_int_tokens(tokens):
try:
for t in tokens:
int(t)
return True
except ValueError:
return False
[docs]
def strtovecs(lines3):
vecs = []
for ln in lines3:
parts = ln.split()
if len(parts) < 3:
raise ValueError("Lattice vector line must have 3 numbers.")
vecs.append([float(parts[0]), float(parts[1]), float(parts[2])])
return vecs
[docs]
def parse_poscar(path):
lines = read_nonempty_lines(path)
if len(lines) < 8:
raise ValueError("POSCAR too short or malformed.")
title = lines[0]
scale = float(lines[1])
lattice = strtovecs(lines[2:5]) # 3 lines
tokens5 = lines[5].split()
idx = 6
if is_all_int_tokens(tokens5): # VASP4-like
counts = [int(t) for t in tokens5]
elems = [f"E{i+1}" for i in range(len(counts))]
next_line = lines[idx].lower()
else: # VASP5-like
elems = tokens5
counts = [int(t) for t in lines[6].split()]
idx = 7
next_line = lines[idx].lower() if idx < len(lines) else ""
selective = False
if next_line.startswith("s"):
selective = True
idx += 1
next_line = lines[idx].lower() if idx < len(lines) else ""
if next_line.startswith("d") or next_line.startswith("f"):
coord_type = "Direct"
elif next_line.startswith("c") or next_line.startswith("k"):
coord_type = "Cartesian"
else:
coord_type = "Direct"
idx -= 1
idx += 1
n_atoms = sum(counts)
positions = []
tags = []
for e, c in zip(elems, counts):
tags.extend([e]*c)
for i in range(n_atoms):
if idx + i >= len(lines):
raise ValueError("Not enough position lines.")
parts = lines[idx + i].split()
if len(parts) < 3:
raise ValueError(f"Position line {i+1} malformed.")
positions.append([float(parts[0]), float(parts[1]), float(parts[2])])
idx += n_atoms
if scale > 0:
lattice = mat_scale(lattice, scale)
elif scale < 0:
target_vol = abs(scale)
vol_now = abs(mat_det_3(lattice))
if vol_now <= 0:
raise ValueError("Invalid lattice (zero or negative volume).")
factor = (target_vol / vol_now) ** (1.0/3.0)
lattice = mat_scale(lattice, factor)
return title, lattice, elems, counts, coord_type, positions, tags
# ----------------- coordinate helpers -----------------
[docs]
def frac_to_cart(frac, lattice):
a, b, c = lattice
return [
frac[0]*a[0] + frac[1]*b[0] + frac[2]*c[0],
frac[0]*a[1] + frac[1]*b[1] + frac[2]*c[1],
frac[0]*a[2] + frac[1]*b[2] + frac[2]*c[2],
]
[docs]
def cart_to_frac(cart, lattice):
# Build A with columns being lattice vectors: A = [[a0,b0,c0],[a1,b1,c1],[a2,b2,c2]]
A = [[lattice[0][0], lattice[1][0], lattice[2][0]],
[lattice[0][1], lattice[1][1], lattice[2][1]],
[lattice[0][2], lattice[1][2], lattice[2][2]]]
detA = mat_det_3(A)
if abs(detA) < 1e-12:
raise ValueError("Singular lattice")
(a,b,c),(d,e,f),(g,h,i) = A
adj = [
[ e*i - f*h, c*h - b*i, b*f - c*e],
[ f*g - d*i, a*i - c*g, c*d - a*f],
[ d*h - e*g, b*g - a*h, a*e - b*d],
]
num = [adj[0][0]*cart[0] + adj[0][1]*cart[1] + adj[0][2]*cart[2],
adj[1][0]*cart[0] + adj[1][1]*cart[1] + adj[1][2]*cart[2],
adj[2][0]*cart[0] + adj[2][1]*cart[1] + adj[2][2]*cart[2]]
return [num[0]/detA, num[1]/detA, num[2]/detA]
[docs]
def wrap_frac(fr):
return [fr[0] - math.floor(fr[0]), fr[1] - math.floor(fr[1]), fr[2] - math.floor(fr[2])]
# ----------------- bonding (minimum-image) -----------------
[docs]
def nearest_image_delta(frac_i, frac_j):
# fractional displacement from i -> nearest image of j (in [-0.5,0.5) range)
d = [frac_j[0]-frac_i[0], frac_j[1]-frac_i[1], frac_j[2]-frac_i[2]]
return [d[0] - round(d[0]), d[1] - round(d[1]), d[2] - round(d[2])]
[docs]
def integer_image_shift(frac_i, frac_j):
# integer image shift to apply to j so its image is nearest to i
raw = [frac_j[0]-frac_i[0], frac_j[1]-frac_i[1], frac_j[2]-frac_i[2]]
return [int(round(raw[0])), int(round(raw[1])), int(round(raw[2]))]
[docs]
def detect_bonds(fracs, tags, lattice, bond_length=None, scale=1.20, min_dist=0.10):
"""
fracs: list of fractional coords (not necessarily wrapped) or wrapped -- we'll accept wrapped or not
tags: element labels for each atom
lattice: 3x3 row-vector lattice
Returns list of (i,j,dx,dy,dz,dist)
"""
n = len(fracs)
bonds = []
# Ensure we have wrapped fracs for consistent comparisons:
fr_wrapped = [wrap_frac(f) for f in fracs]
for i in range(n):
for j in range(i+1, n):
# compute nearest-image fractional displacement and integer shift for j
dfrac = nearest_image_delta(fr_wrapped[i], fr_wrapped[j])
shift = integer_image_shift(fr_wrapped[i], fracs[j]) # use raw j to get integer shift (works)
# cart vector from i to j_image
dcart = frac_to_cart(dfrac, lattice)
dist = norm(dcart)
if dist <= min_dist:
continue
if bond_length is not None:
cutoff = bond_length
else:
el_i = tags[i]; el_j = tags[j]
ri = COVALENT_RADII.get(el_i, FALLBACK_RADIUS)
rj = COVALENT_RADII.get(el_j, FALLBACK_RADIUS)
cutoff = (ri + rj) * scale
if dist <= cutoff:
dx, dy, dz = shift[0], shift[1], shift[2]
bonds.append((i, j, dx, dy, dz, dist))
return bonds
# ----------------- rendering helpers (your original) -----------------
[docs]
def unit_cell_vertices(lattice):
o = [0.0, 0.0, 0.0]
a, b, c = lattice
ab = vec_add(a, b)
ac = vec_add(a, c)
bc = vec_add(b, c)
abc = vec_add(a, bc)
return [o, a, b, c, ab, ac, bc, abc]
[docs]
def draw_edge(ax, p0, p1, **kwargs):
ax.plot([p0[0], p1[0]], [p0[1], p1[1]], [p0[2], p1[2]], **kwargs)
[docs]
def draw_unit_cell(ax, lattice, lw=1.5):
o = [0.0, 0.0, 0.0]
a, b, c = lattice
ab = vec_add(a, b)
ac = vec_add(a, c)
bc = vec_add(b, c)
abc = vec_add(a, bc)
edges = [
(o, a), (o, b), (o, c),
(a, ab), (a, ac),
(b, ab), (b, bc),
(c, ac), (c, bc),
(ab, abc), (ac, abc), (bc, abc),
]
for p0, p1 in edges:
draw_edge(ax, p0, p1, linewidth=lw, color = 'k')
[docs]
def auto_equal_limits(ax, pts):
if not pts:
return
mn = pts[0][:]
mx = pts[0][:]
for p in pts[1:]:
mn = vec_min(mn, p)
mx = vec_max(mx, p)
for i in range(3):
if mx[i] == mn[i]:
mx[i] += 1e-6
cx = [(mn[i] + mx[i]) * 0.5 for i in range(3)]
span = max(mx[i] - mn[i] for i in range(3))
half = 0.55 * span
ax.set_xlim(cx[0] - half, cx[0] + half)
ax.set_ylim(cx[1] - half, cx[1] + half)
ax.set_zlim(cx[2] - half, cx[2] + half)
[docs]
def color_wheel(i):
palette = ['C0','C1','C2','C3','C4','C5','C6','C7','C8','C9']
return palette[i % len(palette)]
# ----------------- main plotting with bonds -----------------
[docs]
def plot_poscar(path, point_size=120, wrap=True, show_labels=True,
bond_length=None, scale=1.20, min_dist=0.10, highlight_si=True):
title, lattice, elems, counts, coord_type, positions, tags = parse_poscar(path)
# ensure fractional positions if needed (for wrapping), then cartesian for plotting
if coord_type.lower().startswith("d"):
fracs = [wrap_frac(p) if wrap else p[:] for p in positions]
carts = [frac_to_cart(fr, lattice) for fr in fracs]
else:
carts = [p[:] for p in positions] # already cartesian
if wrap:
# Convert to frac, wrap, back to cart using your Cramer's rule
A = [[lattice[0][0], lattice[1][0], lattice[2][0]],
[lattice[0][1], lattice[1][1], lattice[2][1]],
[lattice[0][2], lattice[1][2], lattice[2][2]]]
detA = mat_det_3(A)
if abs(detA) < 1e-12:
raise ValueError("Singular lattice for wrapping.")
def invA_mul(r):
(a,b,c),(d,e,f),(g,h,i) = A
adj = [
[ e*i - f*h, c*h - b*i, b*f - c*e],
[ f*g - d*i, a*i - c*g, c*d - a*f],
[ d*h - e*g, b*g - a*h, a*e - b*d],
]
num = [adj[0][0]*r[0] + adj[0][1]*r[1] + adj[0][2]*r[2],
adj[1][0]*r[0] + adj[1][1]*r[1] + adj[1][2]*r[2],
adj[2][0]*r[0] + adj[2][1]*r[1] + adj[2][2]*r[2]]
return [num[0]/detA, num[1]/detA, num[2]/detA]
fracs = [wrap_frac(invA_mul(r)) for r in carts]
carts = [frac_to_cart(fr, lattice) for fr in fracs]
# If fracs not defined (cartesian input with no wrap), compute fracs for bonding anyway
if coord_type.lower().startswith("d"):
fracs_for_bond = fracs
else:
# compute fractional coords for cart positions (not wrapped)
fracs_for_bond = [cart_to_frac(r, lattice) for r in carts]
# Detect bonds (minimum-image)
bonds = detect_bonds(fracs_for_bond, tags, lattice, bond_length=bond_length, scale=scale, min_dist=min_dist)
# Collect points for limits: cell vertices + atoms
verts = unit_cell_vertices(lattice)
all_pts = verts + carts
# Prepare plot
fig = plt.figure(figsize=(9, 8))
ax = fig.add_subplot(111, projection='3d')
ax.set_title(title)
# Draw cell
draw_unit_cell(ax, lattice, lw=1.5)
# Group atoms by element for legend/color
by_elem = {}
for pos, el in zip(carts, tags):
by_elem.setdefault(el, []).append(pos)
# Plot atoms
for i, el in enumerate(elems):
pts = by_elem.get(el, [])
if not pts:
continue
xs = [p[0] for p in pts]
ys = [p[1] for p in pts]
zs = [p[2] for p in pts]
ax.scatter(xs, ys, zs, s=point_size, depthshade=True, label=f"{el} ({len(pts)})", marker='o')
if show_labels and len(pts) <= 50:
for p in pts:
ax.text(p[0], p[1], p[2], el, fontsize=8, ha='center', va='center')
# Draw bonds: thin gray for all, red thick for Si-Si (if highlight_si)
for (i, j, dx, dy, dz, dist) in bonds:
p_i = carts[i]
# compute image of j in Cartesian: j + dx*a + dy*b + dz*c
a, b, c = lattice[0], lattice[1], lattice[2]
p_j = carts[j] #add(carts[j], [dx*a[0] + dy*b[0] + dz*c[0],
# dx*a[1] + dy*b[1] + dz*c[1],
# dx*a[2] + dy*b[2] + dz*c[2]])
# draw base bond in light gray
if (np.linalg.norm([dx,dy,dz]) == 0.0):
print('dist',dist,p_i,p_j)
ax.plot([p_i[0], p_j[0]], [p_i[1], p_j[1]], [p_i[2], p_j[2]],
linewidth=3.0, color='0.5', zorder=0)
# highlight Si-Si special
if highlight_si and tags[i].strip().lower() == 'si' and tags[j].strip().lower() == 'si':
ax.plot([p_i[0], p_j[0]], [p_i[1], p_j[1]], [p_i[2], p_j[2]],
linewidth=2.6, color='red', zorder=1)
# label midpoint with distance
xm = 0.5*(p_i[0] + p_j[0]); ym = 0.5*(p_i[1] + p_j[1]); zm = 0.5*(p_i[2] + p_j[2])
ax.text(xm, ym, zm, f"{dist:.3f}Å", color='red', fontsize=7, ha='center', va='center')
# write bonds to file
with open("bonds.txt", "w", encoding="utf-8") as f:
f.write("# i j el_i el_j dx dy dz distance_A\n")
for (i, j, dx, dy, dz, dist) in bonds:
f.write(f"{i+1} {j+1} {tags[i]} {tags[j]} {dx} {dy} {dz} {dist:.6f}\n")
print(f"Wrote {len(bonds)} bonds to bonds.txt")
auto_equal_limits(ax, all_pts)
ax.set_xlabel("x (Å)")
ax.set_ylabel("y (Å)")
ax.set_zlabel("z (Å)")
ax.legend(loc="upper left", bbox_to_anchor=(0.02, 0.98))
plt.tight_layout()
plt.show()
# ----------------- CLI wrapper -----------------
if __name__ == "__main__":
parser = argparse.ArgumentParser(description="Pure-matplotlib POSCAR visualizer with bonds (no numpy)")
parser.add_argument("poscar", help="POSCAR file")
parser.add_argument("--point_size", type=float, default=120.0, help="Atom point size")
parser.add_argument("--bond_length", type=float, default=3.5, help="Absolute bond cutoff (Å). If not set, covalent radii * --scale is used")
parser.add_argument("--scale", type=float, default=1.20, help="Multiplier for covalent radii sum (if --bond_length not set)")
parser.add_argument("--min_dist", type=float, default=0.10, help="Ignore tiny distances (Å)")
parser.add_argument("--no-highlight-si", dest="highlight_si", action="store_false", help="Do not highlight Si-Si bonds")
args = parser.parse_args()
plot_poscar(args.poscar, point_size=args.point_size, wrap=True, show_labels=True,
bond_length=args.bond_length, scale=args.scale, min_dist=args.min_dist, highlight_si=args.highlight_si)
