binding_energies.py

Go to the documentation of this file.

"""
binding_energies.py — Defect cluster binding energies for bcc Fe / EUROFER97.
 
All equations and tables cite:
  Ghoniem, N.M. (2026), "A Cluster Dynamics Model for Radiation Damage
  Evolution in Ferritic-Martensitic Steels" (Rate_Equations.pdf).
 
Conventions
-----------
- Vacancy cluster of size m (pure void):  c_{m,0}
- He-vacancy cluster:  c_{m,ℓ}  — m vacancies + ℓ He atoms
- Interstitial cluster of size n:  c_n
- All energies in eV; lengths in m unless stated.
"""
 
import numpy as np
 
_kB     = 8.617333262e-5    # eV K^-1
_J_eV   = 6.241509074e18    # J → eV
 
 
# ── Helper ────────────────────────────────────────────────────────────────────
 
def atomic_radius(Omega):
    """Wigner-Seitz radius r_0 = (3Ω/4π)^{1/3} [m]."""
    return (3.0 * Omega / (4.0 * np.pi)) ** (1.0 / 3.0)
 
 
def capture_radius(n, Omega):
    """r_n = r_0 · n^{1/3} [m].  Eq. 86 (equivalent form)."""
    return atomic_radius(Omega) * np.asarray(n, dtype=float) ** (1.0 / 3.0)
 
 
# ── 1.  Vacancy cluster (void) binding energies ───────────────────────────────
 
# Atomistic fitting amplitudes A(m) for the void capillary correction.
# Table 18 of Rate_Equations.pdf.  Index = He content ℓ = 0..4 (pure void case).
_A_void = {0: 1.2353, 1: 2.9064, 2: 3.4147, 3: 2.1504, 4: -0.1590}  # eV
_lambda_void = 0.5756   # decay constant [vac^-1], = ln(100)/8
 
 
def E_b_void(m, E_f_v, gamma_s, Omega, lambda_void=None, A_void_0=None):
    """
    Binding energy of the m-th vacancy in a spherical void cluster.
 
    Capillary formula with atomistic correction (Eqs. 66-67):
      E_b(m) = E_f_v − A_cap · [m^{2/3} − (m−1)^{2/3}]
               + A_void(0) · exp(−λ·(m−1))     [Eq. 67 correction]
 
    where  A_cap = 4π·γ_s·r_0²  [eV],  r_0 = (3Ω/4π)^{1/3}.
 
    Special case m=1: return E_f_v (monomer formation energy).
 
    Parameters
    ----------
    m           : int or array
    E_f_v       : float [eV]
    gamma_s     : float [J/m²]
    Omega       : float [m³]
    lambda_void : float or None [vac^-1]
        Blend decay constant λ.  ``None`` (default) uses the module constant
        ``_lambda_void`` = ln(100)/8, preserving legacy behaviour.  Supplied
        from the ``dissociation`` sheet key ``lambda`` by ReactionRates so the
        workbook value is actually honoured — before this argument existed the
        workbook carried ``lambda`` and ``A_void_0`` but nothing read them.
    A_void_0    : float or None [eV]
        Atomistic correction amplitude A(0).  ``None`` uses ``_A_void[0]``.
 
    Returns
    -------
    E_b : float or array [eV]
    """
    r0     = atomic_radius(Omega)
    A_cap  = 4.0 * np.pi * gamma_s * r0**2 * _J_eV   # eV  — Eq. 66
    A_atm  = _A_void[0]      if A_void_0    is None else float(A_void_0)
    _lam   = _lambda_void    if lambda_void is None else float(lambda_void)
 
    m = np.asarray(m, dtype=float)
    capillary = A_cap * (m**(2.0/3.0) - (m - 1.0)**(2.0/3.0))
    correction = A_atm * np.exp(-_lam * np.maximum(m - 1.0, 0.0))
 
    Eb = np.where(m <= 1.0,
                  E_f_v,
                  E_f_v - capillary + correction)
    return float(Eb) if Eb.ndim == 0 else Eb
 
 
# ── 2.  He–vacancy (bubble) binding energies ─────────────────────────────────
 
