Point Defect Model

Many intermetallic compounds are stable over a range of compositions because atoms can swap between sublattices or leave vacancies behind. In the CALPHAD community this is captured by the Compound Energy Formalism (CEF), where each sublattice carries its own site-fraction variables and the Gibbs energy is expanded over all “endmember” configurations. landau provides an equivalent but physically transparent formulation derived directly from the semi-grand canonical partition function, making it straightforward to connect first-principles point-defect calculations to thermodynamic phase diagrams.

The notebook first sets up what all models share, then works through four sublattice models in the same order each — equations, setup in landau, graphs:

1. the exact model, the reference for the other three;

2. the low-temperature expansion (LTE) vs exact;

3. the capped sublattice vs exact;

4. the block-exclusion sublattice vs exact.

Common setting

Crystal model

Consider a crystalline phase whose primitive unit cell is partitioned into \(L\) crystallographic sublattices. Sublattice \(l\) contains \(n_l\) sites per unit cell; the fraction of all atomic sites belonging to sublattice \(l\) is

\[\eta_l = \frac{n_l}{\sum_{l'} n_{l'}}.\]

For a binary A–B alloy each site on sublattice \(l\) can be in three states: occupied by the host species, occupied by an antisite, or vacant.

Excess semi-grand potential of a defect

The excess semi-grand potential of a single point defect relative to perfect host crystal with the same number of atoms is

\[[\phi_i^l] = [g_i^l] - \Delta\mu\,[n_i^l],\]

where \([g_i^l] = E_i^l - T S_i^l\) is the defect formation free energy, \([n_i^l]\) is the excess number of B atoms contributed by defect \(i\) on sublattice \(l\) at that fixed number of atoms, and \(\Delta\mu = \mu_B - \mu_A\) is the semi-grand chemical potential. An antisite exchanges one atom for one of the other species, \([n] = \pm 1\). A vacancy takes an atom off its site, so the remaining atoms are richer in the other species: in B2, a vacancy fraction \(x\) on the A sublattice gives \(c = 1/(2 - x) \approx \tfrac12 + \tfrac{x}{4} = \tfrac12 + \eta\,[n]\,x\) with \(\eta = \tfrac12\), so \([n] = +\tfrac12\) for an A-site vacancy and, likewise, \(-\tfrac12\) for a B-site vacancy.

Sublattice potential

Defects on different sublattices are independent, so the partition function factorises. With the per-site partition variables \(z_i^l = e^{-\beta[\phi_i^l]}\) and \(s_l = \sum_i z_i^l\), every model below gives the semi-grand potential per atom and the concentration in the form

\[\phi(\beta, \Delta\mu) = \phi_\text{host}(\beta, \Delta\mu) - \frac{1}{\beta} \sum_l \eta_l\, L(s_l), \tag{1}\]

\[c(\beta, \Delta\mu) = c_\text{host} + \sum_l \eta_l \sum_i [n_i^l]\, x_i^l, \qquad x_i^l = \frac{z_i^l}{D(s_l)}, \tag{2}\]

