"""
amorphgen.utils.radii
----------------------
Atomic radii data, bonding classification, minsep calculation,
and density estimation for random structure generation.
This module contains:
- Shannon ionic radii (CN=4 and CN=6)
- Metallic radii (Goldschmidt CN=12)
- Pauling electronegativities
- Element classification (nonmetal, metalloid)
- Elemental solid densities
- Bond type classification
- Minsep calculation from radii
- Cell volume / density estimation
"""
from __future__ import annotations
import itertools
import logging
from collections import Counter
from functools import lru_cache
import numpy as np
from ase.data import covalent_radii, atomic_numbers, atomic_masses
logger = logging.getLogger(__name__)
# ==============================================================================
# Radii data tables
# ==============================================================================
# Shannon effective ionic radii (A).
# Keyed as {symbol: {oxidation_state: {CN: radius}}}.
# CN=6 is the default; CN=4 used when target_cn is specified.
# Source: Shannon, R.D. Acta Cryst. A32, 751-767 (1976).
# The anion-formers also carry the cation states in which they centre an
# oxoanion (C in CO3 2-, N in NO3 -, S in SO4 2-, Cl in ClO4 -, I in IO3 -) or
# bond as the H+ of a hydroxide. At the oxoanion's own coordination these
# reproduce its bond: C4+ CN3 -0.08 + O2- 1.40 = 1.32 A against 1.29 A in a
# carbonate, S6+ CN4 0.12 + 1.40 = 1.52 against 1.47. The negative radii are
# Shannon's (a bond shorter than r(O2-) on his r(O2-) = 1.40 scale).
SHANNON_IONIC_RADII = {
# Cations
"Li": {1: {4: 0.59, 6: 0.76}},
"Na": {1: {4: 0.99, 6: 1.02}},
"K": {1: {6: 1.38}},
"Rb": {1: {6: 1.52}},
"Cs": {1: {6: 1.67}},
"Be": {2: {4: 0.27, 6: 0.45}},
"Mg": {2: {4: 0.57, 6: 0.72}},
"Ca": {2: {6: 1.00}},
"Sr": {2: {6: 1.18}},
"Ba": {2: {6: 1.35}},
"B": {3: {4: 0.11, 6: 0.27}},
"Al": {3: {4: 0.39, 6: 0.535}},
"Ga": {3: {4: 0.47, 6: 0.620}},
"In": {3: {4: 0.62, 6: 0.800}},
"Tl": {1: {6: 1.500}, 3: {6: 0.885}},
"Si": {4: {4: 0.26, 6: 0.400}},
"Ge": {4: {4: 0.39, 6: 0.530}},
"Sn": {2: {6: 0.930}, 4: {4: 0.55, 6: 0.690}},
"Pb": {2: {6: 1.190}, 4: {6: 0.775}},
"Ti": {3: {6: 0.670}, 4: {4: 0.42, 6: 0.605}},
"Zr": {4: {4: 0.59, 6: 0.720}},
"Hf": {4: {4: 0.58, 6: 0.710}},
"V": {3: {6: 0.640}, 4: {6: 0.580}, 5: {4: 0.355, 6: 0.540}},
"Nb": {4: {6: 0.680}, 5: {4: 0.48, 6: 0.640}},
"Ta": {5: {6: 0.640}},
"Cr": {3: {6: 0.615}, 4: {6: 0.550}, 6: {4: 0.26, 6: 0.440}},
"Mo": {4: {6: 0.650}, 6: {4: 0.41, 6: 0.590}},
"W": {4: {6: 0.660}, 6: {4: 0.42, 6: 0.600}},
"Ir": {3: {6: 0.680}, 4: {6: 0.625}, 5: {6: 0.570}},
"Ru": {3: {6: 0.680}, 4: {6: 0.620}, 5: {6: 0.565}},
"Rh": {3: {6: 0.665}, 4: {6: 0.600}},
"Pd": {2: {4: 0.64, 6: 0.860}, 4: {6: 0.615}},
"Ag": {1: {6: 1.150}, 2: {6: 0.940}, 3: {6: 0.750}},
"Re": {4: {6: 0.630}, 6: {6: 0.550}, 7: {6: 0.530}},
"Os": {4: {6: 0.630}, 6: {6: 0.545}, 7: {6: 0.525}, 8: {6: 0.390}},
"Pt": {2: {4: 0.60, 6: 0.800}, 4: {6: 0.625}},
"Au": {3: {6: 0.850}, 5: {6: 0.570}},
"Hg": {2: {4: 0.96, 6: 1.020}},
"Mn": {2: {4: 0.66, 6: 0.830}, 3: {6: 0.645}, 4: {4: 0.39, 6: 0.530}},
"Fe": {2: {4: 0.63, 6: 0.780}, 3: {4: 0.49, 6: 0.645}},
"Co": {2: {4: 0.58, 6: 0.745}, 3: {6: 0.610}},
"Ni": {2: {4: 0.55, 6: 0.690}},
"Cu": {1: {4: 0.60, 6: 0.770}, 2: {4: 0.57, 6: 0.730}},
"Zn": {2: {4: 0.60, 6: 0.740}},
"Cd": {2: {4: 0.78, 6: 0.950}},
"Y": {3: {6: 0.900}},
"La": {3: {6: 1.032}},
"Ce": {3: {6: 1.010}, 4: {6: 0.870}},
# Lanthanides (Ln3+ CN6; plus Eu2+/Yb2+ and Pr4+/Tb4+ where relevant)
"Pr": {3: {6: 0.990}, 4: {6: 0.850}},
"Nd": {3: {6: 0.983}},
"Pm": {3: {6: 0.970}},
"Sm": {3: {6: 0.958}},
"Eu": {2: {6: 1.170}, 3: {6: 0.947}},
"Gd": {3: {6: 0.938}},
"Tb": {3: {6: 0.923}, 4: {6: 0.760}},
"Dy": {3: {6: 0.912}},
"Ho": {3: {6: 0.901}},
"Er": {3: {6: 0.890}},
"Tm": {3: {6: 0.880}},
"Yb": {2: {6: 1.020}, 3: {6: 0.868}},
"Lu": {3: {6: 0.861}},
"Bi": {3: {6: 1.030}},
"Sc": {3: {6: 0.745}},
# Actinide dioxide formers (An4+, Shannon CN6/CN8). Without these, ThO2 /
# UO2 / PuO2 had no tabulated 4+ radius, so (a) the MO2 dioxide test could
# not fire and they fell to metal_oxide, and (b) the radius fell back to the
# oversized Cordero value (Th 2.06 A) — together under-predicting their
# density by ~60%. With these entries they route to fluorite_dioxide (large
# 4+ cation) and use the proper ionic radius (ThO2 rho -> ~10.4 vs 10.0).
"Th": {4: {6: 0.940, 8: 1.050}},
"U": {4: {6: 0.890, 8: 1.000}},
"Pu": {4: {6: 0.860, 8: 0.960}},
# Metalloid in high oxidation state oxides (e.g. As2O5).
# As(3+) is not tabulated by Shannon and Sb(III/V) Shannon radii
# are too small (lone-pair distortion), giving unphysical sphere
# volumes; for Sb compounds AmorphGen falls back to the Cordero
# covalent radius which matches experimental Sb2O3 density better.
"As": {5: {4: 0.335, 6: 0.460}},
# Anions (CN-independent)
"O": {-2: {6: 1.400}},
"S": {-2: {6: 1.840}, 4: {6: 0.37}, 6: {4: 0.12, 6: 0.29}},
"Se": {-2: {6: 1.980}, 4: {6: 0.50}, 6: {4: 0.28, 6: 0.42}},
"Te": {-2: {6: 2.210}, 4: {6: 0.970}, 6: {6: 0.560}},
"F": {-1: {6: 1.330}},
"Cl": {-1: {6: 1.810}, 5: {3: 0.12}, 7: {4: 0.08, 6: 0.27}},
"Br": {-1: {6: 1.960}, 5: {3: 0.31}, 7: {4: 0.25, 6: 0.39}},
"I": {-1: {6: 2.200}, 5: {3: 0.44, 6: 0.95}, 7: {4: 0.42, 6: 0.53}},
"N": {-3: {6: 1.460}, 3: {6: 0.16}, 5: {3: -0.104, 6: 0.13}},
"H": {-1: {6: 1.400}, 1: {1: -0.38, 2: -0.18}},
"C": {-4: {6: 1.400}, 4: {3: -0.08, 4: 0.15, 6: 0.16}}, # C4- approximate; Shannon does not list it
"P": {-3: {6: 2.120}, 5: {4: 0.17, 6: 0.380}},
}
# Metallic radii (A, CN=12 Goldschmidt radii).
# Source: Greenwood & Earnshaw, Chemistry of the Elements (1997).
METALLIC_RADII = {
"Li": 1.52, "Na": 1.86, "K": 2.27, "Rb": 2.48, "Cs": 2.65,
"Mg": 1.60, "Ca": 1.97, "Sr": 2.15, "Ba": 2.17,
"Al": 1.43, "Ga": 1.35, "In": 1.67,
"Tl": 1.70, "Sn": 1.58, "Pb": 1.75, "Ti": 1.47, "Zr": 1.60,
"Hf": 1.59, "V": 1.34, "Nb": 1.46, "Ta": 1.46, "Cr": 1.28,
"Mo": 1.39, "W": 1.39, "Mn": 1.27, "Fe": 1.26, "Co": 1.25,
"Ni": 1.24, "Cu": 1.28, "Zn": 1.37, "Cd": 1.52, "Y": 1.80,
"La": 1.87, "Ce": 1.83, "Sc": 1.64, "Bi": 1.56, "Ge": 1.37,
"Si": 1.17, "Ir": 1.36,
"Ru": 1.34, "Rh": 1.34, "Pd": 1.37, "Ag": 1.44, "Re": 1.37,
"Os": 1.35, "Pt": 1.39, "Au": 1.44, "Hg": 1.51,
# Lanthanides (CN=12)
"Pr": 1.82, "Nd": 1.82, "Pm": 1.81, "Sm": 1.80, "Eu": 2.04,
"Gd": 1.80, "Tb": 1.78, "Dy": 1.77, "Ho": 1.77, "Er": 1.76,
"Tm": 1.75, "Yb": 1.94, "Lu": 1.72,
# Metalloids (rhombohedral / metallic phase radii)
"As": 1.39, "Sb": 1.59,
}
# Pauling electronegativities for bonding classification.
# Values are the Allred (1961) revised scale, which retains Pauling's
# original methodology but uses updated thermochemical bond-dissociation
# data. These are the values that appear in modern chemistry textbooks
# and the CRC Handbook.
# Allred, A. L. "Electronegativity values from thermochemical data."
# J. Inorg. Nucl. Chem. 17, 215-221 (1961).
# Originally: Pauling, L. The Nature of the Chemical Bond, 3rd ed.,
# Cornell Univ. Press (1960).
PAULING_EN = {
"H": 2.20, "Li": 0.98, "Be": 1.57, "B": 2.04, "C": 2.55,
"N": 3.04, "O": 3.44, "F": 3.98, "Na": 0.93, "Mg": 1.31,
"Al": 1.61, "Si": 1.90, "P": 2.19, "S": 2.58, "Cl": 3.16,
"K": 0.82, "Ca": 1.00, "Sc": 1.36, "Ti": 1.54, "V": 1.63,
"Cr": 1.66, "Mn": 1.55, "Fe": 1.83, "Co": 1.88, "Ni": 1.91,
"Cu": 1.90, "Zn": 1.65, "Ga": 1.81, "Ge": 2.01, "As": 2.18,
"Se": 2.55, "Br": 2.96, "Rb": 0.82, "Sr": 0.95, "Y": 1.22,
"Zr": 1.33, "Nb": 1.60, "Mo": 2.16, "Cd": 1.69, "In": 1.78,
"Sn": 1.96, "Sb": 2.05, "Te": 2.10, "I": 2.66, "Cs": 0.79,
"Ba": 0.89, "La": 1.10, "Ce": 1.12, "Hf": 1.30, "Ta": 1.50,
"W": 2.36, "Pb": 2.33, "Bi": 2.02, "Tl": 1.62, "Ir": 2.20,
"Ru": 2.20, "Rh": 2.28, "Pd": 2.20, "Ag": 1.93, "Re": 1.90,
"Os": 2.20, "Pt": 2.28, "Au": 2.54, "Hg": 2.00,
"Pr": 1.13, "Nd": 1.14, "Pm": 1.13, "Sm": 1.17, "Eu": 1.20,
"Gd": 1.20, "Tb": 1.10, "Dy": 1.22, "Ho": 1.23, "Er": 1.24,
"Tm": 1.25, "Yb": 1.10, "Lu": 1.27,
}
# Conventional anion charges, used to infer the cation oxidation state
# from a compound's stoichiometry by charge balance (e.g. Ti=+4 in
# TiO_2, Sb=+5 in Sb_2O_5). Limited to clear-cut anion-formers; complex
# / mixed-anion / non-stoichiometric systems fall through to the
# default "highest positive Shannon state" fallback.
ANION_CHARGES = {
"O": -2,
"S": -2,
"Se": -2,
"Te": -2,
"F": -1,
"Cl": -1,
"Br": -1,
"I": -1,
"N": -3,
"H": -1, # treated as hydride; only valid for binary hydrides
}
# Cations that, although Shannon tabulates several oxidation states, take one
# strongly dominant state in oxide / mixed-anion solids. Used ONLY to break a
# tie when more than one charge-balanced assignment exists (never to force an
# unbalanced one). Kept deliberately conservative: d0 / main-group cases with
# little real ambiguity in compounds. Genuinely variable cations (Fe, Mn, Co,
# Cu, Ce, Ag, ...) are omitted, so a compound with two such cations stays
# ambiguous and resolves to None (honest) rather than guessing.
_DOMINANT_OS = {
"Ti": 4, "Zr": 4, "Hf": 4, "Nb": 5, "Ta": 5, "V": 5,
"W": 6, "Mo": 6, "Cr": 3, "Sn": 4, "Pb": 2,
"Al": 3, "Ga": 3, "In": 3, "Sc": 3, "Y": 3,
}
@lru_cache(maxsize=512)
def _solve_oxidation_states(comp_items: tuple) -> dict | None:
"""Jointly assign integer oxidation states to every cation by charge
balance. Returns a {cation: state} dict, or None if the assignment is not
uniquely determined.
Works for arbitrary cation/anion mixtures: the anion charge is summed over
*all* the anions of :func:`anion_elements` (so oxynitrides, oxyfluorides,
oxysulfides balance the same way as simple oxides, and the S of a sulfate
or the N of a nitrate is solved as a cation). Cation states are enumerated
from the Shannon table; at most one cation may be absent from Shannon (it
becomes the single free variable solved by the remaining balance). When
several assignments balance, the one with the fewest :data:`_DOMINANT_OS`
violations wins; if that is still tied the system is genuinely ambiguous
and None is returned.
"""
composition = dict(comp_items)
anions = anion_elements(composition)
total_neg = sum(n * ANION_CHARGES[s] for s, n in composition.items()
if s in anions)
if total_neg == 0:
return None
cation_charge = -total_neg # total positive charge the cations must supply
cations = {s: n for s, n in composition.items() if s not in anions}
if not cations:
return None
known, unknown = {}, []
for s in cations:
states = sorted(k for k in SHANNON_IONIC_RADII.get(s, {}) if k > 0)
if states:
known[s] = states
else:
unknown.append(s)
if len(unknown) > 1:
return None # two un-tabulated cations -> underdetermined
known_syms = list(known)
sizes = 1
for s in known_syms:
sizes *= len(known[s])
if sizes > 20000:
return None # guard against pathological enumeration
solutions = []
for combo in itertools.product(*(known[s] for s in known_syms)):
assign = dict(zip(known_syms, combo))
known_pos = sum(composition[s] * o for s, o in assign.items())
if unknown:
(u,) = unknown
rem = cation_charge - known_pos
if rem <= 0 or rem % composition[u] != 0:
continue
assign[u] = rem // composition[u]
elif known_pos != cation_charge:
continue
solutions.append(assign)
if not solutions:
return None
if len(solutions) == 1:
return solutions[0]
def penalty(sol):
viol = sum(1 for s, o in sol.items()
if s in _DOMINANT_OS and o != _DOMINANT_OS[s])
# an oxoanion centre takes its top state where it can: Ag2SO4 is
# Ag+ and a sulfate, not Ag2+ and a sulfite
centre = sum(o for s, o in sol.items() if s in NONMETALS)
return (viol, -centre, sum(sol.values()))
solutions.sort(key=penalty)
if penalty(solutions[0]) == penalty(solutions[1]):
return None # still ambiguous after the dominant-state tie-break
return solutions[0]
[docs]
def infer_oxidation_state(sym: str, composition: dict) -> int | None:
"""Infer the oxidation state of *sym* in *composition* by charge balance.
Sums the charge of the anions (:func:`anion_elements`, charges from
:data:`ANION_CHARGES`) and assigns the remaining positive charge
equally across cations of the queried element when only one cation
type is present, or for any element that appears as the only
non-anion type. A nonmetal that acts as a cation here
(:func:`cation_nonmetals`: P in Li3PO4, S in Li2SO4) is solved like
any other; the Li of Li2SO4 is +1, not the +5 that counting S as an
anion gave. Returns ``None`` for an anion, for Te (a metalloid on the
anion table, left to the highest-state default whether it is the
anion or the cation), for the covalent hydrogenated networks (a-Si:H,
where charge balance against H- made Si50H50 Si+1), and when the
inference is ambiguous (multiple cation species, non-integer oxidation
state, or no anions at all), letting the caller fall back to a default
Shannon entry.
"""
if ((sym in ANION_CHARGES or sym in NONMETALS)
and sym not in cation_nonmetals(composition)):
return None # not a cation we resolve here
if _hydrogenated_host(composition) is not None:
return None # covalent network, no ionic charges
anions = anion_elements(composition)
total_neg = sum(n * ANION_CHARGES[s] for s, n in composition.items()
if s in anions)
if total_neg == 0:
return None # no anions -> can't balance
cations = {s: n for s, n in composition.items() if s not in anions}
if sym not in cations:
return None
# Single-cation case: balance directly.
if len(cations) == 1:
n_cat = cations[sym]
ox = -total_neg / n_cat
if abs(ox - round(ox)) < 1e-6 and ox > 0:
return int(round(ox))
return None
# Multi-cation case: solve the whole system jointly. This resolves
# multivalent cations whose state is fixed once the other cations are
# pinned (FeTiO3 -> Fe2+/Ti4+, LiNbO3 -> Li+/Nb5+, SrTiO3 -> Sr2+/Ti4+),
# and handles mixed anions (oxynitrides, oxyfluorides) since the balance
# sums over every anion-former. Returns None when genuinely ambiguous.
solution = _solve_oxidation_states(tuple(sorted(composition.items())))
if solution is None:
return None
return solution.get(sym)
# ==============================================================================
# Element classification
# ==============================================================================
#
# References for the three element-type sets used by ``classify_bond``:
#
# IUPAC, "Nomenclature of Inorganic Chemistry — IUPAC Recommendations
# 2005" (Red Book), RSC Publishing, ISBN 0-85404-438-8. Section
# IR-3 — defines the standard nonmetal / metalloid / metal columns.
#
# AmorphGen follows IUPAC's classification exactly. The METALS set is the
# implicit complement of NONMETALS ∪ METALLOIDS — no separate list is
# maintained, so any element outside those two sets (transition metals,
# lanthanides, actinides, ...) is treated as metallic.
#
# Notes on borderline elements (intentionally NOT in METALLOIDS):
#
# Sn — IUPAC post-transition metal. β-Sn (the room-T phase) is
# metallic. Amorphous Sn is metallic with liquid-like CN ≈ 6–8
# and superconducting Tc ≈ 5–7 K, higher than crystalline β-Sn
# (Bückel & Hilsch, Z. Phys. 138, 109 (1954); Bergmann, Phys.
# Rep. 27, 159 (1976)). α-Sn (diamond cubic, semimetal) only
# exists below 13 °C and is not a glass-former, so the metallic
# classification gives the correct behaviour for every case
# AmorphGen would generate (SnO2 rutile CN=6, β-Sn alloys,
# a-Sn).
# Se — IUPAC nonmetal. Although a-Se forms covalent chalcogenide
# chains, the Pauling-Δχ override in classify_bond already
# routes Ge–Se, As–Se etc. to "covalent" without requiring
# metalloid status; the chalcogenide compound-class rule
# likewise fires on the presence of S/Se/Te without checking
# metalloid membership. Promoting Se to metalloid would be
# redundant.
# Standard nonmetals + noble gases. IUPAC Red Book (2005), Section IR-3.
# Se is included here — although it sits on the chalcogen/metalloid border
# and some textbooks promote it to metalloid, IUPAC's recommendation is
# "nonmetal". Se-Se covalency in trigonal Se and a-Se is handled by the
# nonmetal-nonmetal rule in classify_bond (returns "covalent").
NONMETALS = frozenset({
"H", "He", "C", "N", "O", "F", "Ne", "P", "S", "Cl", "Ar",
"Se", "Br", "Kr", "I", "Xe", "At", "Rn",
})
# IUPAC's six recognised metalloids — B, Si, Ge, As, Sb, Te.
# (IUPAC Red Book 2005, Section IR-3.)
METALLOIDS = frozenset({"B", "Si", "Ge", "As", "Sb", "Te"})
# Typical coordination numbers for coordination-aware placement.
# Metalloids are always CN=4. For metals, CN depends on the anion context:
# oxides → CN=4-6, nitrides → CN=4, halides/sulfides → CN=6.
_ALWAYS_TETRAHEDRAL = frozenset({
"Si", "Ge", "B", "P", "As", # metalloids — always CN=4
"Be", # tetrahedral in BeF2, BeO (like Si in SiO2)
# Group-12 post-transition metals form zincblende-structured II-VI
# semiconductors with chalcogens (ZnS, ZnSe, ZnTe, CdS, CdSe, CdTe,
# HgS, HgSe, HgTe — all CN=4 in their stable phases). Without this,
# the chalcogenide branch of auto_target_cn defaults them to octahedral
# CN=6, which is wrong for these specific compounds.
"Zn", "Cd", "Hg",
})
# d-block transition metals that form *interstitial* carbides (rocksalt or
# close-related dense structures: TiC, ZrC, HfC, VC, NbC, TaC, CrC, MoC,
# WC, ...). These compounds want metallic + Cordero radii and a higher
# packing factor (≈ 0.60) than covalent carbides (SiC, B4C, P4C3, ...).
# Pauling Δχ does NOT cleanly separate these from covalent carbides
# because W–C (Δχ = 0.19) and Mo–C (Δχ = 0.39) are interstitial but have
# small Δχ — the discriminator has to be cation group, not electronegativity.
# Group-11 (Cu, Ag, Au) and Group-12 (Zn, Cd, Hg) do not form stable
# carbides and are excluded.
_TRANSITION_METAL_CARBIDE_CATIONS = frozenset({
"Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", # 3d
"Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", # 4d
"Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", # 5d
})
# Rutile-type dioxides (MO2 with a small, octahedrally-coordinated metal cation:
# TiO2, VO2, CrO2, MnO2, NbO2, SnO2, RuO2, IrO2, OsO2, PtO2, ...) pack far more
# densely than light oxides and get a higher packing factor (see
# PACKING_FACTORS["rutile_dioxide"]). They are identified *geometrically* by the
# radius-ratio rule rather than an element list: a 4+ cation radius below
# ~0.70 A (r/r_O < ~0.5) gives octahedral rutile/CaCl2-type packing, while larger
# cations form 8-coordinate fluorite/baddeleyite oxides (ZrO2 0.72, HfO2 0.71,
# CeO2 0.87) that stay `metal_oxide`. Metalloid dioxides (SiO2, GeO2) are
# excluded as covalent network formers; see _classify_compound.
# Elements that form diatomic / small-molecule gases rather than covalent
# solid networks — excluded from the pure-element semiconductor rule.
_DIATOMIC_NONMETALS = frozenset({"H", "N", "O", "F", "Cl", "Br", "I"})
# Halogen anion-formers, shared by the halide / oxyhalide classification and
# the oxyhalide packing-factor interpolation.
_HALOGENS = frozenset({"F", "Cl", "Br", "I"})
# Minimum halogen fraction of the (halogen + O) anion pool for a compound to
# classify as oxyhalide. Below this the halogen is a dopant, not a framework
# anion (e.g. F-doped TiO2 / SnO2:F), and the compound keeps its oxide
# routing — otherwise one F atom would flip rutile_dioxide (pf 0.66) to
# oxyhalide (pf ~0.52), inflating the estimated cell by ~25%. Real oxyhalides
# sit far above this: Sb4O5Cl2 ~0.29, BiOCl 0.50, NaTaOCl4 0.80.
_OXYHALIDE_MIN_HALOGEN_FRAC = 0.10
_RUTILE_DIOXIDE_RMAX = 0.70 # A; Shannon 4+ CN6 cation-radius cutoff
# Cation-radius cutoff (A) separating small-cation nitrides (denser packing,
# placeable at higher pf) from large-cation nitrides (placement-limited by the
# cation-cation distance, kept at the safe default pf).
_SMALL_CATION_NITRIDE_RMAX = 0.70
# Covalent rutile formers whose 4+ radius sits above the cutoff (lone-pair /
# covalency), kept as explicit exceptions.
_RUTILE_DIOXIDE_EXCEPTIONS = frozenset({"Pb"})
# Cations whose oxide is a covalent tetrahedral network rather than an ionic
# solid. Beryllium is the textbook case: Be2+ is so small and polarizing (and
# not even tabulated by Shannon) that BeO is a covalent wurtzite network. The
# ionic sphere model under-predicts its density by ~50% (1.42 vs the measured
# amorphous 3.01 g/cm3); treating it with Cordero covalent radii for BOTH atoms
# + an open-network packing factor recovers ~2.96. Kept as an explicit element
# set (not a radius gate) so it cannot misfire on small-but-ionic cations such
# as Al3+. Distinct from "network_forming_oxide" (SiO2/GeO2/B2O3), whose metalloid
# cations have usable small Shannon radii and stay on the ionic path.
_COVALENT_OXIDE_CATIONS = frozenset({"Be"})
# Glass-forming cations of the corner-sharing chalcogenide glasses
_CHALCOGENIDE_GLASS_FORMERS = frozenset({"Ge", "Si", "As", "Sb", "B", "P"})
# Metalloids that form metal-rich metallic glasses with transition / noble
# metals (Ni80P20, Fe80B20, Pd80Si20); below this fraction the composition is
# an alloy, not a pnictide or boride
_GLASS_METALLOIDS = frozenset({"P", "B", "Si", "Ge", "As", "Sb"})
# Alkali and alkaline-earth metals are excluded from the metal-rich glass rule:
# Li3P, Na3Sb, Mg2Si and the like are Zintl phases / semiconductors, not
# metallic glasses, and keep their pnictide / boride treatment
_S_BLOCK_METALS = frozenset({"Li", "Na", "K", "Rb", "Cs", "Fr",
"Be", "Mg", "Ca", "Sr", "Ba", "Ra"})
_METAL_RICH_ALLOY_MAX_METALLOID_FRAC = 0.35
# Hosts of the hydrogenated group-IV networks (a-Si:H, a-Ge:H, a-C:H, a-SiC:H,
# a-SiGe:H). Their H caps a host atom through a covalent X-H bond; it is not the
# H- anion of the metal hydrides (LiH, MgH2, NaAlH4, TiH2).
_HYDROGENATED_NETWORK_HOSTS = frozenset({"C", "Si", "Ge"})
def _hydrogenated_host(composition) -> dict | None:
"""The H-free host of a hydrogenated group-IV network, else None.
A network is C, Si and/or Ge with H and nothing else, with at most one H
per host atom, so that each host atom keeps at least three of its four
bonds in the network: a-Si:H (Si64H8), a-C:H up to the polymer-like 50 %
H, a-SiC:H. More H than that makes chains and molecules (polyethylene,
SiH4), and any other element (O, N, F, a metal) leaves the compound to
the class rules for that element.
``composition`` may be a mapping of counts (preferred) or a bare set of
symbols, in which case one of each is assumed.
"""
counts = (dict(composition) if hasattr(composition, "items")
else {e: 1 for e in composition})
host = {e: n for e, n in counts.items() if e != "H"}
n_h = counts.get("H", 0)
if (n_h <= 0 or not host or not set(host) <= _HYDROGENATED_NETWORK_HOSTS
or n_h > sum(host.values())):
return None
return host
[docs]
def auto_target_cn(composition: dict) -> tuple[dict | None, int]:
"""
Auto-detect target coordination numbers and tolerance from composition.
Returns (target_cn, cn_tolerance) tuple.
Rules:
- Metalloids (Si, Ge, B): CN=4, tolerance=0 (strict tetrahedral)
- Pure Si/Ge: CN=4, tolerance=0
- Hydrogenated networks (a-Si:H, a-C:H, a-SiC:H): C, Si, Ge CN=4 and
H CN=1 (one X-H bond), tolerance=0
- Metals in nitrides: CN=4, tolerance=0 (wurtzite tetrahedral)
- Metals in halides/sulfides: CN=6, tolerance=0 (strict octahedral)
- Metals in oxides: CN=5, tolerance=1 (flexible 4-6 range)
- Nonmetal cations (:func:`cation_nonmetals`): the ligand count of their
oxoanion, 4 in PO4 3-, SO4 2- and ClO4 -, 3 in CO3 2-, NO3 - and IO3 -,
and 1 for the H of a hydroxide
Parameters
----------
composition : dict
e.g. {"In": 32, "O": 48}
Returns
-------
(target_cn, cn_tolerance) : tuple
target_cn: dict or None
cn_tolerance: int (0=strict, 1=flexible)
"""
centres = cation_nonmetals(composition)
target_cn, tol = _class_target_cn(composition, centres)
if centres:
target_cn = dict(target_cn or {})
for s in centres:
target_cn[s] = _nonmetal_cation_cn(s, composition)
return target_cn, tol
# Ligand count of a nonmetal cation in its highest oxidation state: the O of
# its oxoanion (trigonal CO3 2- and NO3 -; tetrahedral PO4 3-, SO4 2-,
# SeO4 2-, ClO4 -, BrO4 -, IO4 -) or the one O-H bond of a hydroxide. Each
# lone pair of a lower state takes a ligand's place (pyramidal SO3 2-,
# SeO3 2-, ClO3 -, BrO3 -, IO3 -; bent NO2 -).
_NONMETAL_CATION_CN = {"H": 1, "C": 3, "N": 3, "P": 4, "S": 4, "Se": 4,
"Cl": 4, "Br": 4, "I": 4}
def _nonmetal_cation_cn(sym: str, composition: dict) -> int:
"""Target CN of a nonmetal cation (see :data:`_NONMETAL_CATION_CN`)."""
top = _max_positive_os(sym)
# an off-stoichiometry cell can balance to more than the top state (S8O28)
ox = min(infer_oxidation_state(sym, composition) or top, top)
return max(1, _NONMETAL_CATION_CN.get(sym, 4) - (top - ox) // 2)
# Target CN of a one-element low-coordinated solid: 8 - N for the pnictogens
# and chalcogens, 3 for sp2 carbon.
_LOW_COORDINATED_CN = {"P": 3, "As": 3, "Sb": 3, "Bi": 3,
"S": 2, "Se": 2, "Te": 2, "C": 3}
def _class_target_cn(composition: dict,
centres: frozenset) -> tuple[dict | None, int]:
"""The per-class rules of :func:`auto_target_cn`, for every element but
the nonmetal cations ``centres``, which are neither anions nor metal
cations."""
elems = set(composition) - centres
anions = elems & NONMETALS
cations = elems - anions
# Check compound class first (handles semiconductors, chalcogenides)
cls = _classify_compound(composition)
if cls == "tetrahedral_element":
# Pure group IV (Si, Ge, SiGe) — tetrahedral CN=4
target_cn = {s: 4 for s in elems}
return target_cn, 0
if cls == "low_coordinated_element":
# Pure P, As, Sb, Bi (CN 3), S, Se, Te (CN 2): the 8-N rule of a
# covalent element. Pure C takes 3, the sp2 network its packing factor
# assumes. Before this branch As and Sb fell through to the metal rule
# (CN 8) and Se, P and C got no target.
(only_el,) = elems
cn = _LOW_COORDINATED_CN.get(only_el)
return ({only_el: cn} if cn else None), 0
if cls == "hydrogenated_network":
# a-Si:H, a-C:H: tetravalent hosts, each H bonded to one of them.
# Strict, so no host is placed next to an H that has its bond.
return {s: 1 if s == "H" else 4 for s in elems}, 0
if cls == "pnictide":
# III-V compounds (GaAs, InP, InAs) — tetrahedral CN=4
target_cn = {s: 4 for s in elems}
return target_cn, 0
if cls in ("chalcogenide", "chalcogenide_glass"):
# Chalcogenides (ZnS, Sb2Te3, CdTe, GeTe) and chalcogenide glasses
chalcogens = elems & {"S", "Se", "Te"}
target_cn = {}
for s in elems:
if s not in chalcogens:
if s in _ALWAYS_TETRAHEDRAL or s in METALLOIDS:
target_cn[s] = 4
else:
target_cn[s] = 6
return (target_cn if target_cn else None), 0
if cls in {"covalent_carbide", "transition_metal_carbide"}:
# Covalent carbides (SiC, B4C) — metalloid CN=4
# Transition-metal carbides (TiC, WC, ZrC) — transition-metal CN=6 (rocksalt)
target_cn = {}
for s in elems:
if s != "C":
if s in _ALWAYS_TETRAHEDRAL or s in METALLOIDS:
target_cn[s] = 4
else:
target_cn[s] = 6
return (target_cn if target_cn else None), 0
if cls == "hydride":
# Hydrides (LiH, MgH2) — metal CN=6
target_cn = {}
for s in elems:
if s != "H":
target_cn[s] = 6
return (target_cn if target_cn else None), 0
if cls == "boride":
# Borides (TiB2, LaB6) — metal CN=6
target_cn = {}
for s in elems:
if s != "B":
target_cn[s] = 6
return (target_cn if target_cn else None), 0
if cls == "alloy":
# Metal alloys / intermetallics — typical CN=8-12
# Use CN=8 (BCC-like) with tolerance=2 for flexibility
target_cn = {s: 8 for s in elems}
return target_cn, 2
# Dioxides / high-valent oxides — use the crystallographic cation
# coordination (octahedral rutile, 8-coordinate fluorite) instead of the
# generic oxide CN=5. The anion CN follows from stoichiometry (e.g. Ir=6 in
# IrO2 gives O=3). Tolerance 1 allows the amorphous spread (e.g. Ir 5-6).
if cls == "rutile_dioxide":
return {s: 6 for s in cations}, 1 # octahedral MO2: TiO2, IrO2, RuO2...
if cls == "fluorite_dioxide":
return {s: 8 for s in cations}, 1 # 8-coordinate MO2: ZrO2, HfO2, CeO2
if cls == "high_valent_oxide":
return {s: 6 for s in cations}, 1 # octahedral M(V)/M(VI): Nb2O5, WO3...
if not anions:
# Pure metals — typical CN=12 (FCC/HCP) or 8 (BCC)
target_cn = {s: 8 for s in elems}
return target_cn, 2
# Determine anion context
chalcogens = elems & {"S", "Se", "Te"}
pnictogens = elems & {"P", "As", "Sb"}
has_nitrogen = "N" in anions
has_oxide = "O" in anions
has_halide = bool(anions & _HALOGENS)
has_chalcogen = bool(chalcogens)
has_halide_sulfide = has_halide or has_chalcogen
target_cn = {}
if has_oxide and not has_nitrogen and not has_halide_sulfide:
# Pure oxides: metals CN=5 (flexible 4-6), metalloids CN=4 (strict)
has_metal_cation = any(s not in _ALWAYS_TETRAHEDRAL and s not in METALLOIDS for s in cations)
has_high_valent = False
for s in cations:
if s in _ALWAYS_TETRAHEDRAL or s in METALLOIDS:
target_cn[s] = 4
else:
target_cn[s] = 5
# High-valent d0 cations (Sb5+, Nb5+, Ta5+, Mo6+, W6+) are
# octahedral in oxides regardless of the metal/metalloid routing
# above. Matters for MIXED-cation oxides that don't reach the
# high_valent_oxide class (which needs ALL cations >= 5):
# ZnSb2O6 must target Sb=6, not the metalloid default 4 —
# validated by relaxation (forced Sb=6 seeds relax 0.5 eV lower
# with fewer dangling O than the auto Sb=4 seeds).
os_val = infer_oxidation_state(s, composition)
if os_val is not None and os_val >= 5:
target_cn[s] = 6
has_high_valent = True
# tolerance=1 if metal cations present (flexible 4-6) or a
# high-valent octahedral target was set (flexible 5-7);
# pure metalloid oxides (SiO2, GeO2, B2O3): tolerance=0 (strict CN=4)
return target_cn, 1 if (has_metal_cation or has_high_valent) else 0
elif has_nitrogen and not has_oxide:
# Nitrides: all CN=4 strict
for s in cations:
target_cn[s] = 4
return target_cn, 0
elif has_halide and not has_oxide:
# Halides: CN=6 strict, metalloids CN=4
for s in cations:
if s in _ALWAYS_TETRAHEDRAL or s in METALLOIDS:
target_cn[s] = 4
else:
target_cn[s] = 6
return target_cn, 0
# Chalcogenides already handled above via cls == "chalcogenide"
else:
# Mixed or other: metals CN=5 flexible, metalloids CN=4.
# DELIBERATE: the high-valent CN=6 override from the pure-oxide
# branch is NOT applied here. In mixed-anion frameworks the large
# soft anions lower the cation coordination — validated for
# oxyhalides: Ta(V) in a-NaTaOCl4 is ~5-coordinate (trigonal
# bipyramidal; Tang et al., Angew. Chem. 2026, EXAFS+AIMD, and our
# MLIP relaxations agree), unlike octahedral Ta(V) in Ta2O5. Do not
# "fix" this asymmetry by adding the override below.
for s in cations:
if s in _ALWAYS_TETRAHEDRAL or s in METALLOIDS:
target_cn[s] = 4
else:
target_cn[s] = 5
return target_cn, 1
# Elemental solid densities (kg/m3).
# Source: CRC Handbook of Chemistry and Physics, 104th edition (2023).
# Gaseous elements (O, N, Cl, etc.) are excluded.
# (ASE's ase.data.covalent_radii — used elsewhere in this module — is
# sourced from Cordero et al. 2008 for elements 1-96.)
ELEMENTAL_DENSITIES = {
# Solid nonmetals / metalloids (kg/m3)
"B": 2340, "C": 2267, "P": 1823, "S": 2070,
"Se": 4809, "Te": 6240, "I": 4930,
# Metals
"Li": 534, "Be": 1850, "Na": 971, "Mg": 1738, "Al": 2700,
"Si": 2330, "K": 862, "Ca": 1550, "Sc": 2985, "Ti": 4507,
"V": 6110, "Cr": 7190, "Mn": 7440, "Fe": 7874, "Co": 8900,
"Ni": 8908, "Cu": 8960, "Zn": 7134, "Ga": 5904, "Ge": 5323,
"As": 5727, "Sr": 2640, "Y": 4472, "Zr": 6506,
"Nb": 8570, "Mo": 10220, "Cd": 8650, "In": 7310, "Sn": 7287,
"Sb": 6685, "Ba": 3510, "La": 6145, "Ce": 6770,
"Hf": 13310, "Ta": 16654, "W": 19250, "Pb": 11340, "Bi": 9807,
"Tl": 11850,
}
# ==============================================================================
# Scale factors
# ==============================================================================
# Applied to sum of radii to give minimum interatomic distance.
# Scientific basis:
# - Shannon (1976): ionic radii from ~900 oxide/fluoride crystal structures
# - AIRSS buildcell: uses 80-90% of equilibrium distances
SCALE_FACTORS = {
"ionic": 0.80, # M-X bonds (shorter -> better coordination)
"covalent": 0.80, # Directional bonds
"metallic": 0.85, # M-M (keep higher to prevent M-M clustering)
"anion_packing": 0.80, # Small anion packing (O-O, F-F, N-N)
"anion_packing_large": 0.70, # Large anion packing (Cl-Cl, Br-Br, I-I, S-S)
}
# Shannon ionic radius threshold for small vs large anions.
_LARGE_ANION_RADIUS_THRESHOLD = 1.5
# Maximum minsep caps (A). Prevents overly large distances
# that make random placement impossible.
_MAX_SAME_ELEMENT_MINSEP = 2.80 # M-M, X-X same-element
_MAX_IONIC_MINSEP = 3.00 # M-X ionic bonds
_MAX_ANION_MINSEP = 3.00 # X-X anion packing
def _cation_contact_factor(anion_cn: float) -> float:
"""d(M-M) / d(M-X) for two cations sharing an anion at its narrowest angle.
An anion with up to 6 cations can keep them 90 degrees apart (the shared
edge of two octahedra, sqrt(2)). More cations crowd closer: 70.5 degrees
across the edge of a cube (8, as in antifluorite Li2O / Li2S), 60 degrees
in a cuboctahedron or the hexagonal bipyramid of Li3N (beyond 8).
"""
if anion_cn <= 6:
return 2**0.5
if anion_cn <= 8:
return 2 * np.sin(np.radians(70.5288) / 2)
return 1.0
# Anions that bond to themselves when the cations cannot reduce them all to
# their nominal charge (Zintl-Klemm): the anion excess, 1 - supply / demand,
# above which the floor becomes the X-X bond. O needs more: a table cation
# below its real state (Ni2+ for the Ni3+ of LiNiO2, 0.25) must not read as
# a peroxide (Li2O2, 0.5).
_HOMOPOLAR_EXCESS = {"S": 0.10, "Se": 0.10, "P": 0.10, "C": 0.10, "N": 0.10,
"I": 0.10, "O": 0.40}
# ... and the scale on twice the Cordero radius: single bonds (S-S, Se-Se,
# P-P, I-I) at 0.80, the multiple bonds of C2 2-, N3 - and O2 - at 0.70.
_HOMOPOLAR_SCALE = {"C": 0.70, "N": 0.70, "O": 0.70}
_ANION_VALENCE = {**ANION_CHARGES, "C": -4, "P": -3}
# lone-pair cations of the chalcogenide glasses supply 3, not their top 5
_LONE_PAIR_VALENCE = {"As": 3, "Sb": 3, "Bi": 3}
# H-H in a metal hydride stays above ~2.1 A (Switendick), the geminal H of
# BH4 - or ReH9 2- above 1.9 A: 0.8 of the first
_HYDRIDE_HH_MINSEP = 0.8 * 2.1
def _homopolar_anion(counts: dict, anion_syms) -> str | None:
"""The anion that must bond to itself in this compound, or None.
The cations, at their top states (3 for As, Sb, Bi), cannot reduce every
anion to its nominal charge, so the least electronegative anion forms
X-X bonds: S2 2- in FeS2, Se chains in Ge20Se80 and S-rich As-S, C2 2- in
CaC2 and YCI, N3 - in NaN3, O2 2- in Li2O2, I3 - in CsI3. Where other
anions are present only S, Se, P or C qualify; an O-rich doped cell is
not a peroxide.
"""
if not anion_syms or set(counts) <= set(anion_syms):
return None # a pure element, not a compound
supply = sum(_LONE_PAIR_VALENCE.get(e, _max_positive_os(e)) * n
for e, n in counts.items() if e not in anion_syms)
demand = sum(-_ANION_VALENCE.get(a, -2) * counts[a] for a in anion_syms)
x = min(anion_syms, key=lambda a: PAULING_EN.get(a, 4.0))
if x not in _HOMOPOLAR_EXCESS or demand <= 0:
return None
if len(anion_syms) > 1 and x not in ("S", "Se", "P", "C"):
return None
return x if 1 - supply / demand > _HOMOPOLAR_EXCESS[x] else None
# ==============================================================================
# Bonding classification
# ==============================================================================
[docs]
def classify_bond(sym_a: str, sym_b: str, composition=None) -> str:
"""
Classify a pair of elements by bonding type.
Uses an element-type classification rule (nonmetal / metalloid / metal
membership) with a Pauling-electronegativity refinement: pairs that
the type rules would call "ionic" but whose Pauling electronegativity
difference is small (Δχ < 1.0) are reclassified as "covalent".
This catches III–V semiconductors (GaAs Δχ=0.37, AlAs ~0.6) and related
compounds where the chemistry is predominantly covalent and Shannon
ionic radii give unrealistically small inter-atomic distances.
Returns one of: "ionic", "covalent", "metallic".
With ``composition`` (a {symbol: count} mapping), the roles charge
balance gives the two elements in that compound come first. A nonmetal
cation (:func:`cation_nonmetals`: P in a phosphate, S in a sulfate, H in
a hydroxide) bonds to the anions like any other cation, so the pair is
"ionic". Two cations of a compound with anions only meet across an
anion: "cation-cation" when one is a nonmetal cation, and never "ionic"
otherwise, since the Δχ rule is for bonded pairs (Na-B in a borate and
K-Si in a silicate are then "covalent", as Na-Si already is).
"""
if composition is not None:
centres = cation_nonmetals(composition)
anions = {s for s in composition if s in NONMETALS and s not in centres}
if anions and sym_a not in anions and sym_b not in anions:
if sym_a in centres or sym_b in centres:
return "cation-cation"
bond_type = classify_bond(sym_a, sym_b)
return "covalent" if bond_type == "ionic" else bond_type
if sym_a in centres or sym_b in centres:
return "ionic"
a_is_nonmetal = sym_a in NONMETALS
b_is_nonmetal = sym_b in NONMETALS
a_is_metalloid = sym_a in METALLOIDS
b_is_metalloid = sym_b in METALLOIDS
# Helper: small Δχ ⇒ predominantly covalent regardless of type rules.
def _is_covalent_by_en():
en_a = PAULING_EN.get(sym_a)
en_b = PAULING_EN.get(sym_b)
if en_a is None or en_b is None:
return False
return abs(en_a - en_b) < 1.0
if a_is_metalloid and b_is_metalloid:
return "covalent"
if (a_is_metalloid and b_is_nonmetal) or (b_is_metalloid and a_is_nonmetal):
# Was unconditionally "ionic" — but BN, BC etc. are covalent.
return "covalent" if _is_covalent_by_en() else "ionic"
if a_is_metalloid or b_is_metalloid:
# Metal + metalloid (e.g.\ Ga + As). Was unconditionally "ionic" —
# but III-V semiconductors are predominantly covalent.
return "covalent" if _is_covalent_by_en() else "ionic"
if a_is_nonmetal != b_is_nonmetal:
# Metal + nonmetal (e.g.\ Ga + N). Was unconditionally "ionic" —
# check Δχ for borderline cases (GaN Δχ=1.23 stays ionic; carbides
# of low-EN metals like Ti-C Δχ=1.0 may flip).
return "covalent" if _is_covalent_by_en() else "ionic"
if a_is_nonmetal and b_is_nonmetal:
return "covalent"
return "metallic"
# ==============================================================================
# Radius lookup functions
# ==============================================================================
[docs]
def get_ionic_radius(sym: str, cn: int | None = None,
oxidation_state: int | None = None) -> float | None:
"""
Get Shannon ionic radius for an element.
Parameters
----------
sym : str
Element symbol.
cn : int, optional
Coordination number. If provided, use CN-specific radius.
Defaults to CN=6.
oxidation_state : int, optional
Specific oxidation state to look up (e.g.\\ ``5`` for Sb in
Sb_2O_5). If supplied and present in the Shannon table the
matching radius is returned. If supplied but absent, falls
through to the default selection. When ``None`` (default),
anions return their most-negative state and cations return
their **highest** positive state - usually the right choice for
binary oxides / halides / nitrides where the higher state
dominates (TiO_2, V_2O_5, WO_3 etc.).
Returns
-------
float or None
"""
if sym not in SHANNON_IONIC_RADII:
return None
states = SHANNON_IONIC_RADII[sym]
if oxidation_state is not None and oxidation_state in states:
ox = oxidation_state
elif oxidation_state is not None and oxidation_state > 0:
# Requested cation state is missing from Shannon; pick the
# closest available positive state (preferring same-or-higher).
# This handles cases like ZrN, where charge balance suggests
# Zr(III) but Shannon only lists Zr(IV); using Zr(IV) is closer
# than dropping all the way to a covalent fallback.
positive = sorted(k for k in states.keys() if k > 0)
if not positive:
return None
ox = min(positive, key=lambda k: (abs(k - oxidation_state), -k))
elif sym in NONMETALS:
ox = min(states.keys()) # most-negative anion state
else:
positive = {k: v for k, v in states.items() if k > 0}
if not positive:
return None
ox = max(positive.keys()) # highest cation state when unspecified
cn_dict = states[ox]
target_cn = cn if cn is not None else 6
if target_cn in cn_dict:
return cn_dict[target_cn]
elif 6 in cn_dict:
return cn_dict[6]
else:
return next(iter(cn_dict.values()))
[docs]
def get_effective_radius(sym: str) -> float:
"""
Get the most appropriate radius for volume estimation.
- Metals -> metallic radius
- Anion-formers (O, S, Cl, ...) -> Shannon ionic radius
- Metalloids (Si, Ge, ...) -> metallic radius if available
- Fallback -> covalent radius
"""
if sym not in NONMETALS:
r = get_metallic_radius(sym)
if r is not None:
return r
if sym in NONMETALS and sym not in {"He", "Ne", "Ar", "Kr", "Xe", "Rn"}:
r = get_ionic_radius(sym)
if r is not None:
return r
return covalent_radii[atomic_numbers[sym]]
# ==============================================================================
# Minsep calculation
# ==============================================================================
[docs]
def default_minsep(symbols: list[str], scale: float = 0.85,
target_cn: dict | None = None) -> dict:
"""
Build a minsep dict for all element pairs.
Uses bonding-type-aware radii. When target_cn is provided, uses
CN-specific Shannon radii (e.g. CN=4 for Si in SiO2).
Pairs are classified in the compound (``classify_bond`` with the
composition), so a nonmetal that charge balance makes a cation, the P of
a phosphate or the C of a carbonate, is kept at bonding distance from
its anions (P-O 1.26 A, C-O 1.06 A) and away from the other cations.
Two metal cations meet across an anion, at sqrt(2) x their bond (two
octahedra sharing an edge) and at least their metallic contact. Where the
anions average more than 6 cations (Li2O, Li2S, Li3N, Cu2S) the cations
come closer than either, and the floor follows the narrower angle alone:
Li-Li 1.74 A in Li3N, against 2.11 A in the crystal. Beyond 12 (Fe3C,
Ni80P20) the metals touch, at their metallic contact.
Same-element anions are kept at packing distance unless the cations
cannot reduce them all to their nominal charge; then the least
electronegative one gets its X-X bond (S-S in FeS2, Se-Se in Ge20Se80,
C-C in CaC2: :func:`_homopolar_anion`). Two H+ may share an O (water),
hydride H- stays 1.68 A from H-, a covalent metal-P bond takes covalent
radii, and the O of a nitrate or carbonate keep their triangle edge.
A hydrogenated network (a-Si:H, a-C:H) has no anions: see
:func:`_hydrogenated_network_minsep`.
Parameters
----------
symbols : list of str
One symbol per atom, or a ``{symbol: count}`` mapping; the counts
decide which nonmetals act as cations.
scale : float
Fallback scale factor (default 0.85).
target_cn : dict, optional
e.g. {"Si": 4, "O": 2}. Uses CN-specific radii when available.
Returns
-------
dict
Minimum separations keyed as "A-B" with A <= B alphabetically.
"""
counts = Counter(symbols)
host = _hydrogenated_host(counts)
if host is not None:
return _hydrogenated_network_minsep(host, scale, target_cn)
unique = sorted(set(symbols))
minsep = {}
cn_map = target_cn or {}
centres = cation_nonmetals(counts)
anion_syms = [s for s in unique if s in NONMETALS and s not in centres]
has_anion = bool(anion_syms)
# Mean number of cations around an anion, from the cation targets (the
# automatic ones where none are given, 6 where neither is). Above 6 the
# compound is cation-rich (Li2O, Li2S, Li3N, Cu2S): its cation polyhedra
# share faces, and two cations meet across an anion at less than the 90
# degrees the M-M floor below assumes. Above 12 no anion can hold them
# all (Fe3C, Ni80P20): the metals touch each other.
anion_cn = 0.0
if has_anion:
gate_cn = dict(auto_target_cn(counts)[0] or {})
gate_cn.update(cn_map)
anion_cn = (sum(gate_cn.get(s, 6) * n for s, n in counts.items()
if s not in anion_syms)
/ sum(counts[s] for s in anion_syms))
homopolar = _homopolar_anion(counts, anion_syms)
def _cn_radius(sym):
if sym in centres:
# its shortest bond: the highest cation state at the lowest
# tabulated coordination (P5+ CN4, S6+ CN4, C4+ CN3, H+ CN1).
# Lower states bond longer (sulfite, chlorate), so this is a
# floor for them too, whatever target CN the placement uses.
states = SHANNON_IONIC_RADII[sym]
return min(states[max(k for k in states if k > 0)].values())
cn = cn_map.get(sym)
return get_ionic_radius(sym, cn=cn)
for i, s1 in enumerate(unique):
for s2 in unique[i:]:
key = f"{s1}-{s2}" if s1 <= s2 else f"{s2}-{s1}"
bond_type = classify_bond(s1, s2, counts)
sf = SCALE_FACTORS.get(bond_type, scale)
if s1 == s2 and s1 == homopolar:
# the X-X bond of the excess anions, not anion packing
r = covalent_radii[atomic_numbers[s1]]
sf_x = _HOMOPOLAR_SCALE.get(s1, SCALE_FACTORS["covalent"])
minsep[key] = float(2 * r * sf_x)
logger.info(
" minsep %s = %.2f A (anion-excess X-X bond: Cordero "
"%.3f, scale=%.2f)", key, minsep[key], r, sf_x
)
continue
if bond_type == "cation-cation":
# A nonmetal cation meets the other cations only across an
# anion. With another cation it can share a polyhedron edge
# (LiO4-PO4), so the floor is the right-angle contact used for
# M-M below. Two of the same element never share an anion at
# a narrow angle, so theirs is the two bonds end to end:
# P-P 2.67, S-S 2.58, C-C 2.24 A, above the X-X bond, so the
# placement cannot make one and compute_dimers still flags it.
sf = SCALE_FACTORS["metallic"]
r_anion = max((_cn_radius(a) or 0.0) for a in anion_syms)
if s1 == s2 == "H":
# except two H, which share the O of a water molecule
# at 104.5 degrees (H-H 1.52 A): the right-angle contact
# across the smallest anion, 1.23 A
r_small = min((_cn_radius(a) or 9.0) for a in anion_syms)
d_cc = 2**0.5 * (_cn_radius(s1) + r_small) * sf
elif s1 == s2:
d_cc = 2 * (_cn_radius(s1) + r_anion) * sf
else:
d_geom = [2**0.5 * (r + r_anion) * sf
for r in (_cn_radius(s1), _cn_radius(s2))
if r is not None]
d_cc = sum(d_geom) / len(d_geom)
minsep[key] = min(d_cc, _MAX_SAME_ELEMENT_MINSEP)
logger.info(
" minsep %s = %.2f A (cation-cation across the anion, "
"r_anion=%.3f, scale=%.2f)", key, minsep[key], r_anion, sf
)
continue
if bond_type == "ionic":
r1 = _cn_radius(s1)
r2 = _cn_radius(s2)
if r1 is not None and r2 is not None:
d_ionic = min((r1 + r2) * sf, _MAX_IONIC_MINSEP)
minsep[key] = d_ionic
logger.info(
" minsep %s = %.2f A (ionic: Shannon %.3f + %.3f, "
"scale=%.2f)", key, minsep[key], r1, r2, sf
)
continue
elif bond_type == "metallic":
r1 = get_metallic_radius(s1)
r2 = get_metallic_radius(s2)
if r1 is not None and r2 is not None:
d_metallic = (r1 + r2) * sf
if has_anion:
ri1 = _cn_radius(s1)
ri2 = _cn_radius(s2)
r_anion = max(
(_cn_radius(a) or 0.0) for a in anion_syms
) if anion_syms else None
# Geometric estimate for edge-sharing polyhedra:
# d(M-M) = sqrt(2) * d(M-X) * sf
# Uses the same scale factor as ionic/metallic (0.85).
# A cation-rich compound shares narrower edges and
# faces (see _cation_contact_factor).
factor = _cation_contact_factor(anion_cn)
if ri1 is not None and r_anion and r_anion > 0:
d_geom_1 = factor * (ri1 + r_anion) * sf
else:
d_geom_1 = 0.0
if ri2 is not None and r_anion and r_anion > 0:
d_geom_2 = factor * (ri2 + r_anion) * sf
else:
d_geom_2 = 0.0
d_geometric = (d_geom_1 + d_geom_2) / 2 if (
d_geom_1 > 0 and d_geom_2 > 0
) else max(d_geom_1, d_geom_2)
if anion_cn > 12:
# metal-rich: M-M is a metallic contact, as in
# the alloy (Ni-Ni 2.52 A in Ni80P20, Fe-Fe 2.49
# A in Fe3C)
minsep[key] = min(d_metallic,
_MAX_SAME_ELEMENT_MINSEP)
logger.info(
" minsep %s = %.2f A (M-M metal-rich, anion "
"CN=%.1f: metallic %.3f + %.3f, scale=%.2f)",
key, minsep[key], anion_cn, r1, r2, sf
)
continue
if anion_cn > 6 and d_geometric > 0:
# The metallic radius is the neutral atom's: Li+
# sits 2.11 A from Li+ in Li3N, 2.31 A in Li2O,
# against 3.04 A in Li metal.
minsep[key] = min(d_geometric,
_MAX_SAME_ELEMENT_MINSEP)
logger.info(
" minsep %s = %.2f A (M-M cation-rich, anion "
"CN=%.1f: %.3f x d(M-X), met=%.2f unused, "
"scale=%.2f)", key, minsep[key], anion_cn,
factor, d_metallic, sf
)
continue
d_mm = min(max(d_metallic, d_geometric),
_MAX_SAME_ELEMENT_MINSEP)
minsep[key] = d_mm
source = "metallic" if d_metallic >= d_geometric else "geometric"
capped = " (capped)" if d_mm < max(d_metallic, d_geometric) else ""
logger.info(
" minsep %s = %.2f A (M-M ionic context, %s: "
"met=%.2f, geom=%.2f, scale=%.2f)%s",
key, minsep[key], source,
d_metallic, d_geometric, sf, capped
)
else:
minsep[key] = min(d_metallic, _MAX_SAME_ELEMENT_MINSEP)
logger.info(
" minsep %s = %.2f A (metallic: %.3f + %.3f, "
"scale=%.2f)", key, minsep[key], r1, r2, sf
)
continue
elif bond_type == "covalent":
a_is_metalloid = s1 in METALLOIDS
b_is_metalloid = s2 in METALLOIDS
if a_is_metalloid or b_is_metalloid:
# Use metallic radii if available (better for pure elements),
# fall back to covalent radii
r1 = get_metallic_radius(s1) or covalent_radii[atomic_numbers[s1]]
r2 = get_metallic_radius(s2) or covalent_radii[atomic_numbers[s2]]
d_covalent = (r1 + r2) * sf
# In anion context (oxide, nitride, etc.), metalloid-metalloid
# distances are much larger due to bridging anions (Si-O-Si).
# Use geometric estimate: sqrt(2) * d(X-anion)
if has_anion and s1 == s2:
ri = _cn_radius(s1)
r_anion = max(
(_cn_radius(a) or 0.0) for a in anion_syms
) if anion_syms else None
if ri is not None and r_anion and r_anion > 0:
d_geometric = 2**0.5 * (ri + r_anion) * sf
d_xx = min(max(d_covalent, d_geometric),
_MAX_SAME_ELEMENT_MINSEP)
minsep[key] = d_xx
source = "covalent" if d_covalent >= d_geometric else "geometric"
capped = " (capped)" if d_xx < max(d_covalent, d_geometric) else ""
logger.info(
" minsep %s = %.2f A (metalloid anion context, %s: "
"cov=%.2f, geom=%.2f, scale=%.2f)%s",
key, minsep[key], source,
d_covalent, d_geometric, sf, capped
)
continue
minsep[key] = d_covalent
logger.info(
" minsep %s = %.2f A (covalent/metalloid: "
"%.3f + %.3f, scale=%.2f)",
key, minsep[key], r1, r2, sf
)
continue
# Pure single-element nonmetal composition (e.g.\ pure C,
# pure S, pure P). This is a genuine covalent homo-element
# bond, not an anion-anion contact in a compound, so use
# Cordero covalent radii rather than the much larger
# Shannon anion radius. Heteroelement nonmetal pairs and
# same-element pairs in compounds (O-O in TiO_2, Cl-Cl in
# LiCl, ...) still go through anion-packing below.
if s1 == s2 and len(unique) == 1:
r = covalent_radii[atomic_numbers[s1]]
d_cov_homo = 2 * r * sf
minsep[key] = min(d_cov_homo, _MAX_SAME_ELEMENT_MINSEP)
logger.info(
" minsep %s = %.2f A (covalent homo-element: "
"Cordero %.3f, scale=%.2f)",
key, minsep[key], r, sf
)
continue
if "P" in (s1, s2) and "P" in anion_syms and (
s1 not in NONMETALS or s2 not in NONMETALS):
# A metal-P bond of small dchi (Ni2P, CoP, Ni80P20) is
# covalent: the 2.12 A P3- radius put Ni-P at 2.25 A
# against 2.2-2.3 A. Metallic + Cordero radii, as for a
# metalloid (Ni-P 1.85 A).
m = s2 if s1 == "P" else s1
r_m = get_metallic_radius(m) or covalent_radii[atomic_numbers[m]]
r_p = covalent_radii[atomic_numbers["P"]]
minsep[key] = float((r_m + r_p) * sf)
logger.info(
" minsep %s = %.2f A (covalent M-P: %.3f + Cordero "
"%.3f, scale=%.2f)", key, minsep[key], r_m, r_p, sf
)
continue
r1 = _cn_radius(s1)
r2 = _cn_radius(s2)
if r1 is not None and r2 is not None:
avg_r = (r1 + r2) / 2
if avg_r >= _LARGE_ANION_RADIUS_THRESHOLD:
sf_anion = SCALE_FACTORS.get("anion_packing_large", 0.70)
label = "large anion"
else:
sf_anion = SCALE_FACTORS.get("anion_packing", 0.80)
label = "small anion"
minsep[key] = min((r1 + r2) * sf_anion, _MAX_ANION_MINSEP)
if s1 == s2:
# Two O of one nitrate or carbonate sit an edge of
# its triangle apart, sqrt(3) x the bond: 2.17 and
# 2.22 A, just under the packing floor. Keep them at
# 0.85 of that edge where packing would exclude it.
for c in centres:
if _nonmetal_cation_cn(c, counts) != 3:
continue
edge = 3**0.5 * (_cn_radius(c) + r1)
if 0.95 * edge < minsep[key]:
minsep[key] = 0.85 * edge
label = f"{c}{s1}3 triangle-edge"
if s1 == s2 == "H" and "H" in anion_syms:
minsep[key] = min(minsep[key], _HYDRIDE_HH_MINSEP)
label = "Switendick hydride"
logger.info(
" minsep %s = %.2f A (%s packing: Shannon "
"%.3f + %.3f, scale=%.2f)",
key, minsep[key], label, r1, r2, sf_anion
)
continue
# Fallback: covalent radii (with cap)
r1 = covalent_radii[atomic_numbers[s1]]
r2 = covalent_radii[atomic_numbers[s2]]
d_fallback = (r1 + r2) * scale
if s1 == s2:
d_fallback = min(d_fallback, _MAX_SAME_ELEMENT_MINSEP)
else:
d_fallback = min(d_fallback, _MAX_IONIC_MINSEP)
minsep[key] = d_fallback
logger.info(
" minsep %s = %.2f A (fallback: covalent %.3f + %.3f, "
"scale=%.2f)", key, minsep[key], r1, r2, scale
)
return minsep
def _hydrogenated_network_minsep(host: dict, scale: float,
target_cn: dict | None) -> dict:
"""Minsep of a hydrogenated group-IV network (a-Si:H, a-C:H, a-SiC:H).
The H-free host keeps the floors it has on its own (Si-Si 1.87 A as in
a-Si, C-C 1.22 A as in a-C, the C-C anion packing of SiC). Each H is
bonded to a host atom at the covalent radii used for any covalent pair
(C-H 0.86, Si-H 1.18 A, 0.8 of the bond), and two H come no closer than
on one host atom at a right angle (H-H 1.21 A in a-C:H, 1.67 A in
a-Si:H), below the geminal pair of a CH2 (1.78 A) or SiH2 (2.42 A) and
far above the 0.74 A of an H2 molecule. As anions, C and H were kept
2.24 A from everything, so no C-C or C-H bond could be placed.
"""
minsep = default_minsep(host, scale=scale, target_cn=target_cn)
sf = SCALE_FACTORS["covalent"]
r_h = covalent_radii[atomic_numbers["H"]]
x_h = []
for s in sorted(host):
r = get_metallic_radius(s) or covalent_radii[atomic_numbers[s]]
key = f"{s}-H" if s <= "H" else f"H-{s}"
minsep[key] = (r + r_h) * sf
x_h.append(minsep[key])
logger.info(" minsep %s = %.2f A (covalent X-H: %.3f + %.3f, "
"scale=%.2f)", key, minsep[key], r, r_h, sf)
minsep["H-H"] = 2**0.5 * min(x_h)
logger.info(" minsep H-H = %.2f A (two H on one host atom, right "
"angle)", minsep["H-H"])
return minsep
# ==============================================================================
# Density / cell volume estimation
# ==============================================================================
[docs]
def estimate_density(composition: dict, amorphous_factor: float = 0.80) -> float | None:
"""
Estimate density from elemental solid densities.
Returns density in g/cm3, or None if any element lacks density data.
"""
total_mass = 0.0
total_vol = 0.0
for sym, count in composition.items():
dens = ELEMENTAL_DENSITIES.get(sym)
if dens is None or dens == 0:
return None
dens_gcm3 = dens / 1000.0
mass = atomic_masses[atomic_numbers[sym]] * count
total_mass += mass
total_vol += mass / dens_gcm3
if total_vol == 0:
return None
crystal_density = total_mass / total_vol
return crystal_density * amorphous_factor
# Packing factors for density estimation by material class.
#
# The same formula V_cell = sum(sphere_vol) / f_pack is used for all
# classes, but the radius type underneath the spheres differs by
# bonding regime (see :func:`_radius_for_density`):
# - Ionic (oxides, halides, nitrides, hydrides) -> Shannon ionic CN=6
# - typical f_pack range 0.50-0.58
# - Covalent (group-IV, pnictides, chalcogenides, borides, carbides)
# -> Cordero covalent. Covalent radii reflect bond lengths, so
# the sphere is much smaller than the actual atomic volume; f_pack
# is correspondingly smaller (~0.20-0.40).
# - Metallic (alloys) -> Goldschmidt metallic. Metallic radii are
# close to atomic centres in close-packed metals; f_pack ~0.65-0.75.
#
# Values calibrated to give experimental amorphous densities for
# representative systems in each class.
PACKING_FACTORS = {
# Ionic — Shannon ionic radii
"network_forming_oxide": 0.50, # SiO2, GeO2, B2O3 (network formers)
"metal_oxide": 0.52, # In2O3, Al2O3, Ga2O3 (and fluorite ZrO2/HfO2/CeO2)
"rutile_dioxide": 0.66, # TiO2, SnO2, RuO2, IrO2, OsO2 — dense rutile-type MO2
"fluorite_dioxide": 0.63, # ZrO2, HfO2, CeO2 — 8/7-coordinate fluorite/baddeleyite MO2
"high_valent_oxide": 0.60, # V2O5, Nb2O5, Ta2O5, Sb2O5, MoO3, WO3 — cation OS>=5, dense
"halide": 0.58, # Li2ZrCl6, LiF
"oxyhalide": 0.56, # NaTaOCl4, BiOCl, LaOCl, ZrOCl2 — halide framework
# with O substitution/bridging. NOMINAL value only:
# estimate_cell_length interpolates the actual pf
# between metal_oxide (0.52) and halide (0.58) by
# halogen fraction via _oxyhalide_packing_factor
# (NaTaOCl4 -> 0.568, rho 2.68 -> 2.93); this entry
# is the class's documented midpoint / fallback.
"large_cation_nitride": 0.52, # large-cation nitrides (ZrN, HfN, ScN). Kept low
# DELIBERATELY (verified 2026-07-31): placement
# saturates at ~73-82% of crystal density under
# the minsep constraints regardless of pf —
# raising it just triggers auto-expand overshoot
# and LOWERS the final density (ZrN 82%→73% at
# 0.58). Not fixable by SC/min-CN tweaks either
# (same saturation with both disabled). Use
# --target-density or post-generation NPT /
# cell relaxation when the density matters.
"small_cation_nitride": 0.62, # small-cation nitrides (AlN, GaN, Si3N4, TiN) —
# small cations fit the N sublattice, so they pack
# denser and still place. Cation-radius-gated.
"hydride": 0.55, # LiH, MgH2, NaAlH4
# Covalent — Cordero covalent radii
"tetrahedral_element": 0.30, # Si, Ge, SiGe, B (tetrahedral covalent network)
"low_coordinated_element": 0.28, # a-C, a-P, a-As, a-Sb, a-S, a-Se, a-Te (chain/layer, Cordero)
"pnictide": 0.32, # GaAs, InP, InAs (III-V compounds)
"chalcogenide": 0.30, # ZnS, CdTe, GeTe (II-VI / IV-VI)
"chalcogenide_glass": 0.23, # GeS2, GeSe2, As2S3, As2Se3: corner-sharing
# network glasses (sulfides/selenides of the
# network formers); 0.30 over-predicts them by
# 20-50 % (benchmark 2026-09-28)
"boride": 0.60, # TiB2, MgB2, ZrB2 — Goldschmidt cation +
# Cordero B (carbide-style; calibrated
# 2026-07-31: diborides land 80-84 % of
# crystal; was 56-89 % with all-Cordero/0.50)
"small_cation_oxide": 0.35, # BeO — covalent wurtzite oxide network;
# Cordero radii for both atoms (Be 0.96,
# O 0.66) + 0.35 -> rho ~2.96 vs measured
# amorphous 3.01 g/cm3 (ionic model: 1.42)
"hydrogenated_network": 0.30, # a-Si:H, a-Ge:H, a-C:H, a-SiC:H. NOMINAL
# value only: get_packing_factor uses the
# class of the H-free host (tetrahedral_element 0.30,
# low_coordinated_element 0.28 for a-C:H,
# covalent_carbide 0.32), so the estimate
# runs continuously into a-Si / a-C / SiC
# as the H goes to zero; see _NETWORK_H_RADIUS
# Carbides split by cation chemistry (2026-05-11) — see
# _TRANSITION_METAL_CARBIDE_CATIONS and _classify_compound for the routing.
"covalent_carbide": 0.32, # SiC, B4C — Cordero radii, open network
"transition_metal_carbide": 0.60, # TiC, WC, ZrC — rocksalt-like; uses
# Goldschmidt metallic for cation,
# Cordero for C (see _radius_for_density)
# Metallic — Goldschmidt metallic radii. Lowered from 0.70 to 0.60
# (2026-05-10) so random sequential placement has enough headroom:
# random close packing of equal hard spheres caps at ~0.64, and 0.70
# left no room for the placement attempt budget — Cu/Au/Ti/Ni etc.
# all failed single-call generate_random with auto density. At 0.60
# the cell is ~5% larger; predicted density for a-Cu is ~7.2 g/cm³
# vs experimental amorphous ~7.7-8.5 (~10% low — within design-target
# accuracy and easily overridden via target_density for production runs).
"alloy": 0.60, # NiTi, CuZr, brass, pure metals
"default": 0.52, # fallback (ionic)
}
# Radius (A) of H in the density estimate of a hydrogenated network; the host
# atoms keep their Cordero radii. At its Cordero radius (0.31 A) H adds no
# volume, so the H content would not move the density: a-C:H stayed near the
# 3.0 g/cm3 of H-free a-C from 10 to 50 at.% H. But each H replaces a host-host
# bond and brings its share of free volume (the H-decorated vacancies of
# a-Si:H), so the measured density falls steadily with H. At 0.90 A and the
# host packing factor one H adds about 10 A^3, which gives a-Si:H with 11 at.% H
# 2.30 g/cm3 (glow-discharge a-Si:H ~2.2 at ~10 %) and a-C:H 2.19 / 1.84 /
# 1.52 / 1.24 at 20 / 30 / 40 / 50 at.% H (hard a-C:H 1.6-2.2 at 30-40 %,
# polymer-like 1.2-1.6 at 40-50 %; Robertson, Mater. Sci. Eng. R 37, 129
# (2002); Casiraghi et al., PRB 72, 085401 (2005)).
_NETWORK_H_RADIUS = 0.90
# NOTE: ``AMORPHOUS_DENSITY_FACTORS`` and the associated elemental
# density-mixing branch were retired in favour of the unified
# class-aware sphere-packing model. ``estimate_density()`` (which
# implements the mixing rule) is kept as a public helper for users
# who want the mass-averaged crystal density directly.
def _classify_compound(composition: dict) -> str:
"""
Classify a compound by material class for packing factor selection.
Returns one of: "tetrahedral_element", "hydrogenated_network", "pnictide",
"chalcogenide", "network_forming_oxide", "small_cation_oxide",
"metal_oxide", "halide", "oxyhalide", "large_cation_nitride", "carbide", "hydride",
"boride", "alloy", "default".
"""
elems = set(composition.keys())
has_metal = any(s not in NONMETALS and s not in METALLOIDS for s in elems)
has_metalloid = any(s in METALLOIDS for s in elems)
# a nonmetal that charge balance makes a cation (S of a sulfate, Cl of a
# perchlorate, H of a hydroxide) is not an anion of the class tests below
centres = cation_nonmetals(composition)
anions = (elems & NONMETALS) - centres
chalcogens = (elems - centres) & {"S", "Se", "Te"}
pnictogens = elems & {"P", "As", "Sb"}
# Group IV: pure Si, Ge, SiGe (all metalloid, no chalcogen/pnictogen)
all_metalloid = all(s in METALLOIDS for s in elems)
if all_metalloid and not chalcogens and not pnictogens:
return "tetrahedral_element"
# Hydrogenated group-IV networks (a-Si:H, a-Ge:H, a-C:H, a-SiC:H): H caps
# a covalent network (see _hydrogenated_host). As the anion of the
# "hydride" branch below it gave a-Si:H the Si4+ ionic radius (15 g/cm3),
# and a-C:H fell to covalent_carbide with H as the anion (3.5 g/cm3).
if _hydrogenated_host(composition) is not None:
return "hydrogenated_network"
# Pure single-element covalent / semimetal solids that group-IV doesn't
# catch (a-Se, a-Te, a-As, a-Sb, a-P, a-S). A one-element system has no
# electronegativity difference, so it bonds covalently — it must use Cordero
# covalent radii, not the ionic "highest positive" fallback, which gives
# nonsense (e.g. pure As → As(5+) radius → ~150 g/cm3). This is the bonding
# analogue of the radius rule: ΔEN = 0 ⇒ covalent.
if len(elems) == 1:
(only_el,) = elems
if ((only_el in NONMETALS or only_el in METALLOIDS)
and only_el not in _DIATOMIC_NONMETALS):
return "low_coordinated_element"
# Metal-rich metal/metalloid glasses (Ni80P20, Fe80B20, Pd80Si20,
# Pd40Ni40P20): a transition or noble metal with a minority of the glass-
# forming metalloids P, B, Si, Ge, As, Sb and nothing else. These are
# metallic glasses, not pnictides or borides: the metal lattice sets the
# density (Goldschmidt radii, alloy packing). Gated on the metalloid
# fraction so III-V compounds (50 %) and MB2 / MB6 borides (>= 67 %) are
# untouched.
metals = {s for s in elems if s not in NONMETALS and s not in METALLOIDS}
minority = elems - metals
if (metals and minority and minority <= _GLASS_METALLOIDS
and not (metals & _S_BLOCK_METALS)):
f_metalloid = sum(composition[s] for s in minority) / sum(composition.values())
if f_metalloid < _METAL_RICH_ALLOY_MAX_METALLOID_FRAC:
return "alloy"
# Pnictides: III-V compounds (GaAs, InP, InAs, GaSb)
# Must have at least one non-pnictogen element (the group III metal)
# Pure pnictogens (Sb, As) are not III-V compounds
# P is both NONMETAL and pnictogen — exclude from anion check here
non_pnictogen = elems - pnictogens
other_anions = anions - {"P", "As"} # P/As can act as pnictogen anion
if pnictogens and non_pnictogen and not other_anions and not chalcogens:
return "pnictide"
# Chalcogenides (S, Se, Te) — check before oxides. Sulfides and selenides
# whose cations are all network formers (Ge, Si, As, Sb, B, P) are the
# corner-sharing chalcogenide glasses (GeS2, GeSe2, As2S3, As2Se3), far
# more open than the II-VI / IV-VI semiconductors; tellurides (GeTe,
# Sb2Te3, GST) keep the dense chalcogenide factor.
if chalcogens and "O" not in anions:
cations = elems - chalcogens
if (cations and chalcogens <= {"S", "Se"}
and cations <= _CHALCOGENIDE_GLASS_FORMERS):
return "chalcogenide_glass"
return "chalcogenide"
# Oxyhalides (O + a halogen, e.g. NaTaOCl4, BiOCl, LaOCl, ZrOCl2). These
# are halide-framework materials with O substitution/bridging, not open
# metal-oxide networks — the large halide anions set the packing. Routed
# BEFORE the oxide branch, which would otherwise catch any O-bearing
# compound and hand back the loose metal_oxide factor (0.52), sizing the
# cell too big and under-predicting the density. Uses Shannon ionic radii
# (same as oxide/halide) with a halide-leaning packing factor.
# Gated on the halogen being a framework anion, not a dopant: below
# _OXYHALIDE_MIN_HALOGEN_FRAC of the (halogen + O) pool the compound keeps
# its oxide routing (so F-doped TiO2 stays rutile_dioxide).
if "O" in anions and (anions & _HALOGENS):
if _halogen_fraction(composition) >= _OXYHALIDE_MIN_HALOGEN_FRAC:
return "oxyhalide"
# Dopant-level halogen: classify the halogen-free composition instead
# of falling through with the dopant still present — a stray F/Cl in
# the dict would otherwise corrupt the downstream checks that read the
# full composition (the charge-balance solve behind high_valent_oxide
# in particular: V2O5 + 1 F must stay high_valent_oxide, not degrade
# to metal_oxide). Recursion is safe: the stripped composition has no
# halogens, so this branch cannot re-enter.
stripped = {k: v for k, v in composition.items()
if k not in _HALOGENS}
return _classify_compound(stripped)
# Oxides
if "O" in anions:
cations = elems - anions
# Anomalous small-cation covalent oxides (BeO): the cation is too small
# and polarizing for the ionic model, so the oxide is a covalent
# tetrahedral network. Routed before the ionic branches; uses Cordero
# radii (see _radius_for_density) + a low packing factor.
if cations and all(c in _COVALENT_OXIDE_CATIONS for c in cations):
return "small_cation_oxide"
# High-valent oxides (all cations in oxidation state >= 5): V2O5, Nb2O5,
# Ta2O5, Sb2O5, MoO3, WO3, ... Small, highly-charged cations form short
# M-O bonds and pack denser than the generic oxide factor predicts, so
# they get a higher packing factor. Restricted to metal/metalloid
# cations (excludes molecular P2O5 / SO3 etc., and the oxoanion
# centres of a salt: NbOPO4 is high-valent through its Nb).
mm_cations = [c for c in cations if c not in NONMETALS]
if mm_cations:
os_list = [infer_oxidation_state(c, composition) for c in mm_cations]
if all(o is not None and o >= 5 for o in os_list):
return "high_valent_oxide"
# Rutile-type dioxides (MO2) pack much denser than light oxides; route
# them to a higher packing factor. Identified by the radius-ratio rule:
# a metal MO2 whose 4+ cation radius is below the rutile/fluorite cutoff
# (or a known covalent exception, e.g. PbO2). Fluorite dioxides
# (ZrO2, HfO2, CeO2 — larger cations) and metalloid oxides (SiO2, GeO2)
# fall through to metal_oxide / network_forming_oxide.
metal_cations = [c for c in cations
if c not in NONMETALS and c not in METALLOIDS]
if cations and len(metal_cations) == len(cations):
n_cat = sum(composition[c] for c in cations)
if composition.get("O", 0) == 2 * n_cat:
# MO2 dioxide: split rutile (small cation, octahedral) vs
# fluorite/baddeleyite (large cation, 8/7-coordinate) by the 4+
# cation radius, with PbO2 as a covalent rutile exception.
r4 = {c: get_ionic_radius(c, cn=6, oxidation_state=4)
for c in cations}
if all(c in _RUTILE_DIOXIDE_EXCEPTIONS or
(r4[c] is not None and r4[c] < _RUTILE_DIOXIDE_RMAX)
for c in cations):
return "rutile_dioxide"
if all(r4[c] is not None and r4[c] >= _RUTILE_DIOXIDE_RMAX
for c in cations):
return "fluorite_dioxide"
if has_metal:
return "metal_oxide"
elif has_metalloid or (len(elems) > 1
and all(c in NONMETALS for c in cations)):
# metalloid oxides (SiO2, GeO2, B2O3) and non-metal oxides, where
# every cation is a NONMETAL (P2O5, SO3), are covalent networks /
# molecular glasses, not "default".
# Pure oxygen (one element) still falls through to "default".
return "network_forming_oxide"
# Halides
if anions & _HALOGENS:
return "halide"
# Nitrides — small-cation nitrides (AlN, GaN, Si3N4, TiN, ...) pack densely
# and tolerate a higher packing factor; large-cation nitrides (ZrN, HfN,
# ScN) are placement-limited by the cation-cation distance, so they keep the
# safe default. Gate on cation radius (the same radius-rule idea as the
# rutile/fluorite dioxide split).
if "N" in anions:
cations = elems - anions
crs = [get_ionic_radius(c, cn=6,
oxidation_state=infer_oxidation_state(c, composition))
for c in cations]
if cations and all(cr is not None and cr < _SMALL_CATION_NITRIDE_RMAX
for cr in crs):
return "small_cation_nitride"
return "large_cation_nitride"
# Carbides — split by cation chemistry
# Transition-metal carbides (TiC, WC, ZrC, ...) form rocksalt-like dense
# structures with d-block transition metals and need a higher packing
# factor + Goldschmidt metallic radii for the cation.
# Covalent carbides (SiC, B4C, ...) keep Cordero radii + low packing
# factor. See _TRANSITION_METAL_CARBIDE_CATIONS for the cation set.
if "C" in anions:
cations = elems - anions
if cations & _TRANSITION_METAL_CARBIDE_CATIONS:
return "transition_metal_carbide"
return "covalent_carbide"
# Hydrides
if "H" in anions:
return "hydride"
# Borides (B as anion-like element with metals)
if "B" in elems and has_metal and not anions and not chalcogens:
return "boride"
# Metal alloys / intermetallics (no anions)
if not anions and not chalcogens and has_metal:
return "alloy"
return "default"
def _halogen_fraction(composition: dict) -> float:
"""Halogen fraction of the (halogen + O) anion pool.
The single definition shared by the oxyhalide classification gate and
the packing-factor interpolation, so the two can never disagree about
the same composition.
"""
n_hal = sum(composition.get(x, 0) for x in _HALOGENS)
n_o = composition.get("O", 0)
if n_hal + n_o == 0:
return 0.0
return n_hal / (n_hal + n_o)
def _oxyhalide_packing_factor(composition: dict) -> float:
"""Anion-ratio-interpolated packing factor for oxyhalides.
An oxyhalide's packing sits between its parent oxide and a pure halide,
controlled by how much of the anion sublattice is halogen vs oxygen.
Interpolate linearly on the halogen fraction of the anions:
f_pack = pf_base + (pf_halide - pf_base) * n_halogen/(n_halogen + n_O)
where ``pf_base`` is the packing factor of the **halogen-free base
class** (what `_classify_compound` returns with the halogens stripped),
not a hardcoded oxide value. This keeps the estimate continuous at the
dopant gate: a rutile-type MO2 just above the 10% halogen threshold
interpolates from rutile's 0.66 (pf ~0.65), not from metal_oxide's 0.52
— without this, one atom across the boundary jumped the predicted cell
volume by ~26%. Loose-oxide parents (NaTaOCl4 -> metal_oxide base,
pf 0.568; BiOCl -> 0.55) are unchanged, and the halogen-only endpoint
still recovers ``halide`` (0.58).
"""
frac_hal = _halogen_fraction(composition)
if frac_hal == 0.0: # defensive; not reachable for oxyhalide
return PACKING_FACTORS["oxyhalide"]
stripped = {k: v for k, v in composition.items() if k not in _HALOGENS}
base_cls = _classify_compound(stripped) if stripped else "default"
pf_base = PACKING_FACTORS.get(base_cls, PACKING_FACTORS["metal_oxide"])
pf_hal = PACKING_FACTORS["halide"]
return pf_base + (pf_hal - pf_base) * frac_hal
# Classes whose density radii mix bonding regimes: every cation uses
# Goldschmidt metallic, the mapped "anchor" anion-former uses Cordero
# covalent. Add new metal-covalent-anchor families (silicides, oxycarbides)
# here rather than copying the branch in _radius_for_density.
_MIXED_RADIUS_ANCHOR = {"transition_metal_carbide": "C", "boride": "B"}
[docs]
def get_packing_factor(cls: str, composition: dict | None = None) -> float:
"""Packing factor for a material class — composition-aware where needed.
The single dispatch point between the static :data:`PACKING_FACTORS`
table and classes whose packing depends on the composition itself:
``oxyhalide``, which interpolates between its halogen-free base class
and ``halide`` by halogen fraction, and ``hydrogenated_network``, which
takes the packing factor of its H-free host. New mixed-anion
families (oxynitrides, oxysulfides) should add their interpolation here
rather than special-casing :func:`estimate_cell_length`. Note the
obvious oxynitride interpolation is currently a no-op — metal_oxide and
nitride share pf 0.52 — so no oxynitride branch exists yet; add one only
with calibration data that motivates it.
"""
if cls == "oxyhalide" and composition is not None:
return _oxyhalide_packing_factor(composition)
if cls == "hydrogenated_network" and composition is not None:
host = _hydrogenated_host(composition)
if host is not None:
return PACKING_FACTORS[_classify_compound(host)]
return PACKING_FACTORS.get(cls, PACKING_FACTORS["default"])
_MAX_POSITIVE_OS = {"H": 1, "C": 4, "N": 5, "O": 2, "F": 1, "P": 5, "S": 6,
"Cl": 7, "Se": 6, "Br": 7, "Te": 6, "I": 7,
"Sb": 5} # Sb(V); Sb has no Shannon entry to read it from
def _max_positive_os(sym: str) -> int:
if sym in _MAX_POSITIVE_OS:
return _MAX_POSITIVE_OS[sym]
states = [k for k in SHANNON_IONIC_RADII.get(sym, {}) if k > 0]
return max(states) if states else 4
def _balances(counts: dict, charges: dict) -> bool:
"""Whether the elements of ``counts`` outside ``charges`` (anion: charge)
can supply exactly the anions' charge, each element in one of its
tabulated positive states. One state per element, as in
:func:`_solve_oxidation_states`: with mixed valence, 16 Cr at +4 and +6
would balance CrOOH as a hydride."""
demand = sum(-q * counts[a] for a, q in charges.items())
reachable = {0}
for e, n in counts.items():
if e in charges:
continue
states = ({k for k in SHANNON_IONIC_RADII.get(e, {}) if k > 0}
or {_max_positive_os(e)})
reachable = {r + n * k for r in reachable for k in states
if r + n * k <= demand}
if not reachable:
return False
return demand in reachable
[docs]
def anion_elements(composition) -> set:
"""Which elements act as ANIONS in this compound, decided by charge balance.
Membership of the anion table is not enough on its own: tellurium is the
anion in CdTe and the cation in TeO2, sulfur the anion in ZnS and the
cation in a sulfate, hydrogen the anion in LiH and the cation in a
hydroxide. So every element on the anion table starts as a candidate
anion, and while the cations present cannot supply the positive charge the
candidates demand, the least electronegative candidate is promoted to
cation. This reproduces the chemistry without a lookup table of exceptions:
La2O2S balances with sulfur as an anion and keeps it, H2SO4 does not and
promotes first hydrogen then sulfur, leaving the sulfate O as the anion.
A promotion happens only when it brings the compound CLOSER to charge
balance, which is what keeps real cells intact without any tolerance
setting. In a chalcogen-rich glass or a doped cell, promoting the major
anion would overshoot far past neutrality (Ge20S10Se70: a gap of 80 before,
480 after; F-doped silica: 2 before, 254 after), so every anion is kept:
the mixed chalcogen glasses keep Ge-Se and Ge-Te, F-doped silica keeps
Si-O, an O impurity in NaCl keeps Na-Cl and LiPON keeps P-N. Where a
promotion really is the chemistry, it improves the balance and happens:
TeO2, a tellurite, a sulfate, a nitrate, a hydroxide.
Known limitation: an element that is a cation in one site of a polyatomic
group and an anion in another cannot be both. In a thiosulfate the central
sulfur is cation-like and the terminal sulfur is an anion; the rule keeps
both as anions, so the S-O bonds are not counted. Rare in amorphous work.
``composition`` may be a mapping of counts (preferred) or a bare set of
symbols, in which case one of each is assumed.
"""
counts = (dict(composition) if hasattr(composition, "items")
else {e: 1 for e in composition})
anions = {e for e in counts if e in ANION_CHARGES}
if not anions:
return set()
def imbalance(anion_set):
"""Positive charge the cations can supply minus what the anions demand."""
cations = set(counts) - anion_set
return (sum(_max_positive_os(e) * counts[e] for e in cations)
- sum(-ANION_CHARGES[e] * counts[e] for e in anion_set))
# The most electronegative element is always an anion, so it is never a
# candidate: an off-stoichiometry cell (a random Ga16Zn16O48 composition, a
# defective model) must not end up with no anion at all.
def balances(anion_set):
return _balances(counts, {a: ANION_CHARGES[a] for a in anion_set})
for candidate in sorted(anions, key=lambda e: PAULING_EN.get(e, 2.0))[:-1]:
before = imbalance(anions)
if before >= 0 and balances(anions):
break # the cations already cover the anions
trial = anions - {candidate}
if before < 0 and abs(imbalance(trial)) < abs(before):
anions = trial
elif not balances(anions) and balances(trial):
# The highest states overstate what Ag, Tl, Pb, Mn or Cr supply:
# Ag2SO4 counts Ag+3 and ties (-4 against +4), CrOOH covers O and
# H with Cr+6. Only the sulfate (2 Ag+ + S6+) and the hydroxide
# (Cr3+ + H+) balance, so S and H are the cations.
anions = trial
return anions
def _radius_for_density(sym: str, cls: str,
composition: dict | None = None) -> float:
"""Pick the radius type appropriate for the composition's bonding character.
The density estimator uses sphere-packing on per-element spheres;
different bonding regimes need different radius tables so that the
sphere volume is a faithful proxy for the volume the atom occupies
in the corresponding amorphous structure:
* Ionic compounds (oxides, halides, nitrides, hydrides) -> Shannon
ionic radii at CN=6 (cation-anion contact length).
* Covalent compounds (group-IV, pnictides, chalcogenides, borides,
carbides) -> Cordero covalent radii (bond-length-based). In a
hydrogenated network (a-Si:H, a-C:H) the host atoms are Cordero and
H is :data:`_NETWORK_H_RADIUS`.
* Metallic alloys -> Goldschmidt metallic radii (close-packed
atomic centres).
"""
ionic_classes = {"network_forming_oxide", "metal_oxide", "rutile_dioxide",
"fluorite_dioxide", "high_valent_oxide", "halide",
"oxyhalide", "large_cation_nitride",
"small_cation_nitride", "hydride"}
covalent_classes = {"tetrahedral_element", "low_coordinated_element", "pnictide",
"chalcogenide", "chalcogenide_glass", "covalent_carbide",
"small_cation_oxide", "hydrogenated_network"}
if cls == "hydrogenated_network" and sym == "H":
return _NETWORK_H_RADIUS
# Infer oxidation state from charge balance when a composition is
# supplied; falls back to "highest positive" inside get_ionic_radius
# if inference returns None.
ox = infer_oxidation_state(sym, composition) if composition else None
# A non-metal acting as the CATION of an oxide (P in P2O5, S in SO3, Se in
# SeO2, Te in TeO2) must take its highest positive Shannon state; the anion
# default would hand back P3- (2.12 A) and shrink the density threefold.
# Gated on the element actually being a cation here, which charge balance
# decides (cation_nonmetals): the Cl of an oxychloride, the N of an
# oxynitride and the carbide C of an oxycarbide (SiOC) are anions and keep
# their anion radii.
if (ox is None and composition and "O" in composition and sym != "O"
and sym in NONMETALS and sym in cation_nonmetals(composition)):
positive = [k for k in SHANNON_IONIC_RADII.get(sym, {}) if k > 0]
if positive:
ox = max(positive)
else:
# promoted but no positive state tabulated: the anion radius would
# be ~4x too large, so use the covalent radius instead
return covalent_radii[atomic_numbers[sym]]
# Antimony is oxidation-state split: Sb(V) is d0 and octahedral with no lone
# pair, so it uses its small ionic radius like Nb(V)/Ta(V) (dense oxides such
# as Sb2O5). Sb(III) has a stereochemically active lone pair giving open
# structures (Sb2O3) where the small ionic radius badly over-predicts the
# density, so it falls back to the Cordero covalent radius.
if sym == "Sb":
return 0.60 if ox == 5 else covalent_radii[atomic_numbers["Sb"]]
if cls in ionic_classes:
r = get_ionic_radius(sym, cn=6, oxidation_state=ox)
if r is None:
r = covalent_radii[atomic_numbers[sym]]
return r
if cls in covalent_classes:
return covalent_radii[atomic_numbers[sym]]
if cls in _MIXED_RADIUS_ANCHOR:
# Dense metallic-bonded ceramics with one covalent anchor element:
# the anchor (C in carbides, B in borides) uses Cordero, every
# cation uses Goldschmidt metallic. Carbides: pf=0.60 reproduces
# a-TiC ≈ 4.0 g/cm³ (crystal 4.93) and a-WC ≈ 14 (15.6). Borides
# (calibrated 2026-07-31): pf=0.60 lands TiB2/ZrB2/MgB2 at 80-84 %
# of crystal — all-Cordero at pf 0.50 split them 56-89 %, unfixable
# by any single pf. Cage borides (LaB6-type) run ~100 % of crystal;
# use density_scale < 1 there if placement struggles.
if sym == _MIXED_RADIUS_ANCHOR[cls]:
return covalent_radii[atomic_numbers[sym]]
r = get_metallic_radius(sym)
if r is None:
r = covalent_radii[atomic_numbers[sym]]
return r
if cls == "alloy":
r = get_metallic_radius(sym)
if r is None:
r = covalent_radii[atomic_numbers[sym]]
return r
# Default: ionic with covalent fallback
r = get_ionic_radius(sym, cn=6, oxidation_state=ox)
if r is None:
r = covalent_radii[atomic_numbers[sym]]
return r
[docs]
def estimate_cell_length(composition: dict, target_density: float | None = None,
packing_factor: float | None = None,
density_scale: float = 1.0) -> float:
"""Estimate cubic cell length from composition via sphere packing.
A single class-aware sphere-packing model handles all amorphous
compositions:
.. math::
V_{\\text{cell}} = \\frac{\\sum_i N_i \\cdot \\tfrac{4}{3}\\pi r_i^3}
{f_{\\text{pack}} \\cdot s},
where :math:`r_i` is the radius for element :math:`i` chosen
according to the compound's bonding regime (Shannon ionic for ionic
compounds, Cordero covalent for covalent compounds, Goldschmidt
metallic for alloys; see :func:`_radius_for_density`).
:math:`f_{\\text{pack}}` is the class-aware packing fraction stored
in :data:`PACKING_FACTORS`, and :math:`s` is the user-supplied
``density_scale`` (default 1.0).
If ``target_density`` is supplied, the cell is sized directly from
it and ``density_scale`` / ``packing_factor`` are ignored.
"""
if density_scale <= 0:
raise ValueError(f"density_scale must be > 0 (got {density_scale}).")
total_mass = sum(atomic_masses[atomic_numbers[s]] * n
for s, n in composition.items())
if target_density is not None:
vol_cm3 = (total_mass / 6.022e23) / target_density
vol_A3 = vol_cm3 * 1e24
else:
cls = _classify_compound(composition)
if packing_factor is None:
packing_factor = get_packing_factor(cls, composition)
total_sphere_vol = 0.0
for sym, count in composition.items():
# a cation with a negative Shannon radius (H+) sits inside its
# anion's sphere and adds no volume
r = max(_radius_for_density(sym, cls, composition=composition), 0.0)
total_sphere_vol += count * (4.0 / 3.0) * np.pi * r ** 3
vol_A3 = total_sphere_vol / (packing_factor * density_scale)
return vol_A3 ** (1.0 / 3.0)
# ==============================================================================
# Auto-derivation log line
# ==============================================================================