# He virial EOS coefficients (Table 8, Eqs. 64-65)
_B2 = 1.67e-29    # m^3/atom  — 2nd virial coefficient
_B3 = 1.84e-58    # m^6/atom^2 — 3rd virial coefficient
 
 
def He_pressure(m, ell, Omega, T):
    """
    He pressure in a bubble with m vacancies and ℓ He atoms.
    Virial EOS (Eq. 64):
      P·V = N_He·k_B·T · [1 + B2·(N_He/V) + B3·(N_He/V)²]
 
    V = m·Ω  (bubble volume).
 
    Returns
    -------
    P : float [Pa]
    """
    if ell <= 0:
        return 0.0
    V  = float(m) * Omega
    rho = float(ell) / V        # He number density [m^-3]
    kB_SI = _kB / _J_eV         # eV/K → J/K
    P = float(ell) * kB_SI * T / V * (1.0 + _B2 * rho + _B3 * rho**2)
    return P
 
 
def E_b_bubble(m, ell, E_f_v, gamma_s, Omega, T):
    """
    Effective vacancy binding energy for a He-vacancy bubble (m, ℓ).
 
    Capillary model with He-pressure correction (Eqs. 70-73):
      E_b_bub(m, ℓ) = E_f_v − A_cap·[m^{2/3}−(m−1)^{2/3}]
                       + P_He(m,ℓ)·Ω  [He pressure work term]
 
    The pressure term is POSITIVE: removing a vacancy from a gas-filled
    cavity compresses the gas, costing work +P·Ω, so He raises the
    vacancy binding energy and stabilizes bubbles against thermal
    dissociation (equivalently c_v^eq(surface) ∝ exp[(2γ/r − P)Ω/kT]).
    This term creates the stable-bubble branch and critical radius r*
    of critical-bubble theory (Mansur & Coghlan).
 
    Parameters
    ----------
    m       : int
    ell     : int or float  (can be fractional mean loading)
    E_f_v   : float [eV]
    gamma_s : float [J/m²]
    Omega   : float [m³]
    T       : float [K]
 
    Returns
    -------
    E_b : float [eV]
    """
    r0    = atomic_radius(Omega)
    A_cap = 4.0 * np.pi * gamma_s * r0**2 * _J_eV
    m_f   = float(m)
    capillary = A_cap * (m_f**(2.0/3.0) - max(m_f - 1.0, 0.0)**(2.0/3.0))
    P     = He_pressure(m, ell, Omega, T)
    dE_P  = P * Omega * _J_eV   # Pa·m³ → eV
    return E_f_v - capillary + dE_P
 
 
# ── 3.  He binding energy to bubbles ─────────────────────────────────────────
 
# Atomistic fitting amplitudes A^He(α) for ℓ ≤ α.
# Table 19 of Rate_Equations.pdf.
_A_He = {2: 0.55, 3: 0.40, 4: 0.75}   # eV
_mu_He = 0.658                          # decay constant [He^-1]
 
def E_b_He(m, ell, E_s_He, Omega, T, B2=_B2, B3=_B3):
    """
    He binding energy: energy to remove one He from bubble (m, ℓ).
 
    Combined continuum + atomistic formula (Eqs. 76-77, Table 19):
 
      E_b^He(m, ℓ) = E_s^He + P_He(m,ℓ)·Ω   [continuum, Eq. 76]
                     + f_blend · A^He(α) · exp(−μ·(ℓ−α))  [atomistic, Eq. 77]
 
    where α = ℓ_max(m) = floor(alpha_He · m^{2/3}) is the trap-mutation limit,
    and f_blend fades the atomistic correction for large ℓ.
 
    Parameters
    ----------
    m      : int
    ell    : int  — current He content
    E_s_He : float [eV] — He interstitial solution energy (= 2.35 eV, Table 5)
    Omega  : float [m³]
    T      : float [K]
 
    Returns
    -------
    E_b : float [eV]
    """
    P = He_pressure(m, ell, Omega, T)
    E_cont = E_s_He + P * Omega * _J_eV   # continuum term
 
    # atomistic correction for small ℓ
    alpha = int(round(ell))   # closest integer loading
    A_atm = _A_He.get(alpha, 0.0)
    if A_atm != 0.0:
        E_atm = A_atm * np.exp(-_mu_He * max(ell - alpha, 0.0))
    else:
        E_atm = 0.0
 
    return E_cont + E_atm
 
 
