RT Lab toolkit (hydrogen-rtlab)

The Python half of the RT Lab workbench at https://craigmerry.com/rtlab/. This package is built from the simulation/ directory of the toolkit zip that page serves; the only additions are packaging metadata in pyproject.toml and this paragraph. No material is claimed anywhere in it.

A working research scaffold for the hydrogen-framework hypothesis. It contains an interactive HTML workbench and executable numerical tools. It does not contain a discovered room-temperature superconductor, a DFT spectrum for the proposed framework, or a validated synthesis recipe.

What runs

Component	Implemented	Execution requirements
Pairing explorer	Allen–Dynes with both correction factors, spectrum integration, isotope sensitivity, comparisons	Browser, no network dependencies
Atomic candidates	Published-coordinate reconstruction, substitution, H vacancy, isotope, strain, seeded Mg/Al disorder, CIF export	Python; ASE optional for importing CIF
DFT screening	QE input generation, prerequisite checks, relaxation → SCF → DOS/DFPT → force constants → phonon DOS	Locally installed Quantum ESPRESSO and supplied UPFs
Pairing solver	Linearized isotropic imaginary-axis Eliashberg eigenvalue solver, grid controls	NumPy + SciPy
Stability comparison	Balanced competing-phase enthalpy minimization	Compatible user-supplied reference energies
Retention explorer	Single-event Arrhenius sensitivity	Browser; barrier is an assumption
EPW / SSCHA / pathways	Documented integration and decision gates	External, material-specific calculations; see ADVANCED.md

The last row is not an automated or already executed atomistic simulation. Wannier projections, competing structures and reaction pathways cannot be inferred reliably from the concept drawing.

Quick numerical run

From this simulation directory:

python -m pip install -e .
python -m rtlab doctor
python -m rtlab analyze --modes examples/two-mode.json --eliashberg --matsubara 192 --out analysis.json

The example has two synthetic Einstein modes: 40 meV with λ=0.7 and 140 meV with λ=1.6. At μ*=0.10, the Allen–Dynes result is approximately 215.91 K. The finite-grid isotropic Eliashberg model gives approximately 243.22 K with 192 positive Matsubara frequencies. These numbers are numerical model outputs, not predictions for MgAlFeH6, SrNH4B6C6, or any synthesized material. Results produced during the build are in examples/two-mode-result-192.json and examples/two-mode-result-384.json.

Repeat with a larger grid:

OPENBLAS_NUM_THREADS=1 python -m rtlab analyze --modes examples/two-mode.json --eliashberg --matsubara 384 --out analysis-384.json

For these small dense eigenproblems, one BLAS thread is often faster than a large thread pool. This setting does not change the physics. Do not use agreement at two grids as evidence that a real-material model is complete.

Import analysis.json into the workbench's Calculation pipeline view. The web page never reports this as a DFT or experimental result.

Analyze a calculated spectrum

python -m rtlab analyze --spectrum spectrum.dat --frequency-unit meV --mu 0.10 --eliashberg --matsubara 192 --out analysis.json

Accepted units are meV, THz and cm-1. Use exactly two columns, positive monotonically increasing frequency and nonnegative α²F. # comments and the omega_mev,a2f header are supported. α²F uses the dimensionless convention with energy as the horizontal coordinate and λ=2∫α²F(E)/E dE. Changing the horizontal frequency unit under this convention does not multiply the α²F values. Confirm the originating program's convention before converting a file.

EPW outputs can contain multiple smearing columns. Extract one explicitly; the importer deliberately refuses to guess. Remove any zero-frequency row explicitly, and retain the original file and processing history. An imaginary phonon must be investigated, never converted to a positive energy by taking its absolute value. The numerical integral cannot compensate for missing spectral tails or a poorly sampled spectrum.

Structure workflow

python -m rtlab init --out candidate.json
python -m rtlab check candidate.json
python -m rtlab convert candidate.json --out candidate.cif
python -m rtlab variant candidate.json --replace Al:Mg --out parent.json
python -m rtlab variant candidate.json --deuterate --out isotope.json
python -m rtlab variant candidate.json --vacancy-index 3 --out vacancy.json
python -m rtlab variant candidate.json --strain 0.02 --out strain-plus-two-percent.json
python -m rtlab disorder candidate.json --repeats 2 2 2 --seed 7 --out disorder-7.json

--strain 0.02 means 2% linear scaling, not 2% volume scaling. A vacancy changes stoichiometry; compare appropriate chemical potentials, not just total energies. Disorder preserves the Mg/Al counts but represents only one random occupancy, not an optimized special quasirandom structure or a disorder average. Sample multiple seeds and cell sizes.

The builtin MgAlFeH6 cell is reconstructed from Xu et al.'s rounded conventional a≈6.24 Å and H coordinate x≈0.24. It uses primitive vectors (0,a/2,a/2), (a/2,0,a/2), (a/2,a/2,0). Conventional coordinates map to primitive fractional coordinates as (−x+y+z, x−y+z, x+y−z), modulo one. The nine-atom cell has volume 60.742656 ų. It is a starting geometry requiring relaxation, not the authors' full-precision structure or a reproduction of their Tc result.

The geometry check tests cell handedness, finite positions and overlaps against neighboring images. It does not establish bonding, symmetry, stability or correct stoichiometry for an intended material. The neighboring-image search assumes a reasonably reduced cell; reduce pathological skew cells before use. Large structures above 512 atoms skip the pair-distance search and require external geometry validation.

For CIF import:

python -m pip install -e '.[cif]'
python -m rtlab convert authoritative-structure.cif --out candidate.json

The original B–C framework sketch is conceptual. No authoritative SrNH4B6C6 structure file is bundled. Obtain the authors' data, verify the phase and atom ordering, and import it. Do not use the decorative diagram coordinates as a DFT model.

