The Math Behind Atomency

Every calculation the engine runs — written out in full

§ 1

Bond Detection & Distance Thresholds

When atoms are placed on the canvas the engine computes pairwise Euclidean distances and compares them against scaled van der Waals radius sums to decide whether a bond forms and at what order.

Euclidean Atom–Atom Distance
d(A, B) = √ (xB − xA)2 + (yB − yA)2
where x, y are canvas pixel coordinates
Computed for every unique pair (i, j) with i < j. O(n²) over all atoms on the canvas each time resolve() is called.
Proximity Screening — Potential Bond Threshold
r1 = vdW radius of atom A  (Å)
r2 = vdW radius of atom B  (Å)
dpotential = max(100 px,  (r1 + r2) × 1.35)
Pair enters the candidate list only if   d(A,B) ≤ dpotential
Bond Order from Distance Cutoffs
Single bond → d ≤ max(80, rsum × 1.15)
Double bond → d ≤ max(65, rsum × 0.85)  and valence allows
Triple bond → d ≤ max(45, rsum × 0.60)  and valence allows
Ionic if   |ΔEN| > 1.7   (Pauling scale)
Pairs sorted by |ΔEN| descending before bond assignment — ionic bonds form preferentially
Overlap Rejection
dmin = max(60 px,  (r1 + r2) × 0.9)
if d(A,B) < dmin   →   atoms too close, bond rejected
§ 2

Valence Electron Counting & Max Bond Capacity

The engine derives valence electrons from group number without lookup tables, then computes the maximum bonds each atom can form — including expanded octets for period 3+ elements.