def E_b_He_first(E_b_hV_1=2.30):
    """He binding energy to a monovacancy (first He, Table 5). Eq. 76."""
    return E_b_hV_1
 
 
def E_b_He_second(E_b_hV_2=2.00):
    """He binding energy to a monovacancy (second He, Table 5). Eq. 76."""
    return E_b_hV_2
 
 
# ── 4.  Interstitial loop binding energies ────────────────────────────────────
 
def _E_b_loop_i_continuum(n, E_f_i=3.77, G_shear=82e9, b=2.49e-10,
                           nu=0.29, gamma_sf=0.6, Omega=1.18e-29):
    """
    Continuum binding energy of the n-th SIA to an interstitial loop.
 
    E_b^{i,cont}(n) = E_f^i − μ_loop^{(i)}(n)   [Eqs. i_binding, chem_pot_i_loop]
 
    where μ = Gb·Ω/(4π(1−ν)R)·[ln(R/r_c)+1] + γ_sf·Ω/b,
    R = sqrt(n·Ω/(π·b)), r_c = b/3.
    """
    n = np.asarray(n, dtype=float)
    n_safe = np.maximum(n, 2.0)
    r_c    = b / 3.0
    R      = np.sqrt(n_safe * Omega / (np.pi * b))
    ln_fac = np.log(R / r_c) + 1.0
    mu     = (G_shear * b * Omega / (4.0 * np.pi * (1.0 - nu) * R)) * ln_fac \
             + gamma_sf * Omega / b
    mu_eV  = mu * _J_eV
    Eb     = E_f_i - mu_eV
    return Eb
 
 
def A_111_from_E_b_i2(E_b_i2, B_111):
    """
    Loop-binding amplitude ``A_111`` implied by a measured **di-interstitial
    binding energy** ``E_b^i(2)``.
 
    The small-n branch of :func:`E_b_loop_i` is the DFT power-law fit
    ``E_b^fit(n) = A_111 · n^{+B_111}`` (note the exponent is POSITIVE — binding
    *increases* with loop size).  Evaluated at n = 2 this is the binding of the
    second SIA to a mono-interstitial, i.e. the di-interstitial binding energy.
    Inverting,
 
        A_111 = E_b^i(2) · 2^{−B_111}
 
    Anchoring on ``E_b^i(2)`` rather than on ``A_111`` directly is preferable
    for calibration: ``A_111`` has no independent measurement, whereas
    ``E_b^i(2)`` is reported by DFT for bcc Fe (≈0.8 eV).
 
    Note that :func:`E_b_loop_i` returns the *blended* value at n = 2, which
    differs from ``E_b_i2`` by the continuum admixture (1 − w) ≈ 1e-4 at the
    default blend (n_tr = 25, σ_tr = 5) — negligible, but not identically zero.
 
    Parameters
    ----------
    E_b_i2 : float [eV]  — di-interstitial binding energy
    B_111  : float       — power-law exponent
 
    Returns
    -------
    A_111 : float [eV]
    """
    return float(E_b_i2) * 2.0 ** (-float(B_111))
 
 