where \(x_i^l\) is the defect site fraction. Since \(c = -\partial\phi/\partial(\Delta\mu)\), the two functions are tied together by \(D = 1/L'\). The four models differ only in \(L\) and \(D\) (\(m\) is a maximum defect site fraction, \(0 < m < 1\); \(H(m) = -m\ln m - (1 - m)\ln(1 - m)\)):

Model	\(L(s)\)	\(D(s)\)	\(\sum_i x_i = s/D\)
exact	\(\ln(1 + s)\)	\(1 + s\)	\(s/(1 + s) < 1\)
low-temperature expansion	\(s\)	\(1\)	\(s\) (unbounded)
capped, \(s \le m/(1 - m)\)	\(\ln(1 + s)\)	\(1 + s\)	\(s/(1 + s) \le m\)
capped, \(s > m/(1 - m)\)	\(H(m) + m \ln s\)	\(s/m\)	\(m\)
block exclusion	\(m \ln(1 + s/m)\)	\(1 + s/m\)	\(s/(1 + s/m) < m\)

Python API

All classes live in landau.phases.pointdefects:

Class	Role
ConstantPointDefect	Single defect with constant formation energy \(E_i^l\) and entropy \(S_i^l\)
PointDefectSublattice	Exact sublattice (Section 1); also carries the block-exclusion model (Section 4)
LowTemperatureExpansionSublattice	Low-temperature expansion (Section 2)
CappedPointDefectSublattice	Capped sublattice (Section 3)
PointDefectedPhase	Full phase = host AbstractLinePhase + list of sublattices; site_fractions(T, dmu) returns every \(x_i^l\)

The three sublattice classes are siblings under the abstract base AbstractPointDefectSublattice, which holds the shared machinery; each implements \(L\) and \(D\) (as _site_log_partition and _site_normalization). The block-exclusion model needs no class of its own.

Prelude

from dataclasses import replace

from landau.phases import LinePhase, IdealSolution
from landau.phases.pointdefects import (
    ConstantPointDefect,
    PointDefectSublattice,
    LowTemperatureExpansionSublattice,
    CappedPointDefectSublattice,
    PointDefectedPhase,
)
from landau.features import Locus
from landau.calculate import calc_phase_diagram
from landau.plot import plot_phase_diagram, plot_1d_mu_phase_diagram

import numpy as np
from scipy.constants import Boltzmann, eV
kB = Boltzmann / eV

import matplotlib.pyplot as plt

Worked example: B2 intermetallic

We model a B2 structure — two interpenetrating simple-cubic sublattices with one site per unit cell each (\(\eta_\alpha = \eta_\beta = 0.5\)). At stoichiometry the \(\alpha\)-sublattice is occupied by A and the \(\beta\)-sublattice by B.

Defect	Sublattice	\(E_i\) (eV)	\([n_i]\)
B\(_\alpha\) antisite	\(\alpha\) (A sites)	0.28	+1
V\(_\alpha\) vacancy	\(\alpha\) (A sites)	0.45	+½
A\(_\beta\) antisite	\(\beta\) (B sites)	0.30	−1
V\(_\beta\) vacancy	\(\beta\) (B sites)	0.50	−½

The \(\alpha\)-sublattice defects are given lower formation energies than their \(\beta\) counterparts, so the compound is asymmetric about \(c = 0.5\), reaching further toward the B-rich side.

All four models use these defects, and the compound competes against the same terminal solid solutions and liquid. Terminal solid solutions bracket the B2 compound; a liquid phase appears at high temperature. Parameters are chosen so that the A terminal melts at ≈ 1160 K and the B terminal at ≈ 858 K, giving an asymmetric diagram where the B2 compound remains stable up to intermediate temperatures.

# Host: stoichiometric AB line phase
host = LinePhase('AB', fixed_concentration=0.5, line_energy=-0.40, line_entropy=1.0*kB)

# Defects on the alpha sublattice (A sites)
B_alpha = ConstantPointDefect('B_alpha', excess_energy=0.28, excess_entropy=0, excess_solutes=+1)
V_alpha = ConstantPointDefect('V_alpha', excess_energy=0.45, excess_entropy=0, excess_solutes=+0.5)

# Defects on the beta sublattice (B sites)
A_beta = ConstantPointDefect('A_beta', excess_energy=0.30, excess_entropy=0, excess_solutes=-1)
V_beta = ConstantPointDefect('V_beta', excess_energy=0.50, excess_entropy=0, excess_solutes=-0.5)

eta = 0.5  # sites per atom on each sublattice

# Terminal solid solutions
solid_a = LinePhase('A',    fixed_concentration=0, line_energy=-0.33, line_entropy=1.0*kB)
solid_b = LinePhase('B',    fixed_concentration=1, line_energy=-0.30, line_entropy=1.0*kB)
solid   = IdealSolution('solid', solid_a, solid_b)

# Liquid (higher entropy; T_melt ≈ 1160 K for A, ≈ 858 K for B)
liquid_a = LinePhase('A(l)', fixed_concentration=0, line_energy=-0.10, line_entropy=3.3*kB)
liquid_b = LinePhase('B(l)', fixed_concentration=1, line_energy=-0.13, line_entropy=3.3*kB)
liquid   = IdealSolution('liquid', liquid_a, liquid_b)

Comparison helpers

Sections 2–4 compare their model against the exact one with the same four graphs:

1. the total defect site fraction \(\sum_i x_i\) on the \(\alpha\) sublattice vs \(\Delta\mu\) at 1000 K;

2. \(\phi - \phi_\text{host}\) and \(c\) vs \(\Delta\mu\) at 400 K (dilute) and 2000 K (not dilute);

3. the T–c phase diagram;

4. the stable phase borders of both diagrams, overplotted.

Ts_diagram = np.linspace(200, 1400, 25)
mu_diagram = np.linspace(-0.8, 0.8, 200)


def defect_site_fraction(phase, T, dmu, sublattice='alpha'):
    """Total defect site fraction per site of a sublattice with `eta` sites per atom."""
    sub = next(s for s in phase.sublattices if s.name == sublattice)
    return sub.sublattice_fraction / eta * sum(phase.site_fractions(T, dmu)[sublattice].values())


def compare_site_fraction(model, label, T=1000, dmu=np.linspace(-0.6, 0.8, 700)):
    fig, ax = plt.subplots(figsize=(6, 4))
    ax.plot(dmu, defect_site_fraction(ab_exact, T, dmu), 'C0-', lw=4, alpha=0.4, label='exact')
    ax.plot(dmu, defect_site_fraction(model, T, dmu), 'C3--', label=label)
    ax.set_yscale('log')
    ax.set_ylim(1e-3, 3)
    ax.set_xlabel(r'$\Delta\mu$ (eV)')
    ax.set_ylabel(r'total defect site fraction $\sum_i x_i^\alpha$')
    ax.set_title(rf'$\alpha$ sublattice, T = {T} K')
    ax.legend()
    plt.tight_layout()
    return ax


def compare_potential(model, label, temps=(400, 2000), dmu=np.linspace(-0.6, 0.6, 600), raw=False):
    """phi - phi_host and c of the exact model and `model`; `raw` adds `model` without the clamp."""
    fig, axes = plt.subplots(2, len(temps), figsize=(10, 6), sharex=True)
    for j, T in enumerate(temps):
        phi_host = host.semigrand_potential(T, dmu)
        ax_phi, ax_c = axes[0, j], axes[1, j]
        ax_phi.plot(dmu, ab_exact.semigrand_potential(T, dmu) - phi_host, 'C0', label='exact')
        ax_c.plot(dmu, ab_exact.concentration(T, dmu), 'C0', label='exact')
        if raw:
            ax_phi.plot(dmu, sum(s.semigrand_potential_contribution(T, dmu) for s in model.sublattices),
                        'C3:', lw=1, alpha=0.7, label=f'{label} (raw)')
            ax_c.plot(dmu, host.line_concentration + sum(s.concentration_contribution(T, dmu) for s in model.sublattices),
                      'C3:', lw=1, alpha=0.7, label=f'{label} (raw)')
        ax_phi.plot(dmu, model.semigrand_potential(T, dmu) - phi_host, 'C3--', label=label)
        ax_c.plot(dmu, model.concentration(T, dmu), 'C3--', label=label)
        ax_phi.set_title(f'T = {T} K')
        ax_phi.set_ylim(-0.5, 0.03)
        ax_c.axhline(1, color='0.8', lw=0.8)
        ax_c.axhline(0, color='0.8', lw=0.8)
        ax_c.set_xlabel(r'$\Delta\mu$ (eV)')
        ax_c.set_ylim(-0.05, 1.6 if raw else 1.05)
        if j == 0:
            ax_phi.set_ylabel(r'$\phi - \phi_\mathrm{host}$ (eV/atom)')
            ax_phi.legend()
            ax_c.set_ylabel('concentration $c$')
            ax_c.legend(loc='upper left')
    fig.suptitle(f'Exact vs {label} point-defect model')
    plt.tight_layout()


def compare_borders(df_model, label):
    fig, ax = plt.subplots(figsize=(6, 5))
    for diagram, color, name, marker in [(df_exact, 'C0', 'exact', 'o'), (df_model, 'C3', label, 'x')]:
        borders = diagram[(diagram.locus != Locus.INTERIOR) & diagram.stable]
        ax.scatter(borders.c, borders['T'], s=16, c=color, marker=marker, alpha=0.7, linewidths=1, label=name)
    ax.set_xlim(0, 1)
    ax.set_ylim(Ts_diagram.min(), Ts_diagram.max())
    ax.set_xlabel('$c$')
    ax.set_ylabel('$T$ [K]')
    ax.legend()
    ax.set_title(f'Phase borders: exact vs {label}')

1. Exact model

Equations

A site of sublattice \(l\) is in exactly one state: occupied by its host species or by one of the defects \(i\). Summing over the states of a single site gives \(1 + s_l\), so \(L = \ln(1 + s)\) and \(D = 1 + s\):

\[\phi(\beta, \Delta\mu) = \phi_\text{host}(\beta, \Delta\mu) - \frac{1}{\beta} \sum_l \eta_l \ln\!\left(1 + \sum_i e^{-\beta[\phi_i^l]}\right),\]

\[x_i^l = \frac{e^{-\beta[\phi_i^l]}}{1 + \sum_{i'} e^{-\beta[\phi_{i'}^l]}}. \tag{3}\]

