FindMuonWorkchain: Automated DFT+μ Workflow

• FindMuonWorkchain is an automated workflow for first-principles muon site exploration using DFT+μ, incorporating candidate generation, optimization, and field evaluation.
• It integrates supercell convergence and symmetry analysis to filter and cluster candidate muon stopping sites with accurate dipolar and hyperfine field calculations.
• Leveraging Quantum ESPRESSO and MUESR, the workflow provides reproducible μSR predictions validated against benchmark experiments across diverse material systems.

FindMuonWorkchain is a top-level automated AiiDA workflow for performing first-principles muon site finding and interaction calculations using the DFT+$\mu$ framework, as described and implemented by Onuorah et al. (Onuorah et al., 2024). It operationalizes the strategy of modeling the muon as a hydrogen impurity in density functional theory and encompasses the generation, relaxation, and post-processing of candidate muon stopping sites in arbitrary host materials, automating the complex sequence required for quantitative muon spin spectroscopy calculations.

1. Workflow Architecture and Stepwise Execution

FindMuonWorkchain is implemented within the aiida-muon plugin and coordinates the end-to-end DFT+$\mu$ protocol. The high-level orchestration involves two major components: supercell convergence (IsolatedImpurityWorkChain, imported from aiida-impuritysupercellconv) and muon site search, relaxation, and interaction evaluation (in FindMuonWorkchain proper).

The workflow executes as a seven-step pipeline:

1. Initial Position Generation: $N_0$ candidate interstitial muon positions are generated in the primitive cell using a regular grid scheme with user-customizable spacing $d_u$ (default: $1\ \text{\AA}$). Points within $1\ \text{\AA}$ of any host atom are discarded, and space-group symmetry analysis removes redundancies, yielding distinct starting sites.
2. Supercell Convergence: If not specified, IsolatedImpurityWorkChain is invoked to determine the minimal supercell size that renders boundary forces from the inserted muon negligible, using a force-difference recipe with a typical threshold $A_F = 1\times10^{-3}$ Ry/Bohr ($\approx 0.0257$ eV/\AA).
3. Structure Preparation and Relaxation: For each candidate site, a supercell is constructed and a H pseudopotential is placed at the specified fractional coordinate. Full structural relaxation—using aiida-quantumespresso’s PwBaseWorkChain or PwRelaxWorkChain—is performed for each initial configuration.
4. Relaxation Filtering: If fewer than $60\%$ of relaxations converge, the protocol aborts.
5. Site Clustering and Symmetry Analysis: Relaxed structures are analyzed to identify symmetry-inequivalent muon sites based on distance ($A_d$, $\mu$0) and energy ($\mu$1, $\mu$2 eV) tolerances, with additional magnetic-symmetry expansion for magnetic hosts.
6. Spin-Polarized SCF and Contact Field Evaluation (Optional): Spin-polarized self-consistent runs are performed for each unique relaxed configuration, extracting the Fermi-contact field $\mu$3 from the computed spin density.
7. Dipolar Field Calculation: The classical dipolar field $\mu$4 at each site is computed via MUESR code, using a Lorentz-sphere summation of host nuclear or electronic moments.

All outputs—including relaxed supercell structures, total energies $\mu$5, energy differences $\mu$6, Fermi-contact fields $\mu$7, dipolar fields $\mu$8, and combined internal fields $\mu$9—are persisted in the AiiDA database (Onuorah et al., 2024).

2. Algorithms for Muon Site Identification and Supercell Convergence

Candidate muon sites are generated on a uniform grid: for lattice vectors $N_0$0 and grid spacing $N_0$1, the number of points along direction $N_0$2 is $N_0$3. Site exclusion proceeds by removing points within $N_0$4 of any host atom, and symmetry reduction is performed with space-group analysis, resulting in $N_0$5 inequivalent candidates.

Supercell convergence is ensured via the IsolatedImpurityWorkChain. The method iteratively increases the supercell size until the muon-induced force perturbation at the boundary (modeled using a H pseudopotential at a generic interstitial) drops below $N_0$6. At each iteration, the procedure:

• Inserts a H atom at a selected interstitial site.
• Calculates forces with and without the impurity by DFT-SCF.
• Fits the force difference $N_0$7 vs. distance to an exponential decay, $N_0$8.
• Accepts the supercell if the maximum remaining $N_0$9 and the maximal boundary distance exceeds $d_u$0.

If convergence is not achieved within $d_u$1, the workflow exits with an explicit error (Onuorah et al., 2024).

3. Computational Settings and Site Post-processing

The workchain employs Quantum ESPRESSO via aiida-quantumespresso for all DFT tasks. Pseudopotential selection defaults to SSSP PBE efficiency v1.3, with typical planewave cutoffs of 60 Ry (PW) and 480 Ry (charge). The $d_u$2-point mesh spacing defaults to $d_u$3 and smearing is Gaussian with $d_u$4 Ry. Relaxations proceed until all atomic forces are less than $d_u$5 Ry/Bohr; SCF energy convergence is set to $d_u$6 Ry.

Site clustering reflects a multi-stage protocol:

• Pairs of sites within $d_u$7 and $d_u$8 are reduced, keeping the lower-energy configuration.
• Pymatgen’s SpacegroupAnalyzer distinguishes crystal symmetry-equivalent sites within $d_u$9 ($1\ \text{\AA}$0) and $1\ \text{\AA}$1 ($1\ \text{\AA}$2 eV).
• For magnetic hosts, candidate sites are expanded over the magnetic supercell, with re-queueing of missing inequivalent relaxations.