def E_b_loop_i(n, A_111=0.7501, B_111=0.3873, A_100=0.7160, B_100=0.3581,
               n_tr=25.0, sigma_tr=5.0,
               E_f_i=3.77, G_shear=82e9, b_111=2.49e-10, nu=0.29,
               gamma_sf=0.6, Omega=1.18e-29):
    """
    Binding energy of the n-th SIA to an interstitial loop of size n−1.
 
    Small-n DFT fit blended to continuum limit (Eqs. Eb_smalln_fit, Eb_blended):
      E_b^fit(n) = A · n^{+B}                         [Eq. Eb_smalln_fit]
      E_b^cont(n) = E_f^i − μ_loop(n)                 [continuum elasticity]
      w(n)   = 0.5·(1 + tanh((n_tr − n)/σ_tr))        [Eq. weight_function]
      E_b(n) = w·E_b^fit + (1−w)·E_b^cont             [Eq. Eb_blended]
 
    For small n, w→1 and E_b follows the DFT power-law fit.
    For large n, w→0 and E_b approaches E_f^i (formation energy).
 
    Parameters
    ----------
    n            : int or array  (≥ 2)
    A_111, B_111 : float  — ½⟨111⟩ loop fit (Table 18)
    A_100, B_100 : float  — ⟨100⟩ loop fit  (Table 18)
    n_tr         : float  — blend center (25 SIAs)
    sigma_tr     : float  — blend width  (5 SIAs)
    E_f_i        : float  — SIA formation energy [eV]
    G_shear      : float  — shear modulus [Pa]
    b_111        : float  — ½⟨111⟩ Burgers vector [m]
    nu           : float  — Poisson ratio
    gamma_sf     : float  — stacking fault energy [J/m²]
    Omega        : float  — atomic volume [m³]
 
    Returns
    -------
    E_b : float or array [eV]
    """
    n = np.asarray(n, dtype=float)
    n_safe = np.maximum(n, 2.0)
 
    # Small-n DFT power-law fit (exponent is POSITIVE: E_b increases with n)
    Eb_fit = A_111 * n_safe**(B_111)
 
    # Large-n continuum elasticity
    Eb_cont = _E_b_loop_i_continuum(n_safe, E_f_i, G_shear, b_111,
                                     nu, gamma_sf, Omega)
 
    # Smooth blend: w→1 for small n (DFT regime), w→0 for large n (continuum)
    w  = 0.5 * (1.0 + np.tanh((n_tr - n_safe) / sigma_tr))
    Eb = w * Eb_fit + (1.0 - w) * Eb_cont
 
    return float(Eb) if Eb.ndim == 0 else Eb
 
 
def E_b_loop_100(n, A_100=0.7160, B_100=0.3581,
                 n_tr=25.0, sigma_tr=5.0,
                 E_f_i=3.77, G_shear=82e9, b_100=2.8665e-10, nu=0.29,
                 gamma_sf=0.6, Omega=1.18e-29):
    """
    Binding energy of the n-th SIA to a ⟨100⟩ interstitial loop of size n−1.
 
    Sibling of :func:`E_b_loop_i` for the sessile ⟨100⟩ population
    (½⟨111⟩→⟨100⟩ loop-conversion work).  Same small-n DFT fit blended to the
    continuum limit (Eqs. Eb_smalln_fit, Eb_blended), but with the ⟨100⟩ fit
    amplitudes ``A_100, B_100`` (Table 18 — already carried, dormant, by
    :func:`E_b_loop_i`) and the ⟨100⟩ Burgers vector ``b_100 = a`` in the
    continuum tail.
 
      E_b^fit(n)  = A_100 · n^{+B_100}
      E_b^cont(n) = E_f^i − μ_loop(n; b_100)
      w(n)        = 0.5·(1 + tanh((n_tr − n)/σ_tr))
      E_b(n)      = w·E_b^fit + (1−w)·E_b^cont
 
    ⟨100⟩ loops are slightly less bound than ½⟨111⟩ at small n
    (``A_100 < A_111``); emission is therefore slow but nonzero, feeding the
    ⟨100⟩ ``DISSOCIATION`` edge.
 
    Parameters mirror :func:`E_b_loop_i`; ``b_100`` defaults to the bcc-Fe
    lattice parameter ``a = 2.8665 Å`` (= ⟨100⟩ Burgers vector magnitude).
    """
    n = np.asarray(n, dtype=float)
    n_safe = np.maximum(n, 2.0)
    Eb_fit = A_100 * n_safe**(B_100)
    Eb_cont = _E_b_loop_i_continuum(n_safe, E_f_i, G_shear, b_100,
                                     nu, gamma_sf, Omega)
    w  = 0.5 * (1.0 + np.tanh((n_tr - n_safe) / sigma_tr))
    Eb = w * Eb_fit + (1.0 - w) * Eb_cont
    return float(Eb) if Eb.ndim == 0 else Eb
 
 