The site fraction is a softmax over all states of a site. Its denominator captures site competition: a site can only be in one state at a time, so \(\sum_i x_i^l \le 1\) and the concentration is always physically bounded.

Setup in landau

alpha_sublattice = PointDefectSublattice(
    name='alpha', sublattice=0, sublattice_fraction=eta, defects=[B_alpha, V_alpha],
)
beta_sublattice = PointDefectSublattice(
    name='beta', sublattice=1, sublattice_fraction=eta, defects=[A_beta, V_beta],
)
ab_exact = PointDefectedPhase(name='AB', line_phase=host, sublattices=[alpha_sublattice, beta_sublattice])

Graphs

Defect site fractions vs. composition

Equation (3) gives the site fraction \(x_i^l\) of each defect type as a function of \(\Delta\mu\). We scan \(\Delta\mu\) and plot the defect populations against the total B concentration.

T = 1000  # K
dmus = np.linspace(-1.5, 1.5, 500)

# Total B concentration at each chemical potential
c_total = ab_exact.concentration(T, dmus)

# Defect site fractions x_i^l (Eq. 3), keyed by sublattice and defect name
fracs = ab_exact.site_fractions(T, dmus)
fracs_alpha, fracs_beta = fracs['alpha'], fracs['beta']

fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)

linestyles = ['-', '--']
for ax, fracs, title in zip(
    axes,
    [fracs_alpha, fracs_beta],
    [r'$\alpha$ sublattice (A sites)', r'$\beta$ sublattice (B sites)'],
):
    for (name, frac), ls in zip(fracs.items(), linestyles):
        ax.plot(c_total, frac, ls, label=name)
    ax.set_xlabel('B concentration')
    ax.set_ylabel('Defect site fraction')
    ax.set_title(title)
    ax.legend()
    ax.set_xlim(0, 1)
    ax.set_yscale('logit')
    ax.set_ylim(1e-6, 1 - 1e-6)