Valence Electrons by Block
{
V = group numbers-block  (groups 1–2)
V = group − 10p-block  (groups 13–18)
V = from shells arrayd-block  (transition metals)
V = 3f-block  (lanthanides/actinides)
Maximum Bond Count
{
Period 1 (H, He) max_bonds = 2 − V
Period 2 (strict octet) max_bonds = max(0, 8 − V)
Period 3+, V ≥ 5 max_bonds = max(6, 8−V) + ⌊(V−4)/2⌋ × 2
Period 3+, V < 5 max_bonds = max(6, 8 − V)
Transition metals max_bonds = max(0, min(6, V + 2))
Example: Sulfur (V=6, period 3): max_bonds = max(6, 8−6) + ⌊2/2⌋×2 = 6+2 = 8, allowing SF₆. Carbon (V=4, period 2): max_bonds = max(0, 8−4) = 4. This is why CH₄ forms but CF₆ does not.
§ 3

Bond Type Classification

Every bond is classified using the Pauling electronegativity difference scale. Thresholds follow IUPAC conventions.

Pauling ΔEN Classification
ΔEN = |EN(A) − EN(B)|
1
ΔEN < 0.5 → nonpolar covalent
2
0.5 ≤ ΔEN < 1.7 → polar covalent  (δ⁺/δ⁻ arise)
3
ΔEN ≥ 1.7 → ionic  (near-complete electron transfer)
Additional rule: metal + nonmetal pair forces ionic regardless of ΔEN
H₂O: EN(O) = 3.44, EN(H) = 2.20 → ΔEN = 1.24 → polar covalent.   NaCl: EN(Na) = 0.93, EN(Cl) = 3.16 → ΔEN = 2.23 → ionic.
§ 4

Partial Charge (δ⁺ / δ⁻) — Pauling Ionic Character

For each bond the engine calculates the fractional electron transfer using Pauling's empirical exponential formula, then accumulates contributions across all bonds on the atom.

Pauling Partial Ionic Character
δ = 1 − e−ΔEN² / 4
where ΔEN = |EN(A) − EN(B)|
{
A gets +δif EN(A) < EN(B)  (loses electron density)
A gets −δif EN(A) > EN(B)  (gains electron density)
For bond order n:
δcontribution = ±(1 − e−ΔEN²/4) × n
δA = ∑bonds δcontribution,i
Worked Example — Water (H₂O)
EN(O) = 3.44,  EN(H) = 2.20  →  ΔEN = 1.24
δ = 1 − e−1.24² / 4 = 1 − e−0.3844 = 1 − 0.6809 = 0.319
Two O–H bonds (order = 1):
δ(O) = −0.319 × 2 = −0.64
δ(H) = +0.319 each = +0.32
§ 5

Formal Charge & Lone Pair Estimation

Formal charge is computed from standard electron-bookkeeping after bond resolution. Lone pairs follow from remaining valence electrons.

Formal Charge
V = valence electrons of the free atom
L = lone-pair electrons (non-bonding)
B = bonding electrons shared with atom  (= bond order × 2, then halved)
FC = V − L − B 2
Lone Pair Calculation
Vavailable = V − FC
Bused = ∑ bond.order  (over all bonds on atom)
eremaining = Vavailable − Bused
lone_pairs = max(0,  eremaining2 )
Example — oxygen in H₂O (V = 6, two single bonds):
Vavail = 6 − 0 = 6  |  Bused = 1 + 1 = 2  |  erem = 4
lone_pairs = 42 = 2 ✓
§ 6

VSEPR Geometry Classification

Molecular geometry is derived from the electron-domain count using the full VSEPR decision tree. Bond angles are ideal angles adjusted for lone pair compression.

Steric Number
SN = bonding_pairs + lone_pairs
SN=2, LP=0 → linear   180°
SN=3, LP=0 → trigonal planar   120°
SN=3, LP=1 → bent   ~119°
SN=4, LP=0 → tetrahedral   109.5°
SN=4, LP=1 → trig. pyramidal   107°
SN=4, LP=2 → bent   104.5°
SN=5, LP=0 → trig. bipyramidal   90°/120°
SN=5, LP=1 → see-saw   ~86°/102°
SN=5, LP=2 → T-shape   ~87.5°
SN=6, LP=0 → octahedral   90°
SN=6, LP=1 → square pyramidal   ~85°
SN=6, LP=2 → square planar   90°
Lone Pair Compression
θactual = θideal − 2.5° × LP
Each lone pair compresses bond angles by ~2.5° due to greater nuclear repulsion
§ 7

Dipole Moment Calculation

Bond dipoles are computed from a power-law fit calibrated to real diatomic molecules. Net molecular dipole is then calculated by vector summation, accounting for geometry and lone pair contributions.

Bond Dipole — Power-Law Fit
μbond = 1.121 × ΔEN0.924 Debye
ΔEN = |EN(A) − EN(B)|  |  if ΔEN ≤ 0.01 → μ = 0 (nonpolar)
Calibration:
HF: ΔEN = 1.78  →  1.121 × 1.780.924 = 1.898 D  (exp: 1.82 D)
HCl: ΔEN = 0.96  →  1.121 × 0.960.924 = 1.080 D  (exp: 1.08 D ✓)
Net Molecular Dipole — Vector Sum
Linear (θ = 180°):   μnet = |μ₁ − μ₂|
Bent (H₂O, SO₂):   μnet = 2 μbond cos(θ/2)
General — law of cosines:
μnet = √ μ12 + μ22 + 2μ1μ2 cos α
α = angle between the two bond dipole vectors
Lone Pair Correction
Lone pairs on the central atom contribute a dipole along the molecular axis.
When the central atom is more electronegative than terminals, LP dipoles add; otherwise subtract.
LPfactor = 1 + sign × 0.380 × LP × α
μfinal = μraw × LPfactor
Alignment factor α by geometry:
Trigonal pyramidal: α = 1.000  (LP on molecular axis)
Bent, 2 LP:   α = 0.139  (tetrahedral framework)
Bent, 1 LP:   α = 0.350  (trigonal planar framework)
Result capped at 11.0 D  (NaCl gas ≈ 9 D — physical maximum)
Worked Example — H₂O
O–H bond: ΔEN = 1.24
μbond = 1.121 × 1.240.924 = 1.121 × 1.220 = 1.368 D
θ = 104.5°  →  θ/2 = 52.25°  →  cos(52.25°) = 0.612
μraw = 2 × 1.368 × 0.612 = 1.675 D  (before lone pair correction)
Lone pair correction — O has 2 LP, bent geometry, alignment factor = 0.139:
LPfactor = 1 + (+1) × 0.380 × 2 × 0.139 = 1.106
μfinal = 1.675 × 1.106 = 1.851 D exp: 1.85 D ✓
§ 8

Molecular Mass Calculation

Uses IUPAC standard atomic weights — the weighted average over all naturally occurring isotopes — not integer mass numbers.

Molar Mass Sum
M = ∑i wi × ni g/mol
wi = standard atomic weight of element i
ni = count of atoms of element i in the molecule
Result rounded: round(M × 10 000) / 10 000
Example — glucose C₆H₁₂O₆:
M = 6 × 12.011 + 12 × 1.008 + 6 × 15.999
= 72.066 + 12.096 + 95.994
= 180.156 g/mol  ✓
§ 9

Isotopic Mass & Nuclear Stability

For each element the engine generates the full isotope range, calculates masses with binding energy corrections, and classifies stability from the neutron-to-proton ratio.

Isotope Range Generation
Nmin = max(0,  Z − 2)
Nmax = round(Z × 1.6) + 3
For Carbon (Z=6): Nmin=4, Nmax=13 → isotopes C-10 through C-19. Catches all known isotopes without a hardcoded database.
Isotopic Mass — Bethe–Weizsäcker Semi-Empirical Mass Formula
Let A = Z + N (mass number)  |  mp = 1.007 276 u  |  mn = 1.008 665 u
M(Z, A) = Z · mp + N · mn − B(Z, A)c² atomic mass in u
B(Z, A) = aVA − aSA2/3 − aC Z(Z−1)A1/3 − aA(A − 2Z)²A + δ(Z, N) MeV
Coefficients (AME 2020 fit): aV = 15.835  |  aS = 18.330  |  aC = 0.714  |  aA = 23.200  |  aP = 11.200
Pairing term: δ = +aP/√A (even-even)  |  −aP/√A (odd-odd)  |  0 (odd-A)
Conversion: 1 u = 931.494 MeV/c²  |  add Z × me (0.000 549 u) for atomic mass
Verification vs AME 2020 measured masses:
¹²C (Z=6, N=6): SEMF = 12.004 66 u  actual: 12.000 00 u  |err| = 0.005 u
¹⁶O (Z=8, N=8): SEMF = 15.998 91 u  actual: 15.994 92 u  |err| = 0.004 u
⁵⁶Fe (Z=26,N=30): SEMF = 55.935 20 u  actual: 55.934 94 u  |err| = 0.000 26 u ✓
²⁰⁸Pb (Z=82,N=126): SEMF = 207.979 u  actual: 207.977 u  |err| = 0.002 u ✓
²³⁸U (Z=92,N=146): SEMF = 238.035 u  actual: 238.051 u  |err| = 0.016 u
The SEMF is a macroscopic liquid-drop model — it does not capture individual shell effects (magic numbers). For light nuclei (A < 12) errors reach ~0.008 u. For mid-range and heavy nuclei it achieves < 0.003 u, orders of magnitude better than a simple proton+neutron sum with no binding energy correction (which gives errors up to 1.7 u for uranium).
Neutron-to-Proton Stability Criterion
ρ = NZ neutron-to-proton ratio
{
ρexpected = 1.0Z ≤ 20  (light elements)
ρexpected = ln(Z) × 0.3Z > 20  (heavier elements need more neutrons)
stable = |ρ − ρexpected| < 0.15   AND   Z ≤ 83
¹²C (Z=6, N=6): ρ = 1.0, expected = 1.0 → |0| < 0.15 → stable ✓.   ¹⁴C (Z=6, N=8): ρ = 1.333, expected = 1.0 → |0.333| > 0.15 → unstable → β⁻ decay ✓.
§ 10

Chemical Equation Balancing — Gaussian Elimination over ℚ

The balancer constructs a stoichiometry matrix from atom counts and solves it exactly over the rationals using Gaussian elimination with partial pivoting — no floating-point drift, guaranteed integer coefficients.

Stoichiometry Matrix
Rows = elements, Columns = species (reactants negative, products positive)
M · c = 0 find the null vector c
Example:   CH₄ + O₂ → CO₂ + H₂O
CH₄O₂CO₂H₂O −10+10 −400+2 0−2+2+1
Gaussian Elimination with Partial Pivoting
1
For each column, find the first row ≥ current with nonzero entry (pivot)
2
Swap pivot row to current position
3
Divide entire pivot row by pivot entry   (makes pivot = 1)
4
For every other row r:  row[r] = row[r] − row[r][col] × pivot_row
5
Read back-substituted solution from pivot columns
Coefficient Normalization
LCD = LCM of all denominators in the rational solution vector
ci = round(ni × LCDdi )
g = GCD(|c₁|, |c₂|, …)  →  ci,final = ci / g
Result:
CH₄ + 2 O₂ → CO₂ + 2 H₂O   ✓
§ 11

Reaction Enthalpy from Bond Energies

ΔH is estimated using average bond dissociation energies for bonds broken minus bonds formed. Cross-element bond energies use Pauling's geometric mean + electronegativity correction.

Hess's Law — Bond Energy Method
ΔHrxn = ∑ D(bonds broken) − ∑ D(bonds formed) kJ/mol
Bond Dissociation Energies — NIST Lookup + Pauling Fallback
Primary source: NIST Webbook / Luo "Comprehensive Handbook of Chemical Bond Energies" (2007).
The engine uses a 80+ entry lookup table covering all bonds students are likely to draw. For exotic or unknown A–B pairs not in the table, Pauling's formula provides the fallback:
D(A–B) = √Deff(A–A) × Deff(B–B) + 96.485 × (EN(A) − EN(B))² kJ/mol
96.485 kJ/mol = Faraday constant F (1 eV per mol of electrons) — not R×T
Deff values are back-calculated from NIST cross-bond entries, not raw homonuclear BDEs
Bond order multipliers (derived from NIST single→double→triple ratios):
single: × 1.000  |  double: × 1.900  |  triple: × 2.700
The previous purely-derived formula had RMS error of 28.5% across common bonds (C-H off by 15%, N≡N off by 27%). The NIST lookup table gives 0% error for all tabled bonds. Pauling fallback is only used for exotic pairs (e.g. Te–Se, Ti–N) not found in reference literature.
Selected NIST Reference Values (kJ/mol, 298 K)
C–H: 411
N–H: 386
O–H: 459
C=C: 614
C=O: 749
N=O: 607
C≡C: 839
C≡N: 891
N≡N: 945
C–F: 485
C–Cl: 328
O=O: 498
Worked Example — CH₄ Combustion
CH₄ + 2 O₂ → CO₂ + 2 H₂O
Bonds broken:
4 × C–H:   4 × 411 = 1644 kJ
2 × O=O:   2 × 498 =  996 kJ
Total broken: 2640 kJ
Bonds formed:
2 × C=O (in CO₂):   2 × 803 = 1606 kJ  — NIST value for C=O in CO₂
4 × O–H:   4 × 459 = 1836 kJ
Total formed: 3442 kJ
ΔH = 2640 − 3442 = −802 kJ/mol  (literature: −890 kJ/mol, ~10% deviation)
The ~10% deviation is inherent to the average bond energy method — it uses average bond enthalpies from many compounds, not the exact bond energies in these specific molecules. The literature value comes from calorimetry of the actual reaction. This is a known limitation of the Hess-bond method that affects all chemistry textbooks using this approach; it cannot be eliminated without per-molecule quantum calculations.
§ 12

Nuclear Decay — Energetics & Half-Life

Decay mode is selected from branching ratios. Decay energy is estimated using Geiger-Nuttall inspired relationships for alpha and empirical fits for beta. The half-life calculator uses the standard radioactive decay law.

Decay Mode Selection
rand ∈ U[0, 1) — cumulative probability selection:
{
alpha decayif rand < Pα
beta minuselif rand < Pα + Pβ⁻
beta plus / ECelif rand < Pα + Pβ⁻ + Pβ⁺
gammaelse
Decay Equations — Conservation Laws
(A, Z) → (A−4, Z−2) + ⁴He + Qα alpha
(A, Z) → (A, Z+1) + e⁻ + ν̄e + Qβ beta⁻
(A, Z) → (A, Z−1) + e⁺ + νe + Qβ beta⁺
(A, Z*) → (A, Z) + γ + Qγ gamma
Alpha Decay Energy — Geiger-Nuttall Inspired
Qα ≈ 28.3 − 0.026A + 0.000046A2 − 0.03Z + 4 MeV
Qα = max(2 MeV, formula)   — floor prevents unphysical negatives
Radioactive Decay Law
N(t) = N0 × (12) tt½
Equivalent exponential:
N(t) = N0 × e−λt   where   λ = ln 2t½
After n half-lives:   Nn = N02n
Example — Ra-226 (t½ = 1600 yr), after 3200 yr:
n = 3200 / 1600 = 2 half-lives  →  N = N₀ / 2² = 25% remaining ✓
Instability Estimate — Unknown Isotopes
instability = |ρ − ρexpected|
t½,est = 106 − instability × 20 × (108 if Z > 83) seconds
More unstable isotopes get exponentially shorter half-lives. Z > 83 penalty reflects that all trans-bismuth nuclides have fundamentally shorter lives due to Coulomb repulsion overwhelming strong-force binding.
§ 13

Reaction Kinetics — Arrhenius & Collision Theory

The kinetics simulator runs a live particle collision model. Each collision is tested against the Arrhenius probability to determine if it results in a reaction.

Activation Energy from Bond Energies
Ea = 0.4 × ∑ D(bonds broken in rate-limiting step) kJ/mol
0.4 = empirical estimate that ~40% of bond energy is needed to reach the transition state
Catalyst effect:
Ea,cat = Ea × 0.6  (40% reduction)
Arrhenius Equation
k = A · e−Ea / RT
A = pre-exponential (frequency) factor ≈ 10¹⁰ – 10¹³ s⁻¹
R = 8.314 J·mol⁻¹·K⁻¹  |  T = temperature (K)
Normalized collision probability P ∈ [0, 1]:
Prxn = e−Ea / RT A cancels in normalized form
Reverse reaction:
Ea,rev = Ea,fwd − ΔH  →  Prev = e−Ea,rev / RT
Particle Collision Detection
d(A, B) = √(xA−xB)² + (yA−yB)²
collision if   d ≤ 16 px   AND   random() < Prxn
Elastic collision response (1D along collision normal):
v′A = (mA−mB)vA + 2mBvB mA + mB
v′B = (mB−mA)vB + 2mAvA mA + mB
Maxwell-Boltzmann Speed Scaling
vparticle = vbase × √T298
At T=298K: normal speed.   At T=1000K: speeds ×1.83.   At T=100K: ×0.58.
§ 14

Ring Strain & Cyclic Molecule Analysis

Ring strain is calculated from the difference between the ideal geometry bond angle and the constrained angle forced by the ring size, using Baeyer strain theory.

Internal Polygon Angle
θring = (n − 2) × 180°n
cyclopropane (n=3): θ = 60°  |  cyclobutane (n=4): θ = 90°
cyclopentane (n=5): θ = 108°  |  cyclohexane (n=6): θ = 120°
Baeyer Strain Angle
θideal = 109.5°  (tetrahedral sp³ carbon)
β = |θideal − θring| 2
cyclopropane:   β = |109.5 − 60| / 2 = 24.75°   very high strain
cyclobutane:   β = |109.5 − 90| / 2 = 9.75°   moderate
cyclopentane:   β = |109.5 − 108| / 2 = 0.75°   low
cyclohexane:   β = |109.5 − 120| / 2 = 5.25°   minimal (chair)
§ 15

Boiling Point Estimation — Intermolecular Forces

Normal boiling point is estimated from three additive contributions — London dispersion, dipole–dipole interaction, and hydrogen bonding — each calibrated against experimental reference data.

London Dispersion Term
Tb,London = 11.397 × MW2/3 × (1 + 0.3364 × ln Natoms) K
Natoms ≤ 2 (diatomics): log factor = 0
Monoatomics: Tb = 1.205 × MW0.9  (calibrated to He = 4.2 K)
Dipole–Dipole Term
Tb,dipole = 20 × ΔENavg if asymmetric and no H-bonds
= 0 if geometry is linear/tetrahedral symmetric with no central lone pairs
Hydrogen Bonding Term
Tb,Hbond = 20.920 × nHbond0.8628 × ENHbond × (1 + 1.032Nheavy) K
Nheavy = non-hydrogen atom count  |  sub-linear exponent = diminishing returns per H-bond
Total & Special Cases
Tb = Tb,London + Tb,dipole + Tb,Hbond K, at 1 atm
Metallic:  Tb = 500 + 8 × MW
Ionic:  Tb = 1600 + 700 × (1 + 0.8 × ΔENavg) × MW0.25   ≤ 6000 K
Ionic compounds require breaking the entire crystal lattice (Born-Landé) rather than overcoming molecular IMFs — hence the separate formula. Calibrated to NaCl (1686 K), MgO (3873 K), KCl (1693 K), Al₂O₃ (3250 K).
§ 16

Melting Point Estimation

Melting point uses two paths based on molecular flexibility: compact molecules use a crystal-packing model; flexible molecules use a Tb-ratio with a rotatable-bond penalty.

Compact Molecules (no rotatable bonds)
Tm = 20.003 × √MW base crystal-packing model
Symmetry corrections:
Linear symmetric (CO₂-type):   Tm ×= 1.600
Tetrahedral symmetric (CH₄-type):   Tm ×= 1.132
H-bond correction:   Tm ×= (1 + 0.979 × nHbond + 0.129 × ndonors)
Flexible Molecules (rotatable bonds)
ratio₀ = 0.7305 (H-bonds present)  |  0.7500 (no H-bonds)
ratio = ratio₀ × e−0.239 × nrotatable
Tm = Tb × ratio
§ 17

Pressure Adjustments — Clausius-Clapeyron

When the user changes ambient pressure, both boiling and melting points shift. Boiling point uses Clausius-Clapeyron anchored by Trouton's rule. Melting point uses logarithmic pressure scaling.

Boiling Point at Pressure P
Trouton's rule: ΔSvap ≈ 85 J·mol⁻¹·K⁻¹ for most liquids
Tb(P) = Tb,1atm 1 − 0.0978 × ln P
0.0978 = RΔSvap = 8.31485  |  ln(1) = 0 → no change at 1 atm ✓
Melting Point at Pressure P
Tm(P) = Tm,1atm × (1 + s × ln P)
{
s = 0.005ionic  (small ΔV of fusion)
s = 0.010metallic
s = 0.020molecular  (largest ΔV)
§ 18

Critical Point — Modified Guldberg's Rule

Critical temperature and pressure are estimated from the normal boiling point and molecular properties. Used internally to classify supercritical states.

Critical Temperature
Guldberg's rule: Tc/Tb ≈ 1.5–1.6 for non-H-bonding; higher with H-bonds
hbondRatio = 0.08 × nHbond × (1 + 2 / Natoms)
Tc = Tb × (1.55 + hbondRatio) K
CO₂ (no H-bonds): Tc = 195 × 1.55 = 302 K  (actual: 304 K, err <1% ✓).   H₂O: hbondRatio = 0.08×2×(1+2/3) = 0.267 → Tc = 373 × 1.817 = 678 K  (actual: 647 K, ~5% over — inherent to the Guldberg approximation for strongly H-bonded systems).
Critical Pressure — Lydersen-Type
Pc = (0.113 + 0.0032 × Natoms)−2 × (1 + nHbond) atm
Monoatomics: Pc = 2 + 0.2 × MW
§ 19

Decomposition Temperature — Weakest Bond

Decomposition temperature is derived algorithmically from the weakest bond in the molecule using Pauling bond energies and a Boltzmann distribution factor — no lookup table required.

Formula
Dmin = min over all bonds of   D(A–B)  (kJ/mol)
Tdecomp = Dmin × 1000 R × 4.5 K
R = 8.314 J·mol⁻¹·K⁻¹  |  4.5 = empirical Boltzmann factor
H₂O: Dmin(O–H) = 459 kJ/mol → Tdecomp = 459 000 / (8.314 × 4.5) ≈ 12 280 K — correctly predicts water survives until plasma conditions. N–N single bond (163 kJ/mol) → Tdecomp ≈ 4 360 K.
§ 20

Ideal Gas Laws — PV = nRT

For gas-phase molecules, the engine applies PV = nRT using the current temperature and pressure to compute molar volume, density, and mean free path.

Molar Volume & Gas Density
Vm = RTP L/mol
R = 0.08206 L·atm·mol⁻¹·K⁻¹  |  at STP: Vm = 22.4 L/mol ✓
ρ = MWVm g/L
Mean Free Path
λ = kBT √2 · π · d² · P nm
kB = 1.381 × 10⁻²³ J/K  |  d ≈ 0.3 + 0.05 × Natoms nm  (molecular diameter proxy)
P in Pascal (= Patm × 101 325)
Verification — N₂ at STP (d ≈ 0.40 nm): λ ≈ 66 nm. Engine estimate: ~68 nm ✓.
§ 21

Molecular Structure Space — Combinatorial Derivation

The number of distinct molecular structures Atomency can represent is derived from three independent combinatorial factors. The result dwarfs PubChem's catalogue by over 36 orders of magnitude.

Three Factors for N atoms
(1) Element combinations — choosing N atoms from 118 elements with repetition:
E(N) = (117 + N)! 117! · N! multiset coefficient
(2) Structural topologies — labeled spanning trees (Cayley's formula):
T(N) = NN − 2 tree topologies; cyclic adds more
(3) Bond type arrangements — 5 types across N − 1 bonds:
B(N) = 5N − 1 single · double · triple · ionic · aromatic
Total Structure Count — Summed over N = 2 to 15
S = ∑N=215 E(N) × T(N) × B(N)
S ≈ 2.56 × 1044 structures
PubChem catalogue:   1.16 × 108 compounds
S PubChem ≈ 2.2 × 1036 × larger
Conservative lower bound — counts only acyclic tree structures from 2–15 atoms. Ring systems, larger molecules, and isotopic variants all add further to the space. The engine places no hard cap on atom count or element choice.
§ 22

Electron Configuration — Aufbau Principle & Shell Occupancy

The engine reads full electron configurations from the elements dataset — including all Aufbau exceptions for transition metals and lanthanides. It then parses the subshell notation algorithmically to derive shell occupancy, valence electrons, and hybridization without hardcoding per-element rules.

Subshell Capacity — Pauli Exclusion Principle
max(l) = 2(2l + 1) electrons per subshell
s  (l = 0):   2(1) = 2 electrons
p  (l = 1):   2(3) = 6 electrons
d  (l = 2):   2(5) = 10 electrons
f  (l = 3):   2(7) = 14 electrons
Aufbau Filling Order — Energy Diagonal Rule
Subshells fill in order of increasing (n + l), with lower n first for equal sums:
1s→ 2s→ 2p→ 3s→ 3p→ 4s→ 3d→ 4p→ 5s→ 4d→ 5p→ 6s→ 4f→ 5d→ 6p→ 7s→ 5f→ 6d
Notable Aufbau exceptions in the dataset (half-filled/fully-filled d stability):
Cr (Z=24):   [Ar] 4s¹ 3d⁵   not [Ar] 4s² 3d⁴  — half-filled d more stable
Cu (Z=29):   [Ar] 4s¹ 3d¹⁰   not [Ar] 4s² 3d⁹  — fully-filled d more stable
Au (Z=79):   [Xe] 6s¹ 4f¹⁴ 5d¹⁰  — same d-stability exception
Shell Occupancy — Algorithmic Derivation
The engine parses the subshell string and sums electrons per principal quantum number n:
Nshell n = ∑subshells with principal number n esubshell
Maximum capacity per shell (from Pauli):
Nmax,n = 2n²
n=1 (K): max 2
n=2 (L): max 8
n=3 (M): max 18
n=4 (N): max 32
Valence Electrons — Outermost Shell
V = Nshell nmax electrons in the highest occupied principal shell
For transition metals, valence includes both the outer s electrons and the partially filled d subshell.
Worked examples:
Carbon (Z=6):   1s² 2s² 2p²  →  shell 2 has 2+2 = 4 valence e⁻
Chlorine (Z=17):   …3s² 3p⁵  →  shell 3 has 2+5 = 7 valence e⁻
Iron (Z=26):   …4s² 3d⁶  →  n=4: 2e, n=3 d: 6e  →  8 d+s valence e⁻
Shell Diagram Display in the Lab
The side panel renders a compact shell diagram from the parsed occupancy:
Shell label
K   L   M   N
→
Electron count (valence highlighted gold)
2 · 8 · 8 · 1
← Potassium (K, Z=19)