def E_b_loop_v(m, E_f_v, gamma_sf, Omega, b_111):
    """
    Vacancy loop binding energy (Eqs. 104-105).
 
    Continuum (Frank loop) with DFT anchor at m=2:
      E_b^vloop(m) = E_f_v − A_sf · [m^{1/2} − (m−1)^{1/2}]
 
    where A_sf = 2·b_111 · sqrt(π·γ_sf·Ω/√3) · _J_eV.
 
    Parameters
    ----------
    m       : int or array (≥ 1)
    E_f_v   : float [eV]
    gamma_sf: float [J/m²]  — stacking fault energy
    Omega   : float [m³]
    b_111   : float [m]     — ½⟨111⟩ Burgers vector
 
    Returns
    -------
    E_b : float or array [eV]
    """
    A_sf = 2.0 * b_111 * np.sqrt(np.pi * gamma_sf * Omega / np.sqrt(3.0)) * _J_eV
    m = np.asarray(m, dtype=float)
    Eb = np.where(
        m <= 1.0,
        E_f_v,
        E_f_v - A_sf * (m**(1.0/2.0) - (m - 1.0)**(1.0/2.0))
    )
    return float(Eb) if Eb.ndim == 0 else Eb
 
 
# ── 5.  Trap mutation barrier ─────────────────────────────────────────────────
 
# Table 27 of Rate_Equations.pdf
# Key: (m, ℓ)  — (vacancies, He atoms)
_E_TM_table = {
    (1, 5): 1.00,   # He5V1 — Morishita 2003
    (1, 6): 0.50,   # He6V1 — Morishita 2003
    (1, 7): 0.00,   # He7V1 — spontaneous (Gao 2011)
    (0, 4): 0.30,   # He4   self-trapping (Gao 2011)
    (0, 5): 0.10,   # He5   self-trapping (Gao 2011)
    (0, 6): 0.00,   # He6   spontaneous (Gao 2011)
}
 
 
def E_TM(m, ell, default=1e10):
    """
    Trap mutation barrier (Table 27, Eq. 83).
 
    Returns the barrier E_TM(m, ℓ) [eV].
    If (m, ℓ) not in the table, returns `default` (effectively disables TM).
 
    Parameters
    ----------
    m   : int  — vacancies
    ell : int  — He atoms
    default : float [eV]  — barrier for unknown (m, ℓ)
    """
    return _E_TM_table.get((int(m), int(ell)), default)
 
 
def Gamma_TM(m, ell, T, nu0_TM=1e12):
    """
    Trap mutation rate from bubble (m, ℓ) → (m+1, ℓ-1) [s^-1].
 
    Γ_TM(m, ℓ) = ν_0_TM · exp(−E_TM(m,ℓ) / k_B T)   (Eq. 142)
 
    Parameters
    ----------
    m, ell : int
    T      : float [K]
    nu0_TM : float [s^-1]  — attempt frequency (Table 27)
    """
    E = E_TM(m, ell)
    if E >= 1e9:
        return 0.0
    return nu0_TM * np.exp(-E / (_kB * T))
 
 
# ── 6.  Radiation re-solution rate ───────────────────────────────────────────
 
def Gamma_res(ell, phi_dot, b0):
    """
    Radiation re-solution rate for He from bubble (m, ℓ) [s^-1].
 
    Γ_res(m, ℓ) = b_0 · ℓ · φ̇   (Eq. 143)
 
    Parameters
    ----------
    ell     : int or float  — He content
    phi_dot : float [dpa/s] — displacement rate
    b0      : float [dpa^-1 He^-1]  — re-solution parameter (Table 29)
    """
    return b0 * float(ell) * phi_dot
 
 
# ── 7.  Maximum He loading per vacancy cluster ────────────────────────────────
 
def ell_max(m, alpha_He=1.7):
    """
    Maximum He loading for cluster with m vacancies.
    Trap-mutation limit: ℓ_max = floor(alpha_He · m^{2/3}).  Table 29.
    """
    return int(np.floor(alpha_He * float(m)**(2.0 / 3.0)))