fig.suptitle(f'Point defect site fractions — B2 model, T = {T} K')
plt.tight_layout()
plt.show()

On the B-rich side (\(c > 0.5\)), B\(_\alpha\) antisites dominate; on the A-rich side (\(c < 0.5\)), A\(_\beta\) antisites dominate. The site competition captured by the denominator in Eq. (3) ensures that defect fractions remain physically bounded even at large off-stoichiometry.

Defect site fractions vs. Δμ

The same data plotted as a function of the chemical potential Δμ shows the Fermi-like switching of defect populations directly in μ-space, without the non-linear mapping through composition.

fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)

for ax, fracs, title in zip(
    axes,
    [fracs_alpha, fracs_beta],
    [r'$\alpha$ sublattice (A sites)', r'$\beta$ sublattice (B sites)'],
):
    for (name, frac), ls in zip(fracs.items(), linestyles):
        ax.plot(dmus, frac, ls, label=name)
    ax.set_xlabel(r'$\Delta\mu$ (eV)')
    ax.set_ylabel('Defect site fraction')
    ax.set_title(title)
    ax.legend()
    ax.set_yscale('logit')
    ax.set_ylim(1e-6, 1 - 1e-6)

fig.suptitle(rf'Point defect site fractions vs. $\Delta\mu$ — B2 model, T = {T} K')
plt.tight_layout()
plt.show()