4. Interaction Calculations: Dipolar and Hyperfine Fields

FindMuonWorkchain supports rigorous evaluation of both classical dipolar and electronic hyperfine fields at muon sites.

• Dipolar Field: For $1\ \text{\AA}$3 host magnetic moments $1\ \text{\AA}$4, the classical dipolar field at the muon site $1\ \text{\AA}$5 is:

$1\ \text{\AA}$6

where $1\ \text{\AA}$7 is the displacement vector from the muon to each $1\ \text{\AA}$8. Implementation utilizes MUESR for a Lorentz-sphere summation over all relevant moments.

• Fermi-Contact (Isotropic Hyperfine) Field: The Fermi-contact field is extracted as:

$1\ \text{\AA}$9

where $1\ \text{\AA}$0 is the spin density at the muon site.

• Anisotropic Hyperfine Tensor: Although the full tensor is not exhaustively implemented in the primary workflow, its DFT expression is:

$1\ \text{\AA}$1

This suggests potential extension of the workflow to include full tensor evaluation when required.

5. Input Schema, Software Infrastructure, and Example Invocation

Core inputs to FindMuonWorkchain are provided as Python data objects:

• structure (AiiDA StructureData, e.g., from CIF)
• magnetic_structure (optional, mCIF or AiiDA MagneticStructureData)
• du (float, grid spacing, default $1\ \text{\AA}$2)
• charged_supercell (Bool: charged or neutral muon)
• hubbard and custom_hubbard (Bool and Dict: DFT+U control and manual $1\ \text{\AA}$3 values)
• Code labels for Quantum ESPRESSO (pw.x) and MUESR

The workchain interfaces with three main AiiDA computational elements:

Plugin/WorkChain	Role	Key Inputs/Outputs
IsolatedImpurityWorkChain	Supercell convergence	structure, force threshold $1\ \text{\AA}$4 SC, matrix
PwBaseWorkChain / PwRelaxWC	DFT structural relaxation	SC structure, pseudo, cutoffs $1\ \text{\AA}$5 energies, forces
calcfunction muesr_dipolar	Dipolar field evaluation	relaxed structure, mag. config $1\ \text{\AA}$6 $1\ \text{\AA}$7

Editor's term: "SC" = supercell; "PW" = plane wave.

Typical user invocation can occur in Python: $A_d$2 or from the command line: $A_d$3 (Onuorah et al., 2024)

6. Validation and Performance Across Test Cases

Validation is performed on diverse material prototypes:

• LiF and bcc Fe: The force-difference protocol reproduces benchmark supercell choices ($1\ \text{\AA}$8 for LiF, $1\ \text{\AA}$9 for Fe), in line with hyperfine-convergence tests.
• CaF$A_F = 1\times10^{-3}$0 (charged vs neutral): For Mu$A_F = 1\times10^{-3}$1, the linear F–Mu$A_F = 1\times10^{-3}$2–F state (bond length $A_F = 1\times10^{-3}$3) emerges as lowest energy in a $A_F = 1\times10^{-3}$4 SC. Neutral Mu$A_F = 1\times10^{-3}$5 converges to the cubic void center in $A_F = 1\times10^{-3}$6 SC. Displacement patterns match prior DFT+$A_F = 1\times10^{-3}$7 without ZPM correction.
• La$A_F = 1\times10^{-3}$8NiO$A_F = 1\times10^{-3}$9 (AFM insulator): $\approx 0.0257$0 SC, $\approx 0.0257$1 yields 52 initial sites; $\approx 0.0257$2 eV. Four muon–O bound sites (two apical, two planar) are recovered. Dipolar fields at apical sites ($\approx 0.0257$3 mT, $\approx 0.0257$4 mT) closely match experiment ($\approx 0.0257$5 mT, $\approx 0.0257$6 mT); the contact term is negligible. Demonstrates necessity for DFT+U in gap formation and $\approx 0.0257$7 localization on Ni.
• AV$\approx 0.0257$8Sb$\approx 0.0257$9 Kagome metals: Identical $60\%$0 supercells for K, Rb, and Cs hosts. The site p1 (between A and Sb layers) best models the Kubo–Toyabe ZF-$60\%$1SR depolarization, as judged by second-moment $60\%$2 calculations.
• LaCoPO (FM metal): $60\%$3 SC, $60\%$4(Co) $60\%$5 eV, 20 initial positions yield four candidates. Site p1 ($60\%$6 from P) is lowest energy, $60\%$7 mT (rescaled to $60\%$8 mT via $60\%$9). Highlights partial self-consistency effects for itinerant magnets (Onuorah et al., 2024).

7. Significance and Context in Computational Muon Spectroscopy

FindMuonWorkchain enables automated, reproducible DFT+$A_d$0 workflows for the computational quantification of muon stopping sites and their local fields, systematically linking experiment to first principles. A key advance is the fully automated supercell convergence, symmetry-aware clustering, and robust post-processing of nuclear and electronic fields at muon sites. This protocol extends compatibility to a wide variety of hosts, including complex insulators, metals, magnetic and nonmagnetic materials, providing direct validation against experimental $A_d$1SR data. As such, it constitutes a core digital infrastructure for computational muon spin spectroscopy (Onuorah et al., 2024).