"""
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": {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])
return (viol, 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 "covalent_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)
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 == "group_iv":
# Pure group IV (Si, Ge, SiGe) — tetrahedral CN=4
target_cn = {s: 4 for s in elems}
return target_cn, 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
# ==============================================================================
# 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.
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)
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 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:
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).
if ri1 is not None and r_anion and r_anion > 0:
d_geom_1 = 2**0.5 * (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 = 2**0.5 * (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)
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
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)
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
"covalent_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.
"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
"group_iv": 0.30, # Si, Ge, C (tetrahedral covalent network)
"elemental_semiconductor": 0.28, # a-Se, a-Te, a-As, a-Sb, a-P (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)
"covalent_network_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 (group_iv 0.30,
# elemental_semiconductor 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: "group_iv", "hydrogenated_network", "pnictide",
"chalcogenide", "covalent_oxide", "covalent_network_oxide",
"metal_oxide", "halide", "oxyhalide", "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 "group_iv"
# 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 "elemental_semiconductor"
# 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 "covalent_network_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 / covalent_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 "covalent_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 "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
[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.
for candidate in sorted(anions, key=lambda e: PAULING_EN.get(e, 2.0))[:-1]:
before = imbalance(anions)
if before >= 0:
break # the cations already cover the anions
trial = anions - {candidate}
if abs(imbalance(trial)) < abs(before):
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 = {"covalent_oxide", "metal_oxide", "rutile_dioxide",
"fluorite_dioxide", "high_valent_oxide", "halide",
"oxyhalide", "nitride",
"small_cation_nitride", "hydride"}
covalent_classes = {"group_iv", "elemental_semiconductor", "pnictide",
"chalcogenide", "chalcogenide_glass", "covalent_carbide",
"covalent_network_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
# ==============================================================================