Isothermal 1D phase diagram

Before computing the full T–c diagram, we check phase stability at a single temperature (800 K, well below the terminal melting points). plot_1d_mu_phase_diagram draws the semi-grand potential φ(Δμ) of each phase along the chemical-potential axis — solid where the phase is stable (lowest φ), dashed where it is metastable — with the transitions marked.

T_iso = 800  # K
df_1d = calc_phase_diagram([solid, liquid, ab_exact], Ts=T_iso, mu=np.linspace(-1.0, 1.0, 200), keep_unstable=True)

fig, ax = plt.subplots(figsize=(8, 4))
plot_1d_mu_phase_diagram(df_1d, ax=ax)
ax.set_title(f'1D $\\mu$ phase diagram at T = {T_iso} K')
plt.show()

T–c phase diagram

df_exact = calc_phase_diagram([solid, liquid, ab_exact], Ts=Ts_diagram, mu=mu_diagram)
plot_phase_diagram(df_exact, triplepoints=True)

2. Low-temperature expansion vs exact

Equations

The exact \(\ln(1 + \sum_i z_i)\) couples all defects on a sublattice through site competition. When the defects are dilute (\(z_i = e^{-\beta[\phi_i^l]} \ll 1\), i.e. formation energies large compared to \(k_B T\)), the leading term of the low-temperature expansion \(\ln(1 + \sum_i z_i) \to \sum_i z_i\) drops that coupling and leaves each defect an independent excitation, \(L = s\) and \(D = 1\):

\[\phi(\beta, \Delta\mu) = \phi_\text{host} - \frac{1}{\beta} \sum_l \eta_l \sum_i e^{-\beta[\phi_i^l]}, \qquad x_i^l = e^{-\beta[\phi_i^l]}.\]

Because the site fractions are no longer normalised they can leave \([0, 1]\) outside the dilute regime, where the raw \(\phi\) would also run off to \(-\infty\). To stay physical and thermodynamically consistent, PointDefectedPhase clamps the total concentration to \([0, 1]\) and continues \(\phi\) as the tangent line past each saturation point — a line phase at the saturated concentration (\(c = 1\) with slope \(-1\), or \(c = 0\) with slope \(0\)). Because \(c = -\partial\phi/\partial\Delta\mu\), the clamped \(c\) and the linear \(\phi\) stay consistent, and the now-bounded \(\phi\) keeps an out-of-range LTE phase from spuriously dominating a phase diagram.

Setup in landau

LowTemperatureExpansionSublattice consumes the same ConstantPointDefect objects; only the sublattice class changes.

ab_lte = PointDefectedPhase(
    name='AB (LTE)',
    line_phase=host,
    sublattices=[
        LowTemperatureExpansionSublattice(
            name='alpha', sublattice=0, sublattice_fraction=eta, defects=[B_alpha, V_alpha]),
        LowTemperatureExpansionSublattice(
            name='beta', sublattice=1, sublattice_fraction=eta, defects=[A_beta, V_beta]),
    ],
)

Comparison

Total defect site fraction

compare_site_fraction(ab_lte, 'LTE');

The two agree while the defects are dilute. Above \(\sum_i x_i \approx 0.1\) the expansion keeps growing exponentially and runs through 1, where the exact model saturates.

Semi-grand potential and concentration

The dotted curves are the raw expansion, before the clamp.

compare_potential(ab_lte, 'LTE', raw=True)

At 400 K the expansion tracks the exact model near stoichiometry (\(|\Delta\mu| \lesssim 0.2\) eV) and parts from it where an antisite population grows large; the raw site fractions then overshoot \(c = 1\) (or \(0\)), and the clamp turns them into a line phase at the saturated concentration. At 2000 K the defects are not dilute even at stoichiometry and the two models differ everywhere. The clamped result keeps \(c = -\partial\phi/\partial\Delta\mu\) consistent and bounds \(\phi\), but it is still only a stand-in outside the dilute regime.

T–c phase diagram

