Skip to content

Workflow internals

Overview

IsolatedImpurityWorkChain answers the question: "What is the smallest supercell whose boundary is far enough from the muon that the forces on boundary atoms are negligible?"

The algorithm iterates through increasingly large, nearly cubic supercells until two force-based convergence criteria are satisfied simultaneously.

Supercell generation

Nearly cubic supercells

The first supercell is generated by ScGenerators.initialize using pymatgen's CubicSupercellTransformation. This transformation finds a supercell matrix that minimises the deviation from a perfect cube for the available lattice, subject to:

  • min_atoms: number of unit cell atoms + 1 (to ensure at least one muon position)
  • max_atoms: 8 × unit cell atoms (caps the size at ~8× volume)
  • min_length: user-supplied or auto-set to min(lattice.abc) + 1 Å

Each subsequent iteration calls ScGenerators.re_initialize, which increments min_atoms and min_length by one unit relative to the previous supercell, guaranteeing strictly larger cells.

Muon site placement

The muon is inserted at the first Voronoi interstitial site of the unit cell via pymatgen.analysis.defects.generators.VoronoiInterstitialGenerator. A small offset of +0.001 (fractional) is added to each coordinate to break exact crystallographic symmetry and avoid potential issues with symmetry-constrained forces.

The fractional coordinate is then transformed into the supercell frame:

mu_frac_coord_SC = (mu_frac_coord_UC @ inv(sc_matrix)) % 1

The muon is always the last atom in the supercell.

Differential forces

The convergence check uses differential forces rather than raw forces with the muon:

Δf = f(supercell with muon) − f(supercell without muon)

The without-muon supercell is geometrically identical (same atoms in same positions), except the muon H atom is removed. The subtraction cancels any systematic DFT error common to both systems and isolates the muon-induced perturbation.

For the without-muon case a zero-force placeholder [0, 0, 0] is appended at the muon position index so that array lengths match for the subtraction.

Charge handling

When charge_supercell=True, the QE input for the muon-bearing cell includes tot_charge=1.0 (positive muon). The without-muon cell is run neutrally to prevent unphysical charge configurations in the reference calculation.

Convergence criteria

Both criteria must pass in the same iteration for the workflow to stop.

Criterion 1: minimum force

min( |Δf_i| for i ≠ muon, Δf_i ≠ 0 ) < conv_thr

If the smallest differential force on any host atom is already below the threshold, the cell is large enough in the minimum-force sense. This catches the common case where the muon perturbation falls off quickly (e.g. simple metals).

Criterion 2: exponential decay

For each chemical species s with at least mnasf (default 4) atoms:

  1. Compute distances d_i from the muon to each atom of species s.
  2. Filter to atoms with |Δf_i| < max_force_thr and |Δf_i| > 0.
  3. Fit |Δf_i| = A × exp(−τ × d_i) using scipy.optimize.curve_fit.
  4. Compute the minimum distance for convergence: d_conv = ln(A / conv_thr) / τ
  5. Pass if max(d_i) ≥ d_conv.

If the fitting fails (data not exponential) the species is marked as not converged. If the entire fitting step raises an unexpected exception, the workflow exits with code 704.

Context variables

Variable Meaning
ctx.n Current iteration counter (AiiDA Int, incremented by increment_n_by_one)
ctx.structure Effective unit cell (may differ from inputs.structure if pre-relaxation ran)
ctx.sup_struc_mu Current candidate supercell with muon
ctx.sup_struc_without_mu Current candidate supercell without muon
ctx.musite ArrayData holding the Voronoi fractional coordinate ('Voronoi_site')
ctx.sc_mat ArrayData holding the current SC transformation matrix ('SCmat')
ctx.traj_out Dict {'with_muon': TrajectoryData, 'without_muon': TrajectoryData} from the last pair of SCF calculations
ctx.conv_thr Cached convergence threshold (populated on first call to continue_iter)

How forces are extracted

After each pair of SCF calculations, inspect_run_get_forces reads the output_trajectory from each PwBaseWorkChain result and stores them in ctx.traj_out. The trajectory contains the forces (in eV/Å) as its only time-step. check_if_conv_achieved then reads:

traj_out.get_array("forces")[0]           # shape (N_atoms, 3)

For MLIP jobs via aiida-pythonjob, inspect_run_get_forces wraps the raw forces list output from the PythonJob into a minimal TrajectoryData-compatible store so that the same check_if_conv_achieved calcfunction can be reused unchanged.