Prepare and run Quantum ESPRESSO

Install QE separately on a configured workstation or cluster. A compatible release should provide pw.x, ph.x, dos.x, q2r.x and matdyn.x on PATH. No QE executables or pseudopotentials are downloaded or bundled.

python -m rtlab settings --out settings.json

Edit settings.json and provide a local path and SHA256 for every species, for example:

"pseudos": {
  "H": {"path": "/absolute/path/to/H.upf", "sha256": "replace_with_actual_64_character_sha256"}
}

Add the other elements in the same way. sha256sum /absolute/path/to/H.upf computes a checksum. For D, use the same appropriate H pseudopotential with D's mass. Consistent exchange-correlation choices and pseudopotential validation remain the researcher's responsibility. Use compatible norm-conserving pseudopotentials if continuing with EPW. Check spin-polarized and competing magnetic states for Fe compounds; the generated starting inputs are nonmagnetic.

Starter values (80/640 Ry, 8³ k points, 2³ phonon q points, 0.02 Ry smearing) are not converged settings. In particular, 2³ is only an initial phonon grid. Use converged cutoffs appropriate to the supplied UPFs. The DOS stage uses the SCF mesh and is preliminary; production band and projected-DOS calculations require denser meshes and material-specific settings.

python -m rtlab prepare candidate.json --settings settings.json --job jobs/candidate
python -m rtlab run --job jobs/candidate --stage relax
python -m rtlab run --job jobs/candidate --stage scf
python -m rtlab run --job jobs/candidate --stage dos
python -m rtlab run --job jobs/candidate --stage phonon
python -m rtlab run --job jobs/candidate --stage q2r
python -m rtlab run --job jobs/candidate --stage phdos
python -m rtlab collect --job jobs/candidate

prepare refuses missing or checksum-mismatched UPFs and nonempty job directories. It writes SCF input only after relaxation reports converged final coordinates. The final cell must be in angstrom and positions in crystal coordinates; unsupported units cause an explicit error instead of a guessed conversion. QE pressure uses kbar; 1 bar = 0.0001 GPa = 0.001 kbar. The default stress convergence threshold is a numerical threshold, not evidence of exact 1-bar control.

Use MPI explicitly when appropriate:

python -m rtlab run --job jobs/candidate --stage phonon --launcher 'mpirun -np 8' --timeout 86400

This executes locally inside the chosen environment. It does not submit a scheduler job or purchase compute. On a cluster, place the command inside your institution's scheduler script with its required resource directives.

state.json records inputs, outputs and status. A nonzero exit code, missing completion marker or failed convergence is a failure. An unchanged completed stage can be resumed without rerunning. Changed completed inputs/outputs, source structure or settings require a fresh job directory. Preserve the full job directory, including QE scratch data, until dependent calculations finish. Interrupted jobs can report running if the Python process itself was killed; inspect the process and logs before retrying. The runner is for one owner at a time; do not launch simultaneous writers in the same job directory.

Competing phases

python -m rtlab phases examples/phases-toy.json --out comparison.json

The toy energies are synthetic and test reaction balancing only. Real inputs must provide composition counts and total enthalpy in eV for that exact composition unit for the candidate and every reference. Use consistent pressure, pseudopotentials, functional, magnetism and zero-point corrections. For pressure studies include PV in enthalpy. The solver minimizes a nonnegative combination of references subject to element conservation. It reports the candidate's energy per atom relative to that supplied mixture. Missing competitors can change the answer; the output explicitly does not claim a global convex hull. Finite-temperature free energies and open hydrogen chemical potentials need a separate treatment.

Model equations and limits

The numerical solver assumes an isotropic conventional superconductor, an infinite flat electronic band and a supplied positive electron–phonon spectrum. It linearizes the imaginary-axis gap equation at Tc. For positive Matsubara frequencies ω_n=(2n+1)πkBT:

Λ_j = Σ_i λ_i E_i² / [E_i² + (2π kBT j)²]
Z_n = 1 + [Λ_0 + 2 Σ_(j=1..n) Λ_j] / (2n+1)
K_nm = πkBT [Λ_|n−m| + Λ_(n+m+1) − 2μ* θ_n θ_m]
       / sqrt(ω_n Z_n ω_m Z_m)

The normal-state Z expression includes the analytic infinite-frequency tail. The pairing matrix uses N positive frequencies. θ_n=1 below the selected Coulomb energy cutoff and zero otherwise: a symmetric cutoff convention, not a claim to duplicate all EPW conventions. Tc is located where the largest eigenvalue is one. μ* belongs to its specified cutoff; changing cutoff without consistently renormalizing μ* changes the model. The solver rejects a grid that cannot cover the cutoff and reports a bound if Tc is below its resolved temperature range. It does not compute an anisotropic gap, a real-frequency spectral function, quantum nuclear motion, phase coherence or a critical current.

Verification

OPENBLAS_NUM_THREADS=1 python -m unittest discover -s tests -v

Tests check known integrals, isotope energy scaling, the numerical eigenvalue condition and Matsubara-grid sensitivity, atomic stoichiometry and geometry, fixed-composition disorder, balanced phase mixtures, and job prerequisites. QE, EPW and SSCHA were unavailable in the build environment and were not executed. See VALIDATION.json for the numerical run summary and explicit execution limits.

Sources

• Sun & Zhu (2024), B–C clathrate hydride predictions
• Xu et al. (2026), electronically modified hydrides and atomic coordinates
• Quantum ESPRESSO PW input specification
• Quantum ESPRESSO phonon input specification
• EPW theory
• EPW superconductivity reference workflow
• SSCHA and electron–phonon coupling tutorial

The sources motivate the scientific workflow. No published structure, predicted Tc or literature calculation is represented as a newly executed result here.