df_lte = calc_phase_diagram([solid, liquid, ab_lte], Ts=Ts_diagram, mu=mu_diagram)
plot_phase_diagram(df_lte, triplepoints=True)

Phase borders

compare_borders(df_lte, 'LTE')

Both diagrams share the host, solid and liquid, so only the AB region can differ. The A-rich border nearly coincides; on the B-rich side the expansion overestimates the B\(_\alpha\) antisite population, and the compound reaches \(c = 0.75\) instead of \(0.65\) at the eutectic near 870 K. Use LowTemperatureExpansionSublattice when defect concentrations are small (the regime of most DFT-parametrised dilute-defect models) and the exact model otherwise.

3. Capped sublattice vs exact

Equations

A sublattice of \(N\) sites holds at most \(mN\) defects of all types together. A configuration with defects of types \(i_1, \dots, i_n\) on \(n\) distinct sites has weight \(\prod_k z_{i_k}\); summing over the \(\binom{N}{n}\) choices of sites and the types on each gives \(\binom{N}{n} s^n\), so the sublattice partition function is

\[Z_N(s) = \sum_{n=0}^{\lfloor mN \rfloor} \binom{N}{n} s^n .\]

Without the cap the sum runs to \(N\) and gives \((1 + s)^N\) — the exact model. In the thermodynamic limit the sum is dominated by its largest term, so per site

\[L(s) = \lim_{N \to \infty} \frac{1}{N} \ln Z_N(s) = \max_{0 \le x \le m} \left[ H(x) + x \ln s \right].\]

The unconstrained maximum \(x = s/(1 + s)\) is reachable while \(s \le s^* = m/(1 - m)\) and returns the exact \(\ln(1 + s)\); past \(s^*\) the maximum sits on the boundary \(x = m\), a line phase pinned at the cap:

\[\begin{split}L(s) = \begin{cases} \ln(1 + s) & s \le s^* \\ H(m) + m \ln s & s > s^* \end{cases} \qquad D(s) = \begin{cases} 1 + s & s \le s^* \\ s/m & s > s^* \end{cases}\end{split}\]

Above the cap the budget \(m\) is shared in proportion to the Boltzmann weights, \(x_i = z_i m / s\). \(L\) and \(L'\) are continuous at \(s^*\) (both branches give \(-\ln(1 - m)\) and \(1 - m\)) but \(L''\) jumps, so the concentration is continuous with a kink in \(\Delta\mu\), at \(\Delta\mu^* = ([g_i] - k_BT\ln\frac{1 - m}{m})/[n_i]\) for a single defect type. Being a restricted maximum, \(L\) never exceeds the exact \(\ln(1 + s)\): the capped \(\phi\) is never below the exact one, so the cap cannot make the phase spuriously win. The concentration stays bounded, so PointDefectedPhase has nothing to clamp.

Setup in landau

CappedPointDefectSublattice takes the same defects plus max_site_fraction \(= m\), here 0.1 on both sublattices.

m = 0.1
ab_capped = PointDefectedPhase(
    name='AB (capped)',
    line_phase=host,
    sublattices=[
        CappedPointDefectSublattice(
            name='alpha', sublattice=0, sublattice_fraction=eta, defects=[B_alpha, V_alpha],
            max_site_fraction=m),
        CappedPointDefectSublattice(
            name='beta', sublattice=1, sublattice_fraction=eta, defects=[A_beta, V_beta],
            max_site_fraction=m),
    ],
)

Comparison

Total defect site fraction

ax = compare_site_fraction(ab_capped, f'capped, $m = {m}$')
ax.axhline(m, color='0.6', lw=0.8, ls=':');

Identical to the exact model up to \(s^*\), then pinned at \(m\) — the kink at \(\Delta\mu^*\).

Semi-grand potential and concentration

compare_potential(ab_capped, 'capped')

Each sublattice holds at most \(\eta m = 0.05\) defects per atom, so \(c\) stays within \(0.5 \pm 0.05\) at any \(\Delta\mu\) and \(T\); \(\phi\) lies on or above the exact curve everywhere.

T–c phase diagram

df_capped = calc_phase_diagram([solid, liquid, ab_capped], Ts=Ts_diagram, mu=mu_diagram)
plot_phase_diagram(df_capped, triplepoints=True)

Phase borders

compare_borders(df_capped, 'capped')

Where the exact compound field stays inside \(0.5 \pm \eta m\) (the A-rich side, up to near its melting point) the borders coincide. On the B-rich side the cap pins the border at \(c = 0.55\), where the exact field reaches \(0.65\).

4. Block-exclusion sublattice vs exact

Equations

Partition the \(\eta N\) sites of a sublattice into blocks of \(k = 1/m\) sites. A block is either defect-free (weight 1) or holds exactly one defect of type \(i\) at one of its \(k\) positions (weight \(k z_i\)). Blocks are independent, so the per-block partition function is \(1 + k s\) and per site

\[L(s) = \tfrac{1}{k} \ln(1 + k s) = m \ln(1 + s/m), \qquad D(s) = 1 + s/m, \qquad x_i = \frac{z_i}{1 + s/m},\]

so \(\sum_i x_i = s/(1 + s/m) \to m\) as \(s \to \infty\). This is the maximum-entropy model for “each defect excludes \(1/m\) sites”: exact on the block lattice, mean-field for a real exclusion zone. It is a different constraint from the count cap of Section 3, not a smoothing of it.

Limits. \(m = 1\) recovers the exact model. For \(s \ll m\), \(L = s - s^2/(2m) + \dots\), so it agrees with the exact \(s - s^2/2 + \dots\) at first order and differs at second: the effective second virial coefficient is \(1/m\) times the exact one. The configurational entropy per site at defect fraction \(x\) is \(m H(x/m) + x \ln(1/m) = -x \ln x + x - x^2/(2m) + \dots\), the same ideal dilute entropy as the exact \(H(x) = -x \ln x + x - x^2/2 + \dots\); for \(m = 0.1\) it is 0.8% below at \(x = 0.01\), 5% at \(x = 0.04\) and 29% at \(x = m\).

Ordering. \(f(s) = \ln(1 + s) - m \ln(1 + s/m)\) has \(f(0) = 0\) and \(f' \ge 0\), so \(L \le L_\text{exact}\) for all \(s\): \(\phi\) is never below the exact one and the exclusion cannot make the phase spuriously win.

Cost. \(D\) has slope \(1/m\) from \(s = 0\), so the saturation is paid everywhere: at \(m = 0.1\) the total site fraction is 0.083 at \(s = 0.5\) and reaches \(m\) only asymptotically, and \(L\) is already 17% below exact at \(s = 0.05\). Use it when the physics is exclusion-like; use the capped sublattice when it is a hard count cap.

Setup in landau

No new class: a PointDefectSublattice on \(\eta m\) sites per atom (one per block), with every defect’s entropy raised by \(k_B \ln(1/m)\) (the \(k\) positions in a block), is exactly this model: \(\eta m \ln(1 + \sum_i z_i/m) = \eta\, m \ln(1 + s/m)\), and the concentration contribution \(\eta m \sum_i [n_i] (z_i/m)/(1 + s/m) = \eta \sum_i [n_i] z_i / (1 + s/m)\). The same \(m = 0.1\) as Section 3.

def blocked(defect):
    # k = 1/m positions per block
    return replace(defect, excess_entropy=defect.excess_entropy + kB * np.log(1 / m))


ab_blocked = PointDefectedPhase(
    name='AB (blocked)',
    line_phase=host,
    sublattices=[
        PointDefectSublattice(
            name='alpha', sublattice=0, sublattice_fraction=eta * m,
            defects=[blocked(B_alpha), blocked(V_alpha)]),
        PointDefectSublattice(
            name='beta', sublattice=1, sublattice_fraction=eta * m,
            defects=[blocked(A_beta), blocked(V_beta)]),
    ],
)

Comparison

Total defect site fraction

defect_site_fraction scales by sublattice_fraction / eta, which here converts the site fraction per block into one per original site.

ax = compare_site_fraction(ab_blocked, f'blocked, $m = {m}$')
ax.axhline(m, color='0.6', lw=0.8, ls=':');

Unlike the cap, the exclusion is felt before the saturation: the curve falls below the exact one from the dilute end on and approaches \(m\) only asymptotically.

Semi-grand potential and concentration

compare_potential(ab_blocked, 'blocked')

The concentration is bounded by the same \(0.5 \pm \eta m\) as the capped model, but reaches it smoothly, without a kink.

T–c phase diagram

df_blocked = calc_phase_diagram([solid, liquid, ab_blocked], Ts=Ts_diagram, mu=mu_diagram)
plot_phase_diagram(df_blocked, triplepoints=True)

Phase borders

compare_borders(df_blocked, 'blocked')

Because the exclusion lowers the defect populations everywhere, the compound field narrows on both sides — \(0.48 \le c \le 0.54\) against \(0.44 \le c \le 0.65\) for the exact model — including the A-rich side, which the cap left untouched.

Connection to the Compound Energy Formalism

The exact model is equivalent to the CEF at the level of ideal mixing, with the identifications:

• Endmember energies \(G_{s_1:s_2:\cdots}\) \(\leftrightarrow\) combinations of \(\phi_\text{host}\) and defect formation energies \(E_i^l\).

• Site fractions \(y_s^l\) \(\leftrightarrow\) defect site fractions \(x_i^l\) (Eq. 3).

The practical difference: the CEF parametrises endmember energies (which can be hard to access for unstable endmembers), while the point defect model parametrises individual defect formation energies — the natural output of DFT supercell calculations.

Parametrisation: endmember energies vs. point defect energies

The table below shows how to convert between the two conventions for a B2 phase where the \(\alpha\)-sites host A and the \(\beta\)-sites host B:

CEF endmember	Composition	PointDefectedPhase parameters
A:A (all A)	\(c=0\)	\(E_\text{host} + \eta_\beta\, E_{A_\beta}\)
A:B (stoich.)	\(c=0.5\)	\(E_\text{host}\)
B:B (all B)	\(c=1\)	\(E_\text{host} + \eta_\alpha\, E_{B_\alpha}\)

Here \(\eta_\beta E_{A_\beta}\) is the excess energy per atom when every \(\beta\)-site carries an A-antisite (\(\eta_\beta = 0.5\) \(\beta\)-sites per atom, each at cost \(E_{A_\beta}\)).

To verify: at \(T \to 0\) with \(\Delta\mu\) sufficiently negative to fully populate \(A_\beta\) antisites, the semi-grand potential of the defected phase should match the A:A endmember energy.

T_check = 10    # K — near 0 K, entropy term negligible
dmu_check = -0.35  # eV — negative enough that A_beta antisites are fully populated

# Semi-grand potential and concentration of defected phase in A-rich limit
phi_A_limit = ab_exact.semigrand_potential(T_check, dmu=dmu_check)
c_A_limit   = ab_exact.concentration(T_check, dmu=dmu_check)

# CEF A:A endmember energy from the table above (per atom)
E_host   = host.line_energy
E_A_beta = A_beta.excess_energy
eta_beta = beta_sublattice.sublattice_fraction
E_AA_endmember = E_host + eta_beta * E_A_beta  # = -0.40 + 0.5*0.30 = -0.25 eV

print(f'c at dmu={dmu_check}: {c_A_limit:.6f}  (expected ~0)')
print(f'phi(T={T_check} K, dmu={dmu_check}) = {phi_A_limit:.4f} eV')
print(f'E_AA endmember (CEF)               = {E_AA_endmember:.4f} eV')
print(f'Difference (entropy at {T_check} K)         = {phi_A_limit - E_AA_endmember:.4f} eV')

c at dmu=-0.35: 0.000000  (expected ~0)
phi(T=10 K, dmu=-0.35) = -0.2509 eV
E_AA endmember (CEF)               = -0.2500 eV
Difference (entropy at 10 K)         = -0.0009 eV

Further Reading

• Dissertation §11.1 (Poul 2026): Full derivation of the exact model (Section 1) from the semi-grand canonical partition function.

• Khatri, Koju & Mishin (2024) J. Phase Equilib. Diffus. 45, 375–393 — law-of-mass-action / Kröger–Vink approach; the semi-grand potential model above is equivalent at the ideal-mixing level.

• Westwood et al. (2015) Comput. Phys. Commun. 196, 145–151 — pycpd Python module implementing the dilute (low-temperature) limit.