diff --git a/AGENT.md b/AGENT.md
index c56d452..1d41e62 100644
--- a/AGENT.md
+++ b/AGENT.md
@@ -39,7 +39,7 @@ C(T) = (V · 2π/ℏ) · g · W² · Σᵢⱼ pᵢ |⟨χᵢ|Q-Q₀|χⱼ⟩|²
Where:
- `V`: supercell volume (cm³)
- `g`: degeneracy factor (spin/orbital)
-- `W`: electron-phonon coupling matrix element (eV)
+- `W`: electron-phonon coupling matrix element (eV/(amu^0.5·Å))
- `pᵢ`: Boltzmann occupation of initial vibrational state `i`
- `χᵢ, χⱼ`: vibrational wavefunctions (eigenfunctions of 1D Schrödinger equation)
- `δ(εᵢ - εⱼ)`: energy-conserving delta function (Gaussian-broadened)
@@ -88,7 +88,7 @@ CarrierCapture.py/
│ └── cli/ # Command-line interface
│ ├── main.py # Click CLI entry point
│ └── commands/ # Subcommands (fit, solve, capture, scan, viz)
-├── tests/ # 88 tests (pytest)
+├── tests/ # 169 tests (pytest)
├── examples/
│ ├── notebooks/ # Jupyter tutorials
│ └── data/ # Example DFT data
@@ -142,8 +142,8 @@ Manages two-state capture calculation.
```python
pot_i: Potential # Initial state (excited)
pot_f: Potential # Final state (ground)
-W: float # Electron-phonon coupling (eV)
-g: int # Degeneracy factor
+W: float # Electron-phonon coupling (eV/(amu^0.5·Å))
+degeneracy: int # Degeneracy factor
overlap_matrix: ndarray # ⟨χᵢ|Q|χⱼ⟩
capture_coefficient: ndarray # C(T) in cm³/s
```
@@ -169,20 +169,21 @@ cc.capture_coefficient # Array of C(T) values
High-throughput parameter scanning.
-**Key Functions**:
+**Key Classes**:
```python
-from carriercapture.analysis.parameter_scan import scan_parameters
+from carriercapture.analysis.parameter_scan import ParameterScanner, ScanParameters
-results = scan_parameters(
+params = ScanParameters(
dQ_range=(0, 25, 25), # (min, max, n_points)
dE_range=(0, 2.5, 10),
- hw_i=0.008,
- hw_f=0.008,
- W=0.068,
+ hbar_omega_i=0.008,
+ hbar_omega_f=0.008,
+ W=0.068, # eV/(amu^0.5·Å), required
volume=1e-21,
temperature=300.0,
- n_jobs=-1 # Parallel execution
)
+scanner = ParameterScanner(params)
+results = scanner.run_harmonic_scan(n_jobs=-1) # Parallel execution
# Access results
results.capture_coefficients # 2D array [dQ, dE]
@@ -333,13 +334,16 @@ def solve(
**Structure**:
```
tests/
-├── test_potential.py # Potential class tests
-├── test_schrodinger.py # Solver validation (analytical solutions)
-├── test_config_coord.py # Capture coefficient workflows
-├── test_parameter_scan.py # High-throughput scanning
-├── test_visualization.py # Plotting functions
-├── test_io.py # File I/O
-└── test_cli.py # Command-line interface
+├── test_potential.py # Potential class tests
+├── test_schrodinger.py # Solver validation (analytical solutions)
+├── test_config_coord.py # Capture coefficient workflows
+├── test_transfer_coord.py # Marcus theory (TransferCoordinate)
+├── test_sommerfeld.py # Sommerfeld factor
+├── test_parameter_scan.py # High-throughput scanning
+├── test_advanced_fitting.py # Potential fitting methods
+├── test_visualization.py # Plotting functions
+├── test_interactive_dashboard.py# Dash dashboard
+└── test_doped_integration.py # doped interface
```
**Run tests**:
@@ -351,7 +355,7 @@ pytest tests/ -x --pdb # Stop on first failure, debug
```
**Current Status** (2026-01-18):
-- 88 tests passing
+- 169 tests (doped-integration tests skip without the optional doped package)
- Core modules: >90% coverage
- Python 3.9-3.12 supported
@@ -436,7 +440,7 @@ OCC_CUTOFF = 1e-5 # Max occupation for partition function convergence
hw = 0.008 eV # 8 meV phonon
dQ = 10.5 # amu^0.5·Å shift
dE = 0.5 eV # Energy difference
-W = 0.068 eV # Electron-phonon coupling
+W = 0.068 eV/(amu^0.5·Å) # Electron-phonon coupling
volume = 1e-21 cm³ # Supercell volume
temperature = 300 K # Room temperature
nev_initial = 180 # Initial state eigenvalues
@@ -508,12 +512,16 @@ E_n = ℏω * (n + 1/2)
**File**: `benchmarks/benchmark_sn_zn.py`
-**Results** (vs CarrierCapture.jl):
-- Initial eigenvalues: 0.005% difference ✓
-- Final eigenvalues: 0.02% difference ✓
-- Capture coefficient (300K): 1.5% difference ✓
+**Results** (three-tier comparison vs CarrierCapture.jl):
+- Tier 1: with Julia's grid/integration conventions emulated, eigenvalues and
+ C(300K) match the Julia reference to ~1e-12 (algorithmic equivalence) ✓
+- Tier 2: native eigenvalues match analytic ℏω(n+½) to ≤1.5e-4 ✓
+- Tier 3: native C(300K) differs from Julia by ~1.5%, entirely due to
+ CarrierCapture.jl's finite-difference spacing (ΔQ/N vs ΔQ/(N−1)) and
+ rectangle-rule overlaps; Python's native numerics are the more accurate ✓
-**Conclusion**: Python matches Julia within ~1-2% (floating-point precision)
+**Conclusion**: identical physics; the small native offset is a Julia grid
+convention, not floating-point noise.
---
@@ -696,7 +704,7 @@ pot.eigenvectors # Array of χₙ(Q), shape: (nev, len(Q))
### Capture Calculation
```python
-cc = ConfigCoordinate(pot_i, pot_f, W=0.068, g=1)
+cc = ConfigCoordinate(pot_i, pot_f, W=0.068, degeneracy=1)
cc.calculate_overlap(Q0=5.0, sigma=0.025, cutoff=0.25)
cc.calculate_capture_coefficient(
volume=1e-21,
diff --git a/README.md b/README.md
index f9787d4..bbbf23d 100644
--- a/README.md
+++ b/README.md
@@ -44,7 +44,7 @@ CarrierCapture.py started life as an automated rewrite of [CarrierCapture.jl](ht
### 🔬 Scientific Validation
- Validated against CarrierCapture.jl
-- Comprehensive test suite (88 tests)
+- Comprehensive test suite (169 tests)
- Tutorial notebooks with real examples
---
@@ -65,14 +65,14 @@ pip install -e ".[dev]"
### Optional Dependencies
```bash
-# Interactive dashboard
-pip install carriercapture[viz]
-
# doped integration (for defect calculations)
pip install carriercapture[doped]
-# All extras (recommended for development)
-pip install -e ".[all]"
+# Jupyter notebook support
+pip install carriercapture[notebook]
+
+# Development tools (tests, linting)
+pip install -e ".[dev]"
```
---
@@ -118,7 +118,7 @@ carriercapture solve excited.json -n 180 -O excited_solved.json
carriercapture capture config.yaml -V 1e-21 --temp-range 100 500 50
# High-throughput parameter scan
-carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 \
+carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 -W 0.05 \
--dE-min 0 --dE-max 2.5 --dE-points 10 \
-j -1 -o scan_results.npz
@@ -159,7 +159,9 @@ pot_final.solve(nev=60)
### Tutorial Notebooks
- **[01_harmonic_sn_zn.ipynb](examples/notebooks/01_harmonic_sn_zn.ipynb)** - Basic workflow with harmonic oscillators
+- **[02_anharmonic_dx_center.ipynb](examples/notebooks/02_anharmonic_dx_center.ipynb)** - Anharmonic potentials (DX center)
- **[03_parameter_scan.ipynb](examples/notebooks/03_parameter_scan.ipynb)** - High-throughput screening
+- **[04_interactive_viz.ipynb](examples/notebooks/04_interactive_viz.ipynb)** - Interactive visualization
Full examples in [`examples/`](examples/) directory with detailed [README](examples/README.md).
@@ -195,7 +197,7 @@ $$C(T) = \frac{V \cdot 2\pi}{\hbar} \cdot g \cdot W^2 \cdot \sum_{i,j} p_i |\lan
Where:
- `V`: supercell volume
- `g`: degeneracy factor
-- `W`: electron-phonon coupling matrix element
+- `W`: electron-phonon coupling matrix element (eV/(amu^0.5·Å))
- `pᵢ`: thermal occupation of initial state `i`
- `χᵢ, χⱼ`: vibrational wavefunctions
- `δ`: energy-conserving delta function (Gaussian broadened)
@@ -222,9 +224,9 @@ pytest tests/ --cov=src/carriercapture --cov-report=html
```
**Test Statistics:**
-- 88 tests passing (53 Phase 3 tests skipped)
+- 169 tests (doped-integration tests skip without the optional doped package)
- Core modules: >90% coverage
-- All tests pass on Python 3.11-3.12
+- Supported Python versions: 3.9-3.12
- CI/CD with GitHub Actions
---
@@ -255,16 +257,26 @@ CarrierCapture.py has been validated against the original [CarrierCapture.jl](ht
- Phonon energy: ℏω = 8 meV
- Configuration coordinate shift: ΔQ = 10.5 amu^0.5·Å
- Energy offset: ΔE = 0.5 eV
-- Electron-phonon coupling: W = 0.068 eV
-
-**Results**:
-| Observable | Python Value | Julia Value | Relative Diff | Status |
-|------------|--------------|-------------|---------------|--------|
-| Initial eigenvalues (E₀) | 0.504000 eV | 0.504001 eV | 0.005% | ✓ PASS |
-| Final eigenvalues (E₀) | 0.004000 eV | 0.004001 eV | 0.02% | ✓ PASS |
-| Capture coefficient (300K) | 1.339×10⁻¹⁰ cm³/s | 1.359×10⁻¹⁰ cm³/s | 1.5% | ✓ PASS |
-
-**Conclusion**: Python implementation matches Julia results within ~1-2% across all observables. Small differences (~0.01-1.5%) are due to floating-point arithmetic differences between language implementations and are well within acceptable tolerances for physical calculations.
+- Electron-phonon coupling: W = 0.068 eV/(amu^0.5·Å)
+
+**Results** (three-tier comparison):
+
+| Tier | Comparison | Max Relative Diff | Status |
+|------|------------|-------------------|--------|
+| 1 | Julia-convention emulation vs Julia reference (eigenvalues) | 3×10⁻¹² | ✓ PASS |
+| 1 | Julia-convention emulation vs Julia reference (C at 300 K) | 5×10⁻¹² | ✓ PASS |
+| 2 | Native eigenvalues vs analytic ℏω(n+½) | 1.5×10⁻⁴ | ✓ PASS |
+| 3 | Native C(300 K) = 1.339×10⁻¹⁰ vs Julia 1.359×10⁻¹⁰ cm³/s | 1.5×10⁻² | ✓ (informational) |
+
+**Conclusion**: When Python is run with CarrierCapture.jl's numerical
+conventions, the two codes agree to near machine precision (Tier 1) — they
+implement identical physics. The ~1.5% native difference (Tier 3) is entirely
+due to two CarrierCapture.jl conventions: its finite-difference kinetic term
+uses grid spacing ΔQ/N while its grid actually has spacing ΔQ/(N−1), and it
+integrates overlaps with the rectangle rule. Python uses the true grid spacing
+and the trapezoid rule, and its native eigenvalues are closer to the analytic
+harmonic result (Tier 2). Both codes converge to the same answer with
+increasing grid density.
### Running the Benchmark
@@ -340,7 +352,7 @@ This project is licensed under the MIT License - see the [LICENSE](LICENSE) file
| Parameter Scanning | ✅ Complete |
| doped Integration | ✅ Complete |
| Documentation | ✅ Complete |
-| Test Coverage | ✅ 88 tests |
+| Test Coverage | ✅ 169 tests |
| PyPI Release | 🔄 Planned |
---
diff --git a/benchmarks/README.md b/benchmarks/README.md
index c10f7d2..a8c105b 100644
--- a/benchmarks/README.md
+++ b/benchmarks/README.md
@@ -36,16 +36,26 @@ This will:
- Phonon energy: ℏω = 8 meV
- Configuration coordinate shift: ΔQ = 10.5 amu^0.5·Å
- Energy offset: ΔE = 0.5 eV
-- Electron-phonon coupling: W = 0.068 eV
-
-**Compared Quantities**:
-- Initial state eigenvalues (first 20)
-- Final state eigenvalues (first 20)
-- Capture coefficient at 300K
-
-**Tolerances**:
-- Eigenvalues: relative tolerance 1e-4
-- Capture coefficient: relative tolerance 1e-2
+- Electron-phonon coupling: W = 0.068 eV/(amu^0.5·Å)
+
+**Three-tier comparison**:
+
+1. **Algorithmic equivalence** (binding, rtol 1e-9 / 1e-6): Python is re-run
+ with CarrierCapture.jl's numerical conventions — its finite-difference
+ kinetic term uses grid spacing ΔQ/N while its `range(Qi, Qf, length=N)`
+ grid actually has spacing ΔQ/(N−1), and it integrates overlaps with the
+ rectangle rule. With those conventions emulated, eigenvalues and C(300K)
+ match the Julia reference to ~1e-12: both codes implement identical physics.
+2. **Native accuracy** (binding, rtol 5e-4): CarrierCapture.py's native
+ eigenvalues (true grid spacing, trapezoid integration) vs the analytic
+ harmonic result E_n = E0 + ℏω(n+½). Measured max 1.5e-4 at npoints=5000.
+3. **Native vs Julia** (informational, rtol 2e-2): the native C(300K) differs
+ from Julia by ~1.5%, entirely attributable to the conventions in Tier 1.
+ Python's native numerics are the more accurate of the two; both converge
+ to the same answer with increasing grid density.
+
+Note: the Tier-1 spacing convention may be worth reporting upstream to
+CarrierCapture.jl.
## Manual Usage
@@ -71,73 +81,29 @@ python benchmarks/benchmark_sn_zn.py
}
```
-**Benchmark Report** (`results/benchmark_report.json`):
-```json
-{
- "test_case": "Sn_Zn in ZnO (Harmonic)",
- "parameters": {...},
- "comparisons": {
- "eigenvalues_initial": {
- "passed": true,
- "max_relative_difference": 1.5e-6,
- "tolerance": 1e-4
- },
- ...
- },
- "overall_passed": true
-}
-```
+**Benchmark Report** (`results/benchmark_report.json`): tiered structure with
+`tier1_algorithmic_equivalence`, `tier2_native_vs_analytic`,
+`tier3_native_vs_julia` (including an `explanation` field), and
+`overall_passed`.
## Expected Results
-If everything works correctly, you should see:
-
```
-============================================================
-CarrierCapture.jl vs CarrierCapture.py Benchmark
-============================================================
-
-Test Case: Sn_Zn in ZnO (Harmonic Approximation)
-
-...
-
-============================================================
-Comparison Results
-============================================================
-
-1. Initial Eigenvalues (first 20 states):
- Max relative diff: 1.23e-06
- Tolerance: 1.00e-04
- Status: ✓ PASS
-
-2. Final Eigenvalues (first 20 states):
- Max relative diff: 2.34e-06
- Tolerance: 1.00e-04
- Status: ✓ PASS
-
-3. Capture Coefficient (300K):
- Python: 1.2345e-12 cm³/s
- Julia: 1.2346e-12 cm³/s
- Relative diff: 8.10e-05
- Tolerance: 1.00e-02
- Status: ✓ PASS
-
-============================================================
-Overall: ✓ ALL TESTS PASSED
-============================================================
+Tier 1: Algorithmic equivalence (Julia conventions emulated)
+ Initial eigenvalues (emulated): max rel diff 2.66e-14 (tol 1e-09) PASS
+ Final eigenvalues (emulated): max rel diff 2.73e-12 (tol 1e-09) PASS
+ C(300K) (emulated): max rel diff 4.92e-12 (tol 1e-06) PASS
+Tier 2: Native accuracy vs analytic E_n = E0 + hw*(n + 1/2)
+ Initial eigenvalues vs analytic: max rel diff 3.55e-05 (tol 5e-04) PASS
+ Final eigenvalues vs analytic: max rel diff 1.49e-04 (tol 5e-04) PASS
+Tier 3: Native C(300K) vs Julia (informational)
+ Python: 1.338989e-10 cm^3/s | Julia: 1.359429e-10 cm^3/s | rel diff 1.50e-02 (tol 2e-02) PASS
+
+Overall: ALL TESTS PASSED
```
## Troubleshooting
-### Julia script fails
-
-The Julia script may need adjustments for the actual CarrierCapture.jl API. Check:
-- Function names (e.g., `potential()` vs `Potential()`)
-- Parameter names and order
-- Module structure (`using CarrierCapture` vs submodules)
-
-Consult [CarrierCapture.jl documentation](https://github.com/WMD-group/CarrierCapture.jl).
-
### Python benchmark fails to find reference data
Make sure you run the Julia reference first:
@@ -145,32 +111,9 @@ Make sure you run the Julia reference first:
julia benchmarks/run_julia_reference.jl
```
-### Tests fail (exceed tolerance)
-
-Small differences are expected due to:
-- Floating-point rounding differences between languages
-- Compiler optimizations
-- BLAS/LAPACK library versions
-
-If differences are > 0.1%, investigate:
-1. Check eigenvalue magnitudes are reasonable (~0.004-0.5 eV)
-2. Verify same parameters used in both implementations
-3. Check grid size and numerical integration settings
-
-## Adding to README
-
-After successful benchmark, add results to main README.md:
+### Tier 1 fails
-```markdown
-## 🔬 Validation Against Julia
-
-Validated against CarrierCapture.jl for the Sn_Zn in ZnO example:
-
-| Observable | Max Relative Diff | Tolerance | Status |
-|------------|-------------------|-----------|--------|
-| Initial eigenvalues | < 1e-5 | 1e-4 | ✓ PASS |
-| Final eigenvalues | < 1e-5 | 1e-4 | ✓ PASS |
-| Capture coefficient (300K) | < 1e-3 | 1e-2 | ✓ PASS |
-
-See `benchmarks/` for benchmark code.
-```
+Tier 1 should agree to ~1e-12; a failure there means a genuine algorithmic
+divergence (not floating-point noise). Verify the same parameters are used in
+both implementations and that the reference JSON was generated by
+`run_julia_reference.jl` unmodified.
diff --git a/benchmarks/benchmark_sn_zn.py b/benchmarks/benchmark_sn_zn.py
index 897dd61..11484b9 100644
--- a/benchmarks/benchmark_sn_zn.py
+++ b/benchmarks/benchmark_sn_zn.py
@@ -4,10 +4,21 @@
===============================================
Runs the Sn_Zn in ZnO example using CarrierCapture.py and compares
-results against Julia reference data to validate numerical accuracy.
-
-Test Case: Sn substituting Zn in ZnO
-Parameters: From examples/notebooks/01_harmonic_sn_zn.ipynb
+results against CarrierCapture.jl reference data in three tiers:
+
+1. Algorithmic equivalence (binding): Python rebuilt with CarrierCapture.jl's
+ numerical conventions must reproduce the Julia reference to near machine
+ precision. This proves both codes implement the same physics.
+2. Native accuracy (binding): Python's native eigenvalues must match the
+ analytic harmonic result E_n = E0 + hw*(n + 1/2).
+3. Native vs Julia (informational): the native results differ from Julia by
+ ~1.5% at this grid, entirely due to two CarrierCapture.jl conventions:
+ its finite-difference kinetic term uses grid spacing dq = (Q_max-Q_min)/N
+ while its grid range(Qi, Qf, length=N) actually has spacing
+ (Q_max-Q_min)/(N-1), and it integrates overlaps with the rectangle rule.
+ Python uses the true grid spacing and the trapezoid rule, and is closer
+ to the analytic eigenvalues. Both codes converge to the same answer as
+ npoints -> infinity.
"""
import json
@@ -15,50 +26,109 @@
from pathlib import Path
import numpy as np
+import scipy.sparse as sp
+from scipy.sparse.linalg import eigsh
# Add src to path for local development
sys.path.insert(0, str(Path(__file__).parent.parent / "src"))
from carriercapture.core.potential import Potential
-from carriercapture.core.config_coord import ConfigCoordinate
+from carriercapture._constants import AMU, HBAR_C, HBAR, K_B
-def compare_arrays(python_vals, julia_vals, name, rtol):
+def solve_julia_convention(pot, nev, npoints, q_range):
"""
- Compare arrays and return detailed comparison.
-
- Parameters
- ----------
- python_vals : array-like
- Python results
- julia_vals : array-like
- Julia reference results
- name : str
- Name of the comparison
- rtol : float
- Relative tolerance
-
- Returns
- -------
- dict
- Comparison results including pass/fail status
+ Solve with CarrierCapture.jl's solve1D_Quantum grid convention.
+
+ CarrierCapture.jl builds its finite-difference kinetic term with
+ dq = (Q_max - Q_min) / N, while its grid range(Qi, Qf, length=N) has
+ true spacing (Q_max - Q_min) / (N - 1). Reproducing that off-by-one
+ here lets us verify algorithmic equivalence to machine precision.
"""
- python_vals = np.array(python_vals)
- julia_vals = np.array(julia_vals)
+ Q = np.linspace(*q_range, npoints)
+ dq = (Q[-1] - Q[0]) / npoints
+ kinetic = (HBAR_C * 1e10) ** 2 / AMU / (2 * dq**2)
+ H = sp.diags(
+ [np.full(npoints - 1, -kinetic), 2 * kinetic + pot(Q), np.full(npoints - 1, -kinetic)],
+ [-1, 0, 1],
+ )
+ vals, vecs = eigsh(H.tocsc(), k=nev, which="SA", tol=0)
+ order = np.argsort(vals)
+ vals, vecs = vals[order], vecs[:, order]
+ vecs /= np.sqrt(dq * np.sum(vecs**2, axis=0))
+ return vals, vecs, Q, dq
+
+
+def capture_julia_convention(params):
+ """C(300K) with Julia's grid convention and rectangle-rule overlaps."""
+ pot_i = Potential.from_harmonic(
+ hw=params["hw"], Q0=0.0, E0=params["dE"],
+ Q_range=tuple(params["Q_range"]), npoints=params["npoints"],
+ )
+ pot_f = Potential.from_harmonic(
+ hw=params["hw"], Q0=params["dQ"], E0=0.0,
+ Q_range=tuple(params["Q_range"]), npoints=params["npoints"],
+ )
+ vals_i, vecs_i, Q, dq = solve_julia_convention(
+ pot_i, params["nev_initial"], params["npoints"], params["Q_range"]
+ )
+ vals_f, vecs_f, _, _ = solve_julia_convention(
+ pot_f, params["nev_final"], params["npoints"], params["Q_range"]
+ )
- abs_diff = np.abs(python_vals - julia_vals)
- rel_diff = abs_diff / np.abs(julia_vals)
+ # Rectangle-rule overlaps S_ij = dq * sum(psi_i * (Q - Q0) * psi_j)
+ # (errstate: subnormal wavefunction tails trip spurious BLAS warnings)
+ operator = Q - params["Q0_crossing"]
+ with np.errstate(all="ignore"):
+ S = dq * (vecs_i * operator[:, None]).T @ vecs_f
+ dE = vals_i[:, None] - vals_f[None, :]
+ delta = np.exp(-(dE**2) / (2 * params["sigma"] ** 2)) / (params["sigma"] * np.sqrt(2 * np.pi))
+ mask = np.abs(dE) < params["cutoff"]
+ S, delta = np.where(mask, S, 0.0), np.where(mask, delta, 0.0)
+
+ beta = 1.0 / (K_B * params["temperature"])
+ occupation = np.exp(-beta * vals_i)
+ occupation /= occupation.sum()
+ prefactor = params["volume"] * 2 * np.pi / HBAR * params["W"] ** 2
+ C = prefactor * np.sum(occupation[:, None] * S**2 * delta)
+ return vals_i, vals_f, C
+
+
+def capture_native(params):
+ """C(300K) with CarrierCapture.py's native (more accurate) numerics."""
+ pot_i = Potential.from_harmonic(
+ hw=params["hw"], Q0=0.0, E0=params["dE"],
+ Q_range=tuple(params["Q_range"]), npoints=params["npoints"],
+ )
+ pot_f = Potential.from_harmonic(
+ hw=params["hw"], Q0=params["dQ"], E0=0.0,
+ Q_range=tuple(params["Q_range"]), npoints=params["npoints"],
+ )
+ pot_i.solve(nev=params["nev_initial"])
+ pot_f.solve(nev=params["nev_final"])
- max_rel_diff = np.max(rel_diff)
- passed = bool(max_rel_diff < rtol)
+ from carriercapture.core.config_coord import ConfigCoordinate
+ cc = ConfigCoordinate(pot_i=pot_i, pot_f=pot_f, W=params["W"])
+ cc.calculate_overlap(Q0=params["Q0_crossing"], cutoff=params["cutoff"], sigma=params["sigma"])
+ cc.calculate_capture_coefficient(
+ volume=params["volume"], temperature=np.array([params["temperature"]])
+ )
+ return pot_i.eigenvalues, pot_f.eigenvalues, cc.capture_coefficient[0]
+
+
+def compare(python_vals, reference_vals, name, rtol, n=None):
+ """Relative comparison of arrays or scalars."""
+ python_vals = np.atleast_1d(np.asarray(python_vals, dtype=float))
+ reference_vals = np.atleast_1d(np.asarray(reference_vals, dtype=float))
+ if n is not None:
+ python_vals, reference_vals = python_vals[:n], reference_vals[:n]
+ rel = np.abs(python_vals - reference_vals) / np.abs(reference_vals)
return {
"name": name,
- "passed": passed,
- "max_relative_difference": float(max_rel_diff),
+ "passed": bool(np.max(rel) < rtol),
+ "max_relative_difference": float(np.max(rel)),
"tolerance": rtol,
- "max_absolute_difference": float(np.max(abs_diff)),
- "mean_relative_difference": float(np.mean(rel_diff))
}
@@ -68,197 +138,84 @@ def main():
print("=" * 60)
print("\nTest Case: Sn_Zn in ZnO (Harmonic Approximation)")
- # Load Julia reference data
ref_path = Path(__file__).parent / "reference_data" / "sn_zn_julia_reference.json"
-
if not ref_path.exists():
- print(f"\n✗ ERROR: Julia reference data not found at {ref_path}")
- print("\nPlease run Julia reference first:")
- print(" julia benchmarks/run_julia_reference.jl")
+ print(f"\nERROR: Julia reference data not found at {ref_path}")
+ print("Run: julia benchmarks/run_julia_reference.jl")
sys.exit(1)
- print(f"\nLoading Julia reference data from:")
- print(f" {ref_path}")
-
with open(ref_path) as f:
- julia_results = json.load(f)
-
- # Extract parameters
- params = julia_results["parameters"]
- print("\nParameters:")
- print(f" ℏω = {params['hw']} eV")
- print(f" ΔQ = {params['dQ']} amu^0.5·Å")
- print(f" ΔE = {params['dE']} eV")
- print(f" W = {params['W']} eV")
- print(f" Volume = {params['volume']} cm³")
- print(f" Temperature = {params['temperature']} K")
- print(f" Grid points = {params['npoints']}")
-
- # Run Python calculations
- print("\nStep 1: Creating harmonic potentials...")
-
- # Initial state (excited): Q0=0.0, E0=0.5 eV
- pot_initial = Potential.from_harmonic(
- hw=params['hw'],
- Q0=0.0,
- E0=params['dE'],
- Q_range=(params['Q_range'][0], params['Q_range'][1]),
- npoints=params['npoints']
- )
-
- # Final state (ground): Q0=10.5, E0=0.0 eV
- pot_final = Potential.from_harmonic(
- hw=params['hw'],
- Q0=params['dQ'],
- E0=0.0,
- Q_range=(params['Q_range'][0], params['Q_range'][1]),
- npoints=params['npoints']
- )
-
- print(" Initial state: E0=0.5 eV, Q0=0.0")
- print(" Final state: E0=0.0 eV, Q0=10.5")
-
- # Solve Schrödinger equation
- print("\nStep 2: Solving Schrödinger equation...")
- pot_initial.solve(nev=params['nev_initial'])
- pot_final.solve(nev=params['nev_final'])
-
- print(f" Initial state: Found {len(pot_initial.eigenvalues)} eigenvalues")
- print(f" E₀ = {pot_initial.eigenvalues[0]:.6f} eV")
- print(f" E₁ = {pot_initial.eigenvalues[1]:.6f} eV")
- print(f" E₂ = {pot_initial.eigenvalues[2]:.6f} eV")
-
- print(f" Final state: Found {len(pot_final.eigenvalues)} eigenvalues")
- print(f" E₀ = {pot_final.eigenvalues[0]:.6f} eV")
- print(f" E₁ = {pot_final.eigenvalues[1]:.6f} eV")
- print(f" E₂ = {pot_final.eigenvalues[2]:.6f} eV")
-
- # Calculate capture coefficient
- print("\nStep 3: Calculating capture coefficient...")
-
- # Use same crossing point as Julia if available
- Q0_crossing = params.get('Q0_crossing', 5.0)
-
- cc = ConfigCoordinate(
- pot_i=pot_initial,
- pot_f=pot_final,
- W=params['W']
+ julia = json.load(f)
+ params = julia["parameters"]
+ n_states = 20
+
+ print(f"\nParameters: hw={params['hw']} eV, dQ={params['dQ']} amu^0.5*Ang, "
+ f"dE={params['dE']} eV, W={params['W']} eV/(amu^0.5*Ang), "
+ f"npoints={params['npoints']}")
+
+ # --- Tier 1: algorithmic equivalence (Julia conventions emulated) ---
+ print("\nTier 1: Algorithmic equivalence (Julia conventions emulated)")
+ em_i, em_f, C_emulated = capture_julia_convention(params)
+ tier1 = [
+ compare(em_i, julia["eigenvalues_initial"], "Initial eigenvalues (emulated)", 1e-9, n_states),
+ compare(em_f, julia["eigenvalues_final"], "Final eigenvalues (emulated)", 1e-9, n_states),
+ compare(C_emulated, julia["capture_coefficient_300K"], "C(300K) (emulated)", 1e-6),
+ ]
+ for c in tier1:
+ print(f" {c['name']}: max rel diff {c['max_relative_difference']:.2e} "
+ f"(tol {c['tolerance']:.0e}) {'PASS' if c['passed'] else 'FAIL'}")
+
+ # --- Tier 2: native accuracy vs analytic harmonic eigenvalues ---
+ print("\nTier 2: Native accuracy vs analytic E_n = E0 + hw*(n + 1/2)")
+ nat_i, nat_f, C_native = capture_native(params)
+ n_arr = np.arange(n_states)
+ analytic_i = params["dE"] + params["hw"] * (n_arr + 0.5)
+ analytic_f = params["hw"] * (n_arr + 0.5)
+ tier2 = [
+ compare(nat_i, analytic_i, "Initial eigenvalues vs analytic", 5e-4, n_states),
+ compare(nat_f, analytic_f, "Final eigenvalues vs analytic", 5e-4, n_states),
+ ]
+ for c in tier2:
+ print(f" {c['name']}: max rel diff {c['max_relative_difference']:.2e} "
+ f"(tol {c['tolerance']:.0e}) {'PASS' if c['passed'] else 'FAIL'}")
+
+ # --- Tier 3: native vs Julia (informational) ---
+ print("\nTier 3: Native C(300K) vs Julia (informational)")
+ C_julia = julia["capture_coefficient_300K"]
+ tier3 = compare(C_native, C_julia, "C(300K) native vs Julia", 2e-2)
+ tier3["python_value"] = float(C_native)
+ tier3["julia_value"] = C_julia
+ tier3["explanation"] = (
+ "The gap is entirely due to CarrierCapture.jl's numerical conventions "
+ "(finite-difference spacing dq = dQ/N vs the true grid spacing dQ/(N-1), "
+ "and rectangle-rule overlap integration), verified by Tier 1. Python's "
+ "native numerics are closer to the analytic eigenvalues (Tier 2). Both "
+ "codes converge to the same answer with increasing npoints."
)
+ print(f" Python: {C_native:.6e} cm^3/s | Julia: {C_julia:.6e} cm^3/s | "
+ f"rel diff {tier3['max_relative_difference']:.2e} "
+ f"(tol {tier3['tolerance']:.0e}) {'PASS' if tier3['passed'] else 'FAIL'}")
- cc.calculate_overlap(Q0=Q0_crossing, sigma=0.025)
- cc.calculate_capture_coefficient(
- volume=params['volume'],
- temperature=np.array([params['temperature']])
- )
-
- C_300K_python = cc.capture_coefficient[0]
-
- print(f" Overlap matrix: {cc.overlap_matrix.shape}")
- print(f" Q0 (crossing) = {Q0_crossing} amu^0.5·Å")
- print(f" C(300K) = {C_300K_python:.6e} cm³/s")
-
- # Compare results
- print("\n" + "=" * 60)
- print("Comparison Results")
- print("=" * 60)
+ overall = all(c["passed"] for c in tier1 + tier2 + [tier3])
- # Compare initial eigenvalues
- n_compare = min(20, len(pot_initial.eigenvalues), len(julia_results["eigenvalues_initial"]))
- eig_initial_comp = compare_arrays(
- pot_initial.eigenvalues[:n_compare],
- julia_results["eigenvalues_initial"][:n_compare],
- "Initial eigenvalues",
- rtol=1e-4
- )
-
- print(f"\n1. Initial Eigenvalues (first {n_compare} states):")
- print(f" Max relative diff: {eig_initial_comp['max_relative_difference']:.2e}")
- print(f" Mean relative diff: {eig_initial_comp['mean_relative_difference']:.2e}")
- print(f" Max absolute diff: {eig_initial_comp['max_absolute_difference']:.2e} eV")
- print(f" Tolerance: {eig_initial_comp['tolerance']:.2e}")
- print(f" Status: {'✓ PASS' if eig_initial_comp['passed'] else '✗ FAIL'}")
-
- # Compare final eigenvalues
- n_compare_f = min(20, len(pot_final.eigenvalues), len(julia_results["eigenvalues_final"]))
- eig_final_comp = compare_arrays(
- pot_final.eigenvalues[:n_compare_f],
- julia_results["eigenvalues_final"][:n_compare_f],
- "Final eigenvalues",
- rtol=1e-4
- )
-
- print(f"\n2. Final Eigenvalues (first {n_compare_f} states):")
- print(f" Max relative diff: {eig_final_comp['max_relative_difference']:.2e}")
- print(f" Mean relative diff: {eig_final_comp['mean_relative_difference']:.2e}")
- print(f" Max absolute diff: {eig_final_comp['max_absolute_difference']:.2e} eV")
- print(f" Tolerance: {eig_final_comp['tolerance']:.2e}")
- print(f" Status: {'✓ PASS' if eig_final_comp['passed'] else '✗ FAIL'}")
-
- # Compare capture coefficient
- C_300K_julia = julia_results["capture_coefficient_300K"]
- capture_rel_diff = abs(C_300K_python - C_300K_julia) / abs(C_300K_julia)
- capture_passed = bool(capture_rel_diff < 1e-2)
-
- capture_comp = {
- "name": "Capture coefficient at 300K",
- "python_value": float(C_300K_python),
- "julia_value": C_300K_julia,
- "relative_difference": float(capture_rel_diff),
- "absolute_difference": float(abs(C_300K_python - C_300K_julia)),
- "tolerance": 1e-2,
- "passed": capture_passed
- }
-
- print(f"\n3. Capture Coefficient (300K):")
- print(f" Python: {capture_comp['python_value']:.6e} cm³/s")
- print(f" Julia: {capture_comp['julia_value']:.6e} cm³/s")
- print(f" Relative diff: {capture_comp['relative_difference']:.2e}")
- print(f" Absolute diff: {capture_comp['absolute_difference']:.2e} cm³/s")
- print(f" Tolerance: {capture_comp['tolerance']:.2e}")
- print(f" Status: {'✓ PASS' if capture_comp['passed'] else '✗ FAIL'}")
-
- # Overall status
- overall_passed = all([
- eig_initial_comp["passed"],
- eig_final_comp["passed"],
- capture_comp["passed"]
- ])
-
- # Save detailed report
report = {
"test_case": "Sn_Zn in ZnO (Harmonic)",
"parameters": params,
- "comparisons": {
- "eigenvalues_initial": eig_initial_comp,
- "eigenvalues_final": eig_final_comp,
- "capture_coefficient": capture_comp
- },
- "overall_passed": overall_passed
+ "tier1_algorithmic_equivalence": tier1,
+ "tier2_native_vs_analytic": tier2,
+ "tier3_native_vs_julia": tier3,
+ "overall_passed": overall,
}
-
report_path = Path(__file__).parent / "results" / "benchmark_report.json"
report_path.parent.mkdir(parents=True, exist_ok=True)
-
with open(report_path, "w") as f:
json.dump(report, f, indent=2)
+ print(f"\nReport saved to {report_path}")
- print(f"\nDetailed report saved to:")
- print(f" {report_path}")
-
- # Final summary
print("\n" + "=" * 60)
- if overall_passed:
- print("Overall: ✓ ALL TESTS PASSED")
- print("\nConclusion: Python implementation matches Julia results")
- print("within numerical precision!")
- else:
- print("Overall: ✗ SOME TESTS FAILED")
- print("\nSome comparisons exceeded tolerance thresholds.")
- print("Check the detailed report for more information.")
+ print("Overall: " + ("ALL TESTS PASSED" if overall else "SOME TESTS FAILED"))
print("=" * 60 + "\n")
-
- # Exit with appropriate code
- sys.exit(0 if overall_passed else 1)
+ sys.exit(0 if overall else 1)
if __name__ == "__main__":
diff --git a/benchmarks/run_julia_reference.jl b/benchmarks/run_julia_reference.jl
index 3c0afb6..b70595f 100644
--- a/benchmarks/run_julia_reference.jl
+++ b/benchmarks/run_julia_reference.jl
@@ -36,7 +36,7 @@ const Q_max = 20.0 # amu^0.5·Å
const npoints = 5000
const nev_initial = 180
const nev_final = 60
-const W = 0.068 # eV - electron-phonon coupling
+const W = 0.068 # eV/(amu^0.5 Angstrom) - electron-phonon coupling
const volume = 1e-21 # cm³
const temperature = 300.0 # K
const cut_off = 0.25 # eV
diff --git a/docs/TEST_RESULTS.md b/docs/TEST_RESULTS.md
deleted file mode 100644
index ffe791f..0000000
--- a/docs/TEST_RESULTS.md
+++ /dev/null
@@ -1,198 +0,0 @@
-# Documentation Build Test Results
-
-**Date**: 2026-01-17
-**MkDocs Version**: 1.6.1
-**Material Theme Version**: 9.7.1
-**Build Status**: ✅ **SUCCESS**
-
----
-
-## Summary
-
-The documentation successfully builds and renders all completed sections:
-
-- ✅ Landing page (index.md)
-- ✅ Getting Started section (4 pages)
-- ✅ User Guide section (7 pages)
-- ✅ API Reference (5 pages)
-
-**Total pages**: 17 pages
-**Build time**: ~2.3 seconds
-**Site size**: ~2.8 MB
-
----
-
-## Build Statistics
-
-### Page Sizes
-
-| Section | Pages | Size Range | Status |
-|---------|-------|------------|--------|
-| Getting Started | 4 | 64-108 KB | ✅ All render correctly |
-| User Guide | 7 | 140-204 KB | ✅ All render correctly |
-| API Reference | 5 | 72-180 KB | ✅ All render correctly |
-
-### Largest Pages
-1. `user-guide/parameter-scanning` - 204 KB
-2. `user-guide/visualization` - 188 KB
-3. `api/core` - 180 KB
-4. `user-guide/capture-coefficients` - 164 KB
-5. `user-guide/doped-integration` - 164 KB
-
----
-
-## Issues Found
-
-### 1. Expected Warnings (Missing Files)
-
-The following warnings are **expected** and will be resolved in future phases:
-
-**Missing Theory Section (Phase 4):**
-- `theory/multiphonon-theory.md`
-- `theory/configuration-coordinates.md`
-- `theory/equations.md`
-- `theory/references.md`
-
-**Missing Tutorial Section (Phase 6):**
-- `tutorials/index.md`
-- `tutorials/01-harmonic-oscillator.md`
-- `tutorials/02-dx-center.md`
-- `tutorials/03-parameter-scan.md`
-- `tutorials/04-interactive-dashboard.md`
-
-**Missing Examples Section (Phase 7):**
-- `examples/gallery.md`
-- `examples/notebooks.md`
-
-**Missing Development Section (Phase 8):**
-- `development/contributing.md`
-- `development/testing.md`
-- `development/architecture.md`
-
-**Missing Changelog (Phase 9):**
-- `changelog.md`
-
-### 2. Anchor Link Warnings
-
-Some internal links to CLI reference sections use URL-encoded anchors that don't exactly match:
-
-**Affected links:**
-- `#capture---calculate-capture-coefficient` (referenced from 3 pages)
-- `#scan---parameter-scan` (referenced from 1 page)
-- `#fit---fit-potential-energy-surface` (referenced from 1 page)
-- `#solve---solve-schr%C3%B6dinger-equation` (referenced from 1 page)
-- `#viz---interactive-dashboard` (referenced from 1 page)
-- `#plot---static-plots` (referenced from 1 page)
-
-**Impact**: Links may not jump to exact section, but page still loads correctly.
-
-**Fix**: Update anchor format in `docs/api/cli.md` headers or update links to match generated anchors.
-
----
-
-## Verification Checklist
-
-### ✅ Build Process
-- [x] `mkdocs build` completes without errors
-- [x] All markdown files converted to HTML
-- [x] Site directory created successfully
-- [x] No critical errors or failures
-
-### ✅ Navigation
-- [x] Sidebar navigation renders
-- [x] Top navigation tabs work
-- [x] Search bar present
-- [x] Dark/light mode toggle present
-- [x] Breadcrumbs functional
-
-### ✅ Content Rendering
-- [x] All pages have substantial content
-- [x] Code blocks formatted correctly
-- [x] Tables render properly
-- [x] Lists and nested lists work
-- [x] Blockquotes render correctly
-
-### ✅ Features
-- [x] Table of contents (TOC) generated
-- [x] Code syntax highlighting
-- [x] Anchor links in headers
-- [x] Search indexing
-- [x] Responsive design (mobile-friendly)
-
-### ⏳ Pending Verification
-- [ ] LaTeX equations (will verify when Theory section added)
-- [ ] Admonitions (info, warning, etc.) - used in some pages
-- [ ] Mermaid diagrams (if any)
-- [ ] External links functionality
-- [ ] Image embedding (no images yet)
-
----
-
-## Server Information
-
-**Local development server**: http://127.0.0.1:8000/CarrierCapture.py/
-**GitHub Pages URL**: https://wmd-group.github.io/CarrierCapture.py/ (not yet deployed)
-
----
-
-## Recommendations
-
-### Immediate Actions
-
-1. **Fix anchor links in api/cli.md**
- - Option A: Change headers to use simpler IDs
- - Option B: Update links to match generated anchors
-
-2. **Add placeholder pages** (optional)
- - Create stub files for missing sections to avoid navigation errors
- - Example: Empty `docs/theory/multiphonon-theory.md` with "Coming soon"
-
-### Before GitHub Pages Deployment
-
-1. **Complete remaining sections** (Phases 4, 6, 7, 8, 9)
-2. **Test all internal links**
-3. **Verify search functionality** with full content
-4. **Check mobile responsiveness**
-5. **Review meta tags and SEO**
-6. **Test GitHub Actions workflow**
-
-### Documentation Quality
-
-1. **Add images/diagrams** where helpful
- - Configuration coordinate diagrams
- - Workflow diagrams
- - Example plots
-
-2. **Add more code examples** with expected output
-
-3. **Cross-reference consistency check**
- - Verify all internal links work
- - Ensure consistent terminology
-
----
-
-## Testing Commands
-
-```bash
-# Build documentation
-mkdocs build
-
-# Serve locally
-mkdocs serve
-
-# Build with strict mode (fail on warnings)
-mkdocs build --strict
-
-# Deploy to GitHub Pages (requires permissions)
-mkdocs gh-deploy
-```
-
----
-
-## Conclusion
-
-✅ **The documentation build is successful and ready for continued development.**
-
-All completed sections (Phases 1, 2, 3, 5) render correctly with proper formatting, navigation, and search functionality. The remaining phases can be implemented without any blocking issues.
-
-**Next steps**: Continue with Phase 4 (Theory), Phase 7 (Examples), Phase 8 (Development), or Phase 9 (Changelog) as planned.
diff --git a/docs/api/cli.md b/docs/api/cli.md
index b84291f..e65aaac 100644
--- a/docs/api/cli.md
+++ b/docs/api/cli.md
@@ -181,7 +181,7 @@ carriercapture capture [CONFIG_FILE] [OPTIONS]
|--------|------|---------|-------------|
| `--pot-i` | path | - | Initial state potential file |
| `--pot-f` | path | - | Final state potential file |
-| `-W, --coupling` | float | - | Electron-phonon coupling (eV) |
+| `-W, --coupling` | float | - | Electron-phonon coupling (eV/(amu^0.5·Å)) |
| `-g, --degeneracy` | int | 1 | Degeneracy factor |
| `-V, --volume` | float | - | Supercell volume (cm³) |
| `--temp-range` | float×3 | 100 500 50 | Temperature range: `T_min T_max n_points` (K) |
diff --git a/docs/api/index.md b/docs/api/index.md
index f62e1c3..37a07f0 100644
--- a/docs/api/index.md
+++ b/docs/api/index.md
@@ -35,7 +35,7 @@ from carriercapture.analysis import ParameterScanner, ScanParameters
from carriercapture.visualization import plot_potential, plot_capture_coefficient
# I/O
-from carriercapture.io import load_potential, save_results
+from carriercapture.io import load_potential_from_file, write_capture_results
# doped integration
from carriercapture.io.doped_interface import load_defect_entry, create_potential_from_doped
diff --git a/docs/api/io.md b/docs/api/io.md
index ddfe5bc..88885c6 100644
--- a/docs/api/io.md
+++ b/docs/api/io.md
@@ -40,11 +40,11 @@ Integration with the [doped](https://github.com/SMTG-Bham/doped) package for def
### Loading Data
```python
-from carriercapture.io import load_potential, read_csv_data
+from carriercapture.io import load_potential_from_file, read_csv
from carriercapture.core import Potential
# Load from CSV file
-Q_data, E_data = read_csv_data('potential_data.csv')
+Q_data, E_data = read_csv('potential_data.csv')
# Create potential from data
pot = Potential(Q_data=Q_data, E_data=E_data)
@@ -54,7 +54,7 @@ pot.fit(fit_type='spline', order=4, smoothness=0.001)
### Saving Results
```python
-from carriercapture.io import save_results
+from carriercapture.io import write_capture_results
from carriercapture.core import ConfigCoordinate
# After calculating capture coefficient
@@ -62,10 +62,10 @@ cc = ConfigCoordinate(...)
cc.calculate_capture_coefficient(...)
# Save to JSON
-save_results(cc, 'results.json', format='json')
+write_capture_results(cc, 'results.json', file_format='json')
-# Save to HDF5
-save_results(cc, 'results.h5', format='hdf5')
+# Save to CSV
+write_capture_results(cc, 'results.csv', file_format='csv')
```
### doped Integration
diff --git a/docs/examples/gallery.md b/docs/examples/gallery.md
index 707f911..8489c86 100644
--- a/docs/examples/gallery.md
+++ b/docs/examples/gallery.md
@@ -226,6 +226,7 @@ from carriercapture.visualization import plot_scan_heatmap
params = ScanParameters(
dQ_range=(0, 25, 25),
dE_range=(0, 2.5, 10),
+ W=0.05,
hbar_omega_i=0.008,
hbar_omega_f=0.008,
temperature=300.0,
@@ -562,7 +563,7 @@ fig.add_annotation(
fig.add_annotation(
xref='paper', yref='paper',
x=0.95, y=0.95,
- text="T = 300K
W = 0.205 eV",
+ text="T = 300K
W = 0.205 eV/(amu^0.5·Å)",
showarrow=False,
bgcolor='white',
bordercolor='black',
diff --git a/docs/getting-started/basic-concepts.md b/docs/getting-started/basic-concepts.md
index dd111a2..8580e1b 100644
--- a/docs/getting-started/basic-concepts.md
+++ b/docs/getting-started/basic-concepts.md
@@ -230,7 +230,7 @@ print(f"C(300K) = {cc.capture_coefficient[0]:.3e} cm³/s")
| Phonon energy | $\hbar\omega$ | eV | Vibrational quantum |
| Displacement | $\Delta Q$ | amu$^{0.5}$·Å | Shift between states |
| Energy difference | $\Delta E$ | eV | Vertical separation |
-| Coupling | $W$ | eV | Electron-phonon interaction |
+| Coupling | $W$ | eV/(amu^0.5·Å) | Electron-phonon interaction |
| Capture coefficient | $C$ | cm³/s | Capture rate constant |
| Huang-Rhys factor | $S$ | - | Effective phonon number |
| Reorganization energy | $\lambda$ | eV | Energy to relax |
diff --git a/docs/getting-started/first-calculation.md b/docs/getting-started/first-calculation.md
index 0b6aa1d..35bc048 100644
--- a/docs/getting-started/first-calculation.md
+++ b/docs/getting-started/first-calculation.md
@@ -144,7 +144,7 @@ Calculate overlaps between initial and final state wavefunctions.
cc = ConfigCoordinate(
pot_i=pot_initial,
pot_f=pot_final,
- W=0.068, # Electron-phonon coupling (eV)
+ W=0.068, # Electron-phonon coupling (eV/(amu^0.5·Å))
name="Sn_Zn"
)
diff --git a/docs/getting-started/installation.md b/docs/getting-started/installation.md
index 4c30e47..da1ff50 100644
--- a/docs/getting-started/installation.md
+++ b/docs/getting-started/installation.md
@@ -72,7 +72,7 @@ pip install -e ".[dev]"
For the Dash dashboard and advanced plotting:
```bash
-pip install carriercapture[viz]
+pip install carriercapture[doped]
```
This includes:
diff --git a/docs/getting-started/quick-start.md b/docs/getting-started/quick-start.md
index 6f671c9..9154f22 100644
--- a/docs/getting-started/quick-start.md
+++ b/docs/getting-started/quick-start.md
@@ -46,7 +46,7 @@ pot_final.solve(nev=60) # 60 eigenvalues
cc = ConfigCoordinate(
pot_i=pot_initial,
pot_f=pot_final,
- W=0.068 # Electron-phonon coupling (eV)
+ W=0.068 # Electron-phonon coupling (eV/(amu^0.5·Å))
)
# Calculate overlap matrix
diff --git a/docs/index.md b/docs/index.md
index 7476b19..ce88c25 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -41,7 +41,7 @@
### Installation
```bash
-pip install carriercapture
+pip install carriercapture # PyPI release coming soon; install from source for now
```
For development or from source:
@@ -91,7 +91,7 @@ carriercapture solve excited.json -n 180 -O excited_solved.json
carriercapture capture config.yaml -V 1e-21 --temp-range 100 500 50
# High-throughput parameter scan
-carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 \
+carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 -W 0.05 \
--dE-min 0 --dE-max 2.5 --dE-points 10 \
-j -1 -o scan_results.npz
@@ -147,7 +147,7 @@ Where:
- $V$: supercell volume
- $g$: degeneracy factor
-- $W$: electron-phonon coupling matrix element
+- $W$: electron-phonon coupling matrix element (eV/(amu^0.5·Å))
- $p_i$: thermal occupation of initial state $i$
- $\chi_i, \chi_j$: vibrational wavefunctions
- $\delta$: energy-conserving delta function (Gaussian broadened)
@@ -169,7 +169,7 @@ Where:
| Parameter Scanning | ✅ Complete |
| doped Integration | ✅ Complete |
| Documentation | ✅ Complete |
-| Test Coverage | ✅ 88 tests |
+| Test Coverage | ✅ 169 tests |
| PyPI Release | 🔄 Planned |
---
diff --git a/docs/user-guide/capture-coefficients.md b/docs/user-guide/capture-coefficients.md
index dd2bd0d..583856f 100644
--- a/docs/user-guide/capture-coefficients.md
+++ b/docs/user-guide/capture-coefficients.md
@@ -26,7 +26,7 @@ import numpy as np
cc = ConfigCoordinate(
pot_i=pot_initial, # Initial state (e.g., neutral defect)
pot_f=pot_final, # Final state (e.g., charged defect)
- W=0.205, # Electron-phonon coupling (eV)
+ W=0.205, # Electron-phonon coupling (eV/(amu^0.5·Å))
degeneracy=1 # Degeneracy factor
)
@@ -127,7 +127,7 @@ from carriercapture.core import ConfigCoordinate
cc = ConfigCoordinate(
pot_i=pot_i,
pot_f=pot_f,
- W=0.205, # Coupling strength (eV)
+ W=0.205, # Coupling strength (eV/(amu^0.5·Å))
degeneracy=1 # g = 1 for non-degenerate states
)
```
@@ -138,7 +138,7 @@ cc = ConfigCoordinate(
|-----------|------|-------------|
| `pot_i` | Potential | Initial state (before capture) |
| `pot_f` | Potential | Final state (after capture) |
-| `W` | float | Electron-phonon coupling (eV) |
+| `W` | float | Electron-phonon coupling (eV/(amu^0.5·Å)) |
| `degeneracy` | int | Degeneracy factor $g$ |
**Determining W:**
@@ -484,6 +484,45 @@ for g, C in results.items():
---
+## Sommerfeld Factor for Charged Defects
+
+The capture coefficient above assumes a neutral defect. For a **charged**
+defect, the Coulomb interaction between the free carrier and the defect
+enhances (attractive) or suppresses (repulsive) the carrier density at the
+defect site. Correct the coefficient with the Sommerfeld factor $s(T)$:
+
+$$C_{\text{charged}}(T) = s(T) \cdot C(T)$$
+
+```python
+from carriercapture import sommerfeld_parameter
+
+s = sommerfeld_parameter(
+ temperature=cc.temperature,
+ Z=-1, # defect charge / carrier charge: Z < 0 attractive, Z > 0 repulsive
+ m_eff=0.2, # carrier effective mass (units of m_e)
+ eps0=10.0, # relative static dielectric constant
+)
+
+C_charged = s * cc.capture_coefficient
+```
+
+Two methods are available:
+
+- `method="Integrate"` (default): Maxwell-Boltzmann thermal average of the
+ exact Coulomb enhancement factor
+- `method="Analytic"`: low-temperature limits of Pässler,
+ phys. stat. sol. (b) 78, 625 (1976) — attractive
+ $s = 4\sqrt{\theta/\pi}$, repulsive
+ $s = (8/\sqrt{3})\,\theta^{2/3} e^{-3\theta^{1/3}}$,
+ with $\theta = \pi^2 Z^2 E_R / k_B T$ and scaled Rydberg
+ $E_R = m^* \mathrm{Ry} / \varepsilon_0^2$
+
+**Caveats**: assumes a parabolic, non-degenerate band and static screening.
+See Alkauskas et al., Phys. Rev. B 90, 075202 (2014), Sec. II.F for
+discussion.
+
+---
+
## Best Practices
### 1. Convergence Checks
@@ -571,9 +610,8 @@ print(f"RMSE (log scale): {rmse:.3f}")
**Debug:**
```python
-print(f"Volume: {cc.volume:.3e} cm³ (typical: 1e-21)")
-print(f"W: {cc.W:.3f} eV (typical: 0.1-0.5)")
-print(f"Q0: {cc.Q0:.2f} amu^0.5·Å")
+# volume and Q0 are passed to the calculation methods, not stored on cc
+print(f"W: {cc.W:.4f} eV/(amu^0.5·Å)")
print(f"States: {len(pot_i.eigenvalues)} initial, {len(pot_f.eigenvalues)} final")
```
diff --git a/docs/user-guide/cli-usage.md b/docs/user-guide/cli-usage.md
index 467d939..e660bda 100644
--- a/docs/user-guide/cli-usage.md
+++ b/docs/user-guide/cli-usage.md
@@ -83,7 +83,7 @@ Initial potential: 180 states
Final potential: 60 states
Capture parameters:
- W (coupling): 0.205 eV
+ W (coupling): 0.205 eV/(amu^0.5·Å)
g (degeneracy): 1
V (volume): 1.00e-21 cm³
Q0: 10.0 amu^0.5·Å
@@ -121,7 +121,7 @@ potential_final:
file: ground_solved.json
capture:
- W: 0.205 # eV
+ W: 0.205 # eV/(amu^0.5·Å)
degeneracy: 1
volume: 1.0e-21 # cm³
Q0: 10.0 # amu^0.5·Å
@@ -145,7 +145,7 @@ Screen materials across (ΔQ, ΔE) space:
```bash
# Run parameter scan
-carriercapture scan \
+carriercapture scan -W 0.05 \
--dQ-min 0 --dQ-max 25 --dQ-points 25 \
--dE-min 0 --dE-max 2.5 --dE-points 10 \
--hbar-omega-i 0.008 --hbar-omega-f 0.008 \
@@ -326,7 +326,7 @@ carriercapture fit data.dat -o fit.json -v >> workflow.log 2>&1
```bash
# Use all CPU cores
-carriercapture scan \
+carriercapture scan -W 0.05 \
--dQ-min 0 --dQ-max 25 --dQ-points 50 \
--dE-min 0 --dE-max 2.5 --dE-points 20 \
-j -1 \
diff --git a/docs/user-guide/configuration.md b/docs/user-guide/configuration.md
index 8dce06d..a8d3383 100644
--- a/docs/user-guide/configuration.md
+++ b/docs/user-guide/configuration.md
@@ -25,7 +25,7 @@ potential_final:
file: ground_solved.json
capture:
- W: 0.205 # Electron-phonon coupling (eV)
+ W: 0.205 # Electron-phonon coupling (eV/(amu^0.5·Å))
degeneracy: 1 # Degeneracy factor
volume: 1.0e-21 # Supercell volume (cm³)
Q0: 10.0 # Coordinate shift (amu^0.5·Å)
@@ -78,7 +78,7 @@ potential_final:
# Capture calculation parameters
capture:
# Coupling parameters
- W: 0.205 # Electron-phonon coupling (eV)
+ W: 0.205 # Electron-phonon coupling (eV/(amu^0.5·Å))
degeneracy: 1 # Degeneracy factor (g)
# System parameters
@@ -175,7 +175,7 @@ Core capture calculation parameters:
```yaml
capture:
# ---- Required parameters ----
- W: 0.205 # Electron-phonon coupling (eV)
+ W: 0.205 # Electron-phonon coupling (eV/(amu^0.5·Å))
# Typical range: 0.1 - 0.5 eV
volume: 1.0e-21 # Supercell volume (cm³)
@@ -232,7 +232,7 @@ output:
```yaml
capture:
- W: 0.205 # eV
+ W: 0.205 # eV/(amu^0.5·Å)
# Guidelines:
# - Must be calculated for the specific defect transition;
diff --git a/docs/user-guide/parameter-scanning.md b/docs/user-guide/parameter-scanning.md
index 70ca391..9b2aaa7 100644
--- a/docs/user-guide/parameter-scanning.md
+++ b/docs/user-guide/parameter-scanning.md
@@ -28,6 +28,7 @@ import numpy as np
params = ScanParameters(
dQ_range=(0, 25, 25), # ΔQ: 0-25 amu^0.5·Å, 25 points
dE_range=(0, 2.5, 10), # ΔE: 0-2.5 eV, 10 points
+ W=0.05,
hbar_omega_i=0.008, # Initial state phonon (eV)
hbar_omega_f=0.008, # Final state phonon (eV)
temperature=300.0, # Temperature (K)
@@ -66,6 +67,7 @@ from carriercapture.analysis import ScanParameters
params = ScanParameters(
dQ_range=(Q_min, Q_max, n_points), # ΔQ scan range
dE_range=(E_min, E_max, n_points), # ΔE scan range
+ W=0.05,
hbar_omega_i=0.008, # Initial phonon (eV)
hbar_omega_f=0.008, # Final phonon (eV)
temperature=300.0, # Temperature (K)
@@ -151,9 +153,9 @@ hbar_omega_f = 0.008 # eV
hbar_omega_i = 0.010 # Initial state
hbar_omega_f = 0.008 # Final state
-# From DFT phonon calculations:
-from carriercapture.analysis import estimate_phonon_energy
-# hbar_omega = estimate_phonon_energy(phonopy_yaml_file)
+# From a Q-E dataset (e.g. VASP path calculations):
+from carriercapture.io.doped_interface import estimate_phonon_frequency
+# hbar_omega = estimate_phonon_frequency(Q_data, E_data)['hw']
```
### Grid Resolution
@@ -165,6 +167,7 @@ Trade-off between accuracy and computational cost:
params_coarse = ScanParameters(
dQ_range=(0, 25, 10), # 10 ΔQ points
dE_range=(0, 2.5, 5), # 5 ΔE points → 50 total
+ W=0.05,
# ... other params
)
@@ -172,12 +175,14 @@ params_coarse = ScanParameters(
params_medium = ScanParameters(
dQ_range=(0, 25, 25), # 25 ΔQ points
dE_range=(0, 2.5, 10), # 10 ΔE points → 250 total
+ W=0.05,
)
# Fine scan (high-res, ~2000 points)
params_fine = ScanParameters(
dQ_range=(0, 25, 50), # 50 ΔQ points
dE_range=(0, 2.5, 20), # 20 ΔE points → 1000 total
+ W=0.05,
)
# Estimate computation time
@@ -203,6 +208,7 @@ from carriercapture.analysis import ParameterScanner, ScanParameters
params = ScanParameters(
dQ_range=(0, 25, 25),
dE_range=(0, 2.5, 10),
+ W=0.05,
hbar_omega_i=0.008,
hbar_omega_f=0.008,
temperature=300.0,
@@ -286,7 +292,7 @@ print(f"ΔE grid: {dE_grid.shape}") # (10,)
print(f"C matrix: {C_matrix.shape}") # (25, 10)
# Access parameters
-print(f"Temperature: {results.temperature} K")
+print(f"Temperature: {results.parameters.temperature} K")
print(f"Phonon ℏω_i: {results.hbar_omega_i} eV")
print(f"Volume: {results.volume:.2e} cm³")
```
@@ -396,7 +402,7 @@ fig.show()
# Customized
fig = plot_scan_heatmap(
results,
- title=f"Capture Coefficient at {results.temperature}K",
+ title=f"Capture Coefficient at {params.temperature}K",
log_scale=True,
colorscale='Viridis',
width=900,
@@ -496,7 +502,7 @@ results = ScanResult.load('scan_results.h5', format='hdf5')
# Access data immediately
print(f"Loaded scan: {results.capture_coefficients.shape}")
-print(f"Temperature: {results.temperature} K")
+print(f"Temperature: {results.parameters.temperature} K")
```
### Combining Multiple Scans
@@ -510,6 +516,7 @@ for T in temperatures:
params = ScanParameters(
dQ_range=(0, 25, 25),
dE_range=(0, 2.5, 10),
+ W=0.05,
temperature=T,
# ... other params
)
@@ -551,6 +558,7 @@ for T in temperatures:
params = ScanParameters(
dQ_range=(0, 25, 25),
dE_range=(0, 2.5, 10),
+ W=0.05,
temperature=T,
# ... other params
)
@@ -572,6 +580,7 @@ Refine grid near interesting regions:
params_coarse = ScanParameters(
dQ_range=(0, 30, 10), # Coarse: 10 points
dE_range=(0, 3, 6), # Coarse: 6 points
+ W=0.05,
# ... other params
)
scanner_coarse = ParameterScanner(params_coarse)
@@ -592,6 +601,7 @@ dE_range_fine = (max(0, dE_max-0.5), dE_max+0.5, 15)
params_fine = ScanParameters(
dQ_range=dQ_range_fine,
dE_range=dE_range_fine,
+ W=0.05,
# ... other params
)
scanner_fine = ParameterScanner(params_fine)
@@ -677,6 +687,7 @@ carriercapture scan-plot scan_results.npz --log-scale -o heatmap.html
params_coarse = ScanParameters(
dQ_range=(0, 30, 10),
dE_range=(0, 3, 6),
+ W=0.05,
# ... other params
)
results_coarse = ParameterScanner(params_coarse).run_harmonic_scan(n_jobs=-1)
@@ -699,6 +710,7 @@ for nev_i in nev_values:
params_test = ScanParameters(
dQ_range=(15, 15, 1), # Single point
dE_range=(1.0, 1.0, 1),
+ W=0.05,
nev_initial=nev_i,
nev_final=60,
# ... other params
@@ -783,6 +795,7 @@ from carriercapture.visualization import plot_scan_heatmap
params = ScanParameters(
dQ_range=(0, 25, 25), # 25 points: 0-25 amu^0.5·Å
dE_range=(0, 2.5, 10), # 10 points: 0-2.5 eV
+ W=0.05,
hbar_omega_i=0.008, # 8 meV (typical for ZnO)
hbar_omega_f=0.008, # 8 meV
temperature=300.0, # Room temperature
diff --git a/docs/user-guide/potentials.md b/docs/user-guide/potentials.md
index e5e1bcf..83bc1cc 100644
--- a/docs/user-guide/potentials.md
+++ b/docs/user-guide/potentials.md
@@ -43,7 +43,7 @@ pot = Potential(Q_data=Q_data, E_data=E_data, name="My Potential")
Load potential data from files:
```python
-from carriercapture.io import load_potential, read_potential_data
+from carriercapture.io import load_potential_from_file, read_potential_data
# Option 1: Read Q-E data from CSV/DAT
Q, E = read_potential_data('potential.csv')
diff --git a/docs/user-guide/visualization.md b/docs/user-guide/visualization.md
index e113fb4..2baf2f6 100644
--- a/docs/user-guide/visualization.md
+++ b/docs/user-guide/visualization.md
@@ -806,7 +806,7 @@ def create_all_figures(pot_i, pot_f, cc, output_dir='figures'):
fig2.write_html(f'{output_dir}/02_potential_final.html')
# 3. CC diagram
- fig3 = plot_configuration_coordinate(pot_i, pot_f, Q0=cc.Q0, show_crossing=True)
+ fig3 = plot_configuration_coordinate(pot_i, pot_f, Q0=Q0, show_crossing=True)
fig3.write_html(f'{output_dir}/03_cc_diagram.html')
# 4. Overlap matrix
diff --git a/examples/README.md b/examples/README.md
index 08a83fe..76d8842 100644
--- a/examples/README.md
+++ b/examples/README.md
@@ -113,7 +113,7 @@ carriercapture capture config.yaml -V 1e-21 --temp-range 100 500 50 -O results.j
carriercapture plot results.json --show
# 5. Parameter scan
-carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 \
+carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 -W 0.05 \
--dE-min 0 --dE-max 2.5 --dE-points 10 \
-j 4 -o scan_results.npz
diff --git a/examples/notebooks/03_parameter_scan.ipynb b/examples/notebooks/03_parameter_scan.ipynb
index 1d37e86..abc25bb 100644
--- a/examples/notebooks/03_parameter_scan.ipynb
+++ b/examples/notebooks/03_parameter_scan.ipynb
@@ -58,6 +58,7 @@
"params = ScanParameters(\n",
" dQ_range=(0, 25, 25), # (min, max, n_points)\n",
" dE_range=(0, 2.5, 10), # (min, max, n_points)\n",
+ " W=0.05, # e-ph coupling (eV/(amu^0.5*A)) - calculate for your defect\n",
" hbar_omega_i=0.008, # Initial state phonon (eV)\n",
" hbar_omega_f=0.008, # Final state phonon (eV)\n",
" temperature=300.0, # Temperature (K)\n",
@@ -315,6 +316,7 @@
" params_hw = ScanParameters(\n",
" dQ_range=(0, 25, 15), # Coarser grid for speed\n",
" dE_range=(0, 2.5, 8),\n",
+ " W=0.05, # e-ph coupling (eV/(amu^0.5*A)) - calculate for your defect\n",
" hbar_omega_i=hw,\n",
" hbar_omega_f=hw,\n",
" temperature=300.0,\n",
@@ -356,4 +358,4 @@
},
"nbformat": 4,
"nbformat_minor": 4
-}
+}
\ No newline at end of file
diff --git a/mkdocs.yml b/mkdocs.yml
index 6f3865a..1fa9c80 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -78,17 +78,6 @@ nav:
- doped Integration: user-guide/doped-integration.md
- CLI Usage: user-guide/cli-usage.md
- Configuration: user-guide/configuration.md
- - Theory:
- - Multiphonon Theory: theory/multiphonon-theory.md
- - Configuration Coordinates: theory/configuration-coordinates.md
- - Equations: theory/equations.md
- - References: theory/references.md
- - Tutorials:
- - Overview: tutorials/index.md
- - 1. Harmonic Oscillator: tutorials/01-harmonic-oscillator.md
- - 2. DX Center: tutorials/02-dx-center.md
- - 3. Parameter Scan: tutorials/03-parameter-scan.md
- - 4. Interactive Dashboard: tutorials/04-interactive-dashboard.md
- API Reference:
- Overview: api/index.md
- Core: api/core.md
@@ -99,8 +88,3 @@ nav:
- Examples:
- Gallery: examples/gallery.md
- Notebooks: examples/notebooks.md
- - Development:
- - Contributing: development/contributing.md
- - Testing: development/testing.md
- - Architecture: development/architecture.md
- - Changelog: changelog.md
diff --git a/src/carriercapture/__init__.py b/src/carriercapture/__init__.py
index ef5502d..85a851d 100644
--- a/src/carriercapture/__init__.py
+++ b/src/carriercapture/__init__.py
@@ -7,10 +7,12 @@
from .__version__ import __version__
from .core import Potential, ConfigCoordinate, TransferCoordinate
+from .analysis.sommerfeld import sommerfeld_parameter
__all__ = [
"__version__",
"Potential",
"ConfigCoordinate",
"TransferCoordinate",
+ "sommerfeld_parameter",
]
diff --git a/src/carriercapture/_constants.py b/src/carriercapture/_constants.py
index 8275455..2698d0d 100644
--- a/src/carriercapture/_constants.py
+++ b/src/carriercapture/_constants.py
@@ -13,6 +13,9 @@
HBAR_C = 0.19732697e-6 # eV·m - reduced Planck constant times speed of light
HBAR = 6.582119514e-16 # eV·s - reduced Planck constant
K_B = 8.6173303e-5 # eV/K - Boltzmann constant
+M_E_C2 = 0.5109989461e6 # eV - electron rest energy m_e·c²
+FINE_STRUCTURE = 7.2973525693e-3 # dimensionless fine-structure constant
+RYDBERG_EV = 0.5 * M_E_C2 * FINE_STRUCTURE**2 # eV - Rydberg energy (13.6057 eV)
# Conversion factors
EV_TO_HARTREE = 1.0 / 27.21138602 # Hartree to eV conversion
diff --git a/src/carriercapture/analysis/__init__.py b/src/carriercapture/analysis/__init__.py
index ab3c2f2..d3d1f04 100644
--- a/src/carriercapture/analysis/__init__.py
+++ b/src/carriercapture/analysis/__init__.py
@@ -5,9 +5,11 @@
ScanResult,
ParameterScanner,
)
+from .sommerfeld import sommerfeld_parameter
__all__ = [
"ScanParameters",
"ScanResult",
"ParameterScanner",
+ "sommerfeld_parameter",
]
diff --git a/src/carriercapture/analysis/parameter_scan.py b/src/carriercapture/analysis/parameter_scan.py
index 9447d17..cb5fa9c 100644
--- a/src/carriercapture/analysis/parameter_scan.py
+++ b/src/carriercapture/analysis/parameter_scan.py
@@ -6,12 +6,13 @@
and progress reporting.
"""
-from typing import Dict, List, Optional, Tuple, Callable, Any, Union
+from typing import Dict, Optional, Tuple, Any, Union
from pathlib import Path
+import json
import numpy as np
from numpy.typing import NDArray
import warnings
-from dataclasses import dataclass, field
+from dataclasses import dataclass, field, asdict
from carriercapture.core.potential import Potential
from carriercapture.core.config_coord import ConfigCoordinate
@@ -32,6 +33,10 @@ class ScanParameters:
ℏω for initial state (eV). If tuple: (min, max, n_points)
hbar_omega_f : float or tuple
ℏω for final state (eV). If tuple: (min, max, n_points)
+ W : float
+ Electron-phonon coupling matrix element (eV/(amu^0.5·Å)).
+ Must be calculated for the defect transition of interest;
+ since C ∝ W², it scales the whole map uniformly.
temperature : float or NDArray
Temperature(s) for calculation (K)
volume : float
@@ -47,6 +52,7 @@ class ScanParameters:
dE_range: Tuple[float, float, int]
hbar_omega_i: Union[float, Tuple[float, float, int]] = 0.008 # 8 meV default
hbar_omega_f: Union[float, Tuple[float, float, int]] = 0.008
+ W: Optional[float] = None
temperature: Union[float, NDArray[np.float64]] = 300.0
volume: float = 1e-21
degeneracy: int = 1
@@ -98,6 +104,11 @@ def save(self, filepath: Union[str, Path], format: str = "npz") -> None:
"""
filepath = Path(filepath)
+ params_dict = asdict(self.parameters)
+ if isinstance(params_dict.get("temperature"), np.ndarray):
+ params_dict["temperature"] = params_dict["temperature"].tolist()
+ params_json = json.dumps(params_dict)
+
if format == "npz":
np.savez_compressed(
filepath,
@@ -105,7 +116,7 @@ def save(self, filepath: Union[str, Path], format: str = "npz") -> None:
dE_grid=self.dE_grid,
capture_coefficients=self.capture_coefficients,
barrier_heights=self.barrier_heights,
- # Store parameters as dict
+ parameters_json=params_json,
**{f"param_{k}": v for k, v in self.metadata.items()}
)
elif format == "hdf5":
@@ -116,7 +127,7 @@ def save(self, filepath: Union[str, Path], format: str = "npz") -> None:
f.create_dataset('dE_grid', data=self.dE_grid)
f.create_dataset('capture_coefficients', data=self.capture_coefficients)
f.create_dataset('barrier_heights', data=self.barrier_heights)
- # Store metadata as attributes
+ f.attrs['parameters_json'] = params_json
for k, v in self.metadata.items():
f.attrs[k] = v
except ImportError:
@@ -147,19 +158,14 @@ def load(cls, filepath: Union[str, Path], format: str = "npz") -> "ScanResult":
data = np.load(filepath, allow_pickle=True)
metadata = {k.replace('param_', ''): v for k, v in data.items()
if k.startswith('param_')}
-
- # Reconstruct ScanParameters (simplified)
- params = ScanParameters(
- dQ_range=(0, 0, 0), # Placeholder
- dE_range=(0, 0, 0), # Placeholder
- )
+ params_json = str(data['parameters_json']) if 'parameters_json' in data else None
return cls(
dQ_grid=data['dQ_grid'],
dE_grid=data['dE_grid'],
capture_coefficients=data['capture_coefficients'],
barrier_heights=data['barrier_heights'],
- parameters=params,
+ parameters=cls._params_from_json(params_json),
metadata=dict(metadata)
)
elif format == "hdf5":
@@ -167,16 +173,13 @@ def load(cls, filepath: Union[str, Path], format: str = "npz") -> "ScanResult":
import h5py
with h5py.File(filepath, 'r') as f:
metadata = dict(f.attrs)
- params = ScanParameters(
- dQ_range=(0, 0, 0),
- dE_range=(0, 0, 0),
- )
+ params_json = metadata.pop('parameters_json', None)
return cls(
dQ_grid=f['dQ_grid'][:],
dE_grid=f['dE_grid'][:],
capture_coefficients=f['capture_coefficients'][:],
barrier_heights=f['barrier_heights'][:],
- parameters=params,
+ parameters=cls._params_from_json(params_json),
metadata=metadata
)
except ImportError:
@@ -184,6 +187,20 @@ def load(cls, filepath: Union[str, Path], format: str = "npz") -> "ScanResult":
else:
raise ValueError(f"Unknown format: {format}. Use 'npz' or 'hdf5'")
+ @staticmethod
+ def _params_from_json(params_json: Optional[str]) -> ScanParameters:
+ """Reconstruct ScanParameters from the JSON stored by save()."""
+ if not params_json:
+ # File predates parameter persistence
+ return ScanParameters(dQ_range=(0, 0, 0), dE_range=(0, 0, 0))
+ d = json.loads(params_json)
+ for key in ("dQ_range", "dE_range", "hbar_omega_i", "hbar_omega_f"):
+ if isinstance(d.get(key), list):
+ d[key] = tuple(d[key])
+ if isinstance(d.get("temperature"), list):
+ d["temperature"] = np.array(d["temperature"])
+ return ScanParameters(**d)
+
class ParameterScanner:
"""
@@ -204,6 +221,7 @@ class ParameterScanner:
>>> params = ScanParameters(
... dQ_range=(0, 25, 25),
... dE_range=(0, 2.5, 10),
+ ... W=0.05,
... )
>>> scanner = ParameterScanner(params)
>>> results = scanner.run_harmonic_scan(n_jobs=4)
@@ -211,6 +229,12 @@ class ParameterScanner:
"""
def __init__(self, params: ScanParameters, verbose: bool = True):
+ if params.W is None:
+ raise ValueError(
+ "ScanParameters.W (electron-phonon coupling, eV/(amu^0.5·Å)) "
+ "must be set: it must be calculated for the defect transition "
+ "of interest, and C ∝ W² scales the whole scan."
+ )
self.params = params
self.verbose = verbose
@@ -301,52 +325,6 @@ def _create_harmonic_potentials(
return pot_i, pot_f
- def _calculate_W_coupling(
- self,
- hbar_omega_f: float,
- dQ: float,
- dE: float
- ) -> float:
- """
- Calculate electron-phonon coupling W.
-
- Uses activationless Marcus regime formula:
- Q_m = sqrt(E0 / a) where a = (amu/2)(ℏω/(ℏc))^2
- W = 0.068 / (Q0 - Q_m)
-
- Parameters
- ----------
- hbar_omega_f : float
- ℏω for final state (eV)
- dQ : float
- Horizontal shift (amu^0.5·Å)
- dE : float
- Vertical shift (eV)
-
- Returns
- -------
- float
- Electron-phonon coupling W (eV)
- """
- from carriercapture._constants import AMU, HBAR_C
-
- # Calculate force constant a
- a = (AMU / 2) * (hbar_omega_f / (HBAR_C * 1e10)) ** 2
-
- # Marcus activationless point
- if dE > 0 and a > 0:
- Q_m = np.sqrt(dE / a)
- else:
- Q_m = 0.0
-
- # Calculate W
- if abs(dQ - Q_m) > 1e-6:
- W = 0.068 / abs(dQ - Q_m)
- else:
- W = 0.068 # Default value
-
- return W
-
def _calculate_single_point(
self,
hbar_omega_i: float,
@@ -381,14 +359,11 @@ def _calculate_single_point(
hbar_omega_i, hbar_omega_f, dQ, dE
)
- # Calculate W coupling
- W = self._calculate_W_coupling(hbar_omega_f, dQ, dE)
-
# Create ConfigCoordinate
cc = ConfigCoordinate(
pot_i=pot_i,
pot_f=pot_f,
- W=W,
+ W=self.params.W,
degeneracy=self.params.degeneracy
)
@@ -421,10 +396,9 @@ def _calculate_single_point(
from carriercapture.core.potential import find_crossing
crossing_Q, crossing_E = find_crossing(pot_f, pot_i)
barrier_height = crossing_E - dE
- except Exception as e:
- # If can't find crossing, set high barrier
- warnings.warn(f"Could not find crossing at dQ={dQ:.2f}, dE={dE:.2f}: {e}. Using default barrier height of 50.0 eV.")
- barrier_height = 50.0
+ except (ValueError, RuntimeError) as e:
+ warnings.warn(f"Could not find crossing at dQ={dQ:.2f}, dE={dE:.2f}: {e}")
+ barrier_height = np.nan
return capture_coeff, barrier_height
@@ -508,40 +482,20 @@ def run_harmonic_scan(
# Parallel execution
try:
from joblib import Parallel, delayed
+ except ImportError:
+ raise ImportError("joblib not installed. Install with: pip install joblib")
- if show_progress:
- try:
- from rich.progress import Progress
- with Progress() as progress:
- task = progress.add_task("Scanning...", total=len(params_list))
-
- def _wrapped_calc(args):
- result = self._calculate_single_point(*args[:4])
- progress.update(task, advance=1)
- return result + args[4:]
-
- results = Parallel(n_jobs=n_jobs)(
- delayed(_wrapped_calc)(p) for p in params_list
- )
- except ImportError:
- # No rich, just use joblib's verbose
- results = Parallel(n_jobs=n_jobs, verbose=10 if self.verbose else 0)(
- delayed(self._calculate_single_point)(*p[:4]) + (p[4], p[5])
- for p in params_list
- )
- else:
- results = Parallel(n_jobs=n_jobs)(
- delayed(self._calculate_single_point)(*p[:4]) + (p[4], p[5])
- for p in params_list
- )
+ def _calc_with_index(args):
+ return self._calculate_single_point(*args[:4]) + args[4:]
- # Unpack results
- for capture_coeff, barrier_height, i, j in results:
- capture_coeffs[i, j] = capture_coeff
- barrier_heights[i, j] = barrier_height
+ verbose = 10 if (show_progress and self.verbose) else 0
+ results = Parallel(n_jobs=n_jobs, verbose=verbose)(
+ delayed(_calc_with_index)(p) for p in params_list
+ )
- except ImportError:
- raise ImportError("joblib not installed. Install with: pip install joblib")
+ for capture_coeff, barrier_height, i, j in results:
+ capture_coeffs[i, j] = capture_coeff
+ barrier_heights[i, j] = barrier_height
if self.verbose:
n_success = np.sum(~np.isnan(capture_coeffs))
diff --git a/src/carriercapture/analysis/sommerfeld.py b/src/carriercapture/analysis/sommerfeld.py
new file mode 100644
index 0000000..0947854
--- /dev/null
+++ b/src/carriercapture/analysis/sommerfeld.py
@@ -0,0 +1,107 @@
+"""
+Sommerfeld (Coulomb enhancement) factor for capture by charged defects.
+
+The capture coefficients computed by ConfigCoordinate assume a neutral
+defect. For a charged defect, the Coulomb interaction between the free
+carrier and the defect enhances (attractive) or suppresses (repulsive)
+the carrier density at the defect site. The corrected coefficient is:
+
+ C_charged(T) = s(T) * C(T)
+
+following Pässler, phys. stat. sol. (b) 78, 625 (1976) and
+Alkauskas et al., Phys. Rev. B 90, 075202 (2014), Sec. II.F.
+"""
+
+from typing import Union
+
+import numpy as np
+from numpy.typing import NDArray
+from scipy.integrate import quad
+
+from .._constants import K_B, RYDBERG_EV
+
+
+def sommerfeld_parameter(
+ temperature: Union[float, NDArray[np.float64]],
+ Z: int,
+ m_eff: float,
+ eps0: float,
+ method: str = "Integrate",
+) -> Union[float, NDArray[np.float64]]:
+ """
+ Calculate the Sommerfeld factor s(T) for capture by a charged defect.
+
+ Parameters
+ ----------
+ temperature : float or NDArray[np.float64]
+ Temperature(s) (K)
+ Z : int
+ Defect charge divided by carrier charge: Z < 0 for attractive
+ centers, Z > 0 for repulsive, Z = 0 returns 1
+ m_eff : float
+ Carrier effective mass (units of the electron mass)
+ eps0 : float
+ Relative static dielectric constant
+ method : str, default="Integrate"
+ "Integrate": Maxwell-Boltzmann thermal average of the exact
+ Coulomb enhancement factor. "Analytic": low-temperature limits
+ of Pässler (1976).
+
+ Returns
+ -------
+ s : float or NDArray[np.float64]
+ Sommerfeld factor (same shape as temperature)
+
+ Notes
+ -----
+ With the scaled Rydberg E_R = m_eff·Ry/eps0² and θ = π²Z²E_R/(k_B·T),
+ the analytic limits are:
+
+ attractive: s = 4·sqrt(θ/π)
+ repulsive: s = (8/sqrt(3))·θ^(2/3)·exp(-3·θ^(1/3))
+
+ The integral form averages f(E) = η/(1 - exp(-η)) (attractive) or
+ η/(exp(η) - 1) (repulsive), with η(E) = 2π|Z|·sqrt(E_R/E), over a
+ Maxwell-Boltzmann distribution.
+
+ Examples
+ --------
+ >>> sommerfeld_parameter(300, Z=-1, m_eff=0.2, eps0=10, method="Analytic")
+ 7.27...
+ """
+ T = np.asarray(temperature, dtype=float)
+ scalar_input = T.ndim == 0
+ T = np.atleast_1d(T)
+
+ if Z == 0:
+ s = np.ones_like(T)
+ return s[0] if scalar_input else s
+
+ E_R = m_eff * RYDBERG_EV / eps0**2 # scaled Rydberg (eV)
+
+ if method.lower().startswith("a"):
+ theta = np.pi**2 * Z**2 * E_R / (K_B * T)
+ if Z < 0:
+ s = 4.0 * np.sqrt(theta / np.pi)
+ else:
+ s = (8.0 / np.sqrt(3.0)) * theta ** (2.0 / 3.0) * np.exp(-3.0 * theta ** (1.0 / 3.0))
+ elif method.lower().startswith("i"):
+ # s(T) = over MB distribution, in x = E/(k_B T):
+ # s = (2/sqrt(pi)) * ∫ f(x·kT) √x e^(-x) dx
+ def enhancement(x: float, kT: float) -> float:
+ eta = 2.0 * np.pi * abs(Z) * np.sqrt(E_R / (x * kT))
+ if Z < 0:
+ return eta / -np.expm1(-eta)
+ return eta / np.expm1(eta)
+
+ s = np.array(
+ [
+ (2.0 / np.sqrt(np.pi))
+ * quad(lambda x: enhancement(x, K_B * t) * np.sqrt(x) * np.exp(-x), 0, np.inf)[0]
+ for t in T
+ ]
+ )
+ else:
+ raise ValueError(f"Unknown method '{method}': use 'Integrate' or 'Analytic'")
+
+ return s[0] if scalar_input else s
diff --git a/src/carriercapture/cli/commands/capture.py b/src/carriercapture/cli/commands/capture.py
index 364478f..5ba8cae 100644
--- a/src/carriercapture/cli/commands/capture.py
+++ b/src/carriercapture/cli/commands/capture.py
@@ -32,13 +32,13 @@
@click.option(
"-W", "--coupling",
type=float,
- help="Electron-phonon coupling (eV)"
+ help="Electron-phonon coupling (eV/(amu^0.5·Å))"
)
@click.option(
"-g", "--degeneracy",
type=int,
- default=1,
- help="Degeneracy factor"
+ default=None,
+ help="Degeneracy factor. Default: 1"
)
@click.option(
"-V", "--volume",
@@ -59,14 +59,14 @@
@click.option(
"--cutoff",
type=float,
- default=0.25,
- help="Energy cutoff for overlaps (eV)"
+ default=None,
+ help="Energy cutoff for overlaps (eV). Default: 0.25"
)
@click.option(
"--sigma",
type=float,
- default=0.025,
- help="Gaussian delta width (eV)"
+ default=None,
+ help="Gaussian delta width (eV). Default: 0.025"
)
@click.option(
"-o", "--output",
@@ -176,8 +176,8 @@ def capture_cmd(ctx, config_file, pot_i, pot_f, coupling, degeneracy, volume,
pot_f = pot_f_config.get('file')
if coupling is None:
coupling = capture_config.get('W')
- if degeneracy == 1: # Default value
- degeneracy = capture_config.get('degeneracy', 1)
+ if degeneracy is None:
+ degeneracy = capture_config.get('degeneracy')
if not volume:
volume = capture_config.get('volume')
if not temp_range:
@@ -189,10 +189,10 @@ def capture_cmd(ctx, config_file, pot_i, pot_f, coupling, degeneracy, volume,
temp_range = (t_min, t_max, n_points)
if q0 is None:
q0 = capture_config.get('Q0')
- if cutoff == 0.25: # Default value
- cutoff = capture_config.get('cutoff', 0.25)
- if sigma == 0.025: # Default value
- sigma = capture_config.get('sigma', 0.025)
+ if cutoff is None:
+ cutoff = capture_config.get('cutoff')
+ if sigma is None:
+ sigma = capture_config.get('sigma')
# Handle doped integration mode
if doped:
@@ -348,6 +348,14 @@ def capture_cmd(ctx, config_file, pot_i, pot_f, coupling, degeneracy, volume,
if not temp_range:
temp_range = (100, 500, 50) # Default
+ # Apply defaults for options not set on the command line or in a config
+ if degeneracy is None:
+ degeneracy = 1
+ if cutoff is None:
+ cutoff = 0.25
+ if sigma is None:
+ sigma = 0.025
+
# Load potentials (skip if already loaded from doped)
if not doped:
if verbose > 0:
@@ -382,7 +390,7 @@ def capture_cmd(ctx, config_file, pot_i, pot_f, coupling, degeneracy, volume,
click.echo(f"\nInitial potential: {len(potential_i.eigenvalues)} states")
click.echo(f"Final potential: {len(potential_f.eigenvalues)} states")
click.echo(f"\nCapture parameters:")
- click.echo(f" W (coupling): {coupling} eV")
+ click.echo(f" W (coupling): {coupling} eV/(amu^0.5·Å)")
click.echo(f" g (degeneracy): {degeneracy}")
click.echo(f" V (volume): {volume:.2e} cm³")
click.echo(f" Q0: {q0} amu^0.5·Å")
@@ -441,18 +449,9 @@ def capture_cmd(ctx, config_file, pot_i, pot_f, coupling, degeneracy, volume,
click.echo(f"\nSaving results to: {output}")
try:
- # Detect format from extension
- ext = output.suffix.lower()
- format_map = {
- '.json': 'json',
- '.yaml': 'yaml',
- '.yml': 'yaml',
- '.csv': 'csv',
- '.npz': 'npz',
- }
- file_format = format_map.get(ext, 'json')
-
- write_capture_results(cc, output, file_format=file_format)
+ from carriercapture.io.readers import detect_format
+
+ write_capture_results(cc, output, file_format=detect_format(output))
if verbose > 0:
click.echo("✓ Saved successfully")
except Exception as e:
@@ -466,38 +465,19 @@ def capture_cmd(ctx, config_file, pot_i, pot_f, coupling, degeneracy, volume,
# Plot if requested
if plot:
try:
- import matplotlib.pyplot as plt
-
- fig, ax = plt.subplots(figsize=(10, 6))
-
- # Arrhenius plot: log(C) vs 1000/T
- x = 1000.0 / temperature
- y = np.log10(cc.capture_coefficient)
+ from carriercapture.visualization import plot_capture_coefficient
- ax.plot(x, y, 'o-', linewidth=2, markersize=5)
-
- ax.set_xlabel("1000/T (K$^{-1}$)", fontsize=12)
- ax.set_ylabel("log$_{10}$(C) [cm$^3$/s]", fontsize=12)
- ax.set_title("Capture Coefficient (Arrhenius Plot)", fontsize=14)
- ax.grid(True, alpha=0.3)
-
- # Add temperature labels on top axis
- ax2 = ax.twiny()
- temps_label = [100, 200, 300, 400, 500]
- ax2.set_xlim(ax.get_xlim())
- ax2.set_xticks([1000/t for t in temps_label if temp_range[0] <= t <= temp_range[1]])
- ax2.set_xticklabels([f"{t}K" for t in temps_label if temp_range[0] <= t <= temp_range[1]])
-
- plt.tight_layout()
+ fig = plot_capture_coefficient(cc)
if plot_output:
- plt.savefig(plot_output, dpi=300, bbox_inches='tight')
+ if plot_output.suffix.lower() in ('.html', '.htm'):
+ fig.write_html(str(plot_output))
+ else:
+ fig.write_image(str(plot_output)) # requires kaleido
if verbose > 0:
click.echo(f"✓ Plot saved to: {plot_output}")
else:
- plt.show()
+ fig.show()
- except ImportError:
- click.echo("Warning: matplotlib not available for plotting", err=True)
except Exception as e:
click.echo(f"Error during plotting: {e}", err=True)
diff --git a/src/carriercapture/cli/commands/scan.py b/src/carriercapture/cli/commands/scan.py
index 215f014..6b4df7f 100644
--- a/src/carriercapture/cli/commands/scan.py
+++ b/src/carriercapture/cli/commands/scan.py
@@ -118,10 +118,16 @@
is_flag=True,
help="Disable progress bar"
)
+@click.option(
+ "-W", "--coupling",
+ type=float,
+ required=True,
+ help="Electron-phonon coupling W (eV/(amu^0.5·Å)); scales C by W²"
+)
@click.pass_context
def scan_cmd(ctx, dq_min, dq_max, dq_points, de_min, de_max, de_points,
hbar_omega_i, hbar_omega_f, temperature, volume, degeneracy,
- sigma, cutoff, nev_i, nev_f, n_jobs, output, no_progress):
+ sigma, cutoff, nev_i, nev_f, n_jobs, output, no_progress, coupling):
"""
Run high-throughput parameter scan.
@@ -132,22 +138,22 @@ def scan_cmd(ctx, dq_min, dq_max, dq_points, de_min, de_max, de_points,
\\b
Examples:
# Basic scan over ΔQ and ΔE
- $ carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 \\
- --dQ-min 0 --dE-max 2.5 --dE-points 10 \\
+ $ carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 -W 0.05 \\
+ --dE-min 0 --dE-max 2.5 --dE-points 10 \\
-o scan_results.npz
# Parallel scan with 4 cores
- $ carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 \\
+ $ carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 -W 0.05 \\
--dE-min 0 --dE-max 2.5 --dE-points 10 \\
-j 4 -o results.npz
# Use all available cores
- $ carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 50 \\
+ $ carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 50 -W 0.05 \\
--dE-min 0 --dE-max 2.5 --dE-points 20 \\
-j -1 -o results.npz
# Custom phonon frequencies
- $ carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 \\
+ $ carriercapture scan --dQ-min 0 --dQ-max 25 --dQ-points 25 -W 0.05 \\
--dE-min 0 --dE-max 2.5 --dE-points 10 \\
--hbar-omega-i 0.010 --hbar-omega-f 0.010 \\
-o results.npz
@@ -175,6 +181,7 @@ def scan_cmd(ctx, dq_min, dq_max, dq_points, de_min, de_max, de_points,
dE_range=(de_min, de_max, de_points),
hbar_omega_i=hbar_omega_i,
hbar_omega_f=hbar_omega_f,
+ W=coupling,
temperature=temperature,
volume=volume,
degeneracy=degeneracy,
@@ -294,7 +301,7 @@ def scan_plot_cmd(ctx, scan_file, plot_type, log_scale, output, show):
try:
from carriercapture.analysis.parameter_scan import ScanResult
- import plotly.graph_objects as go
+ from carriercapture.visualization import plot_scan_heatmap
# Load results
if scan_file.suffix.lower() in ['.h5', '.hdf5']:
@@ -305,43 +312,7 @@ def scan_plot_cmd(ctx, scan_file, plot_type, log_scale, output, show):
if verbose > 0:
click.echo(f"Loaded grid: {results.dQ_grid.shape[0]} × {results.dE_grid.shape[0]}")
- # Prepare data
- Z = results.capture_coefficients
- if log_scale:
- Z = np.log10(Z + 1e-30) # Add small epsilon to avoid log(0)
- colorbar_title = "log₁₀(C) [cm³/s]"
- else:
- colorbar_title = "C [cm³/s]"
-
- # Create figure
- fig = go.Figure()
-
- if plot_type in ["heatmap", "both"]:
- fig.add_trace(go.Heatmap(
- x=results.dE_grid,
- y=results.dQ_grid,
- z=Z,
- colorscale='Viridis',
- colorbar=dict(title=colorbar_title),
- ))
-
- if plot_type in ["contour", "both"]:
- fig.add_trace(go.Contour(
- x=results.dE_grid,
- y=results.dQ_grid,
- z=Z,
- colorscale='Viridis',
- colorbar=dict(title=colorbar_title),
- ))
-
- fig.update_layout(
- title="Parameter Scan: Capture Coefficient",
- xaxis_title="ΔE (eV)",
- yaxis_title="ΔQ (amu0.5·Å)",
- template="plotly_white",
- width=800,
- height=700,
- )
+ fig = plot_scan_heatmap(results, plot_type=plot_type, log_scale=log_scale)
# Save or show
if output:
diff --git a/src/carriercapture/core/config_coord.py b/src/carriercapture/core/config_coord.py
index e75d68a..ea017c5 100644
--- a/src/carriercapture/core/config_coord.py
+++ b/src/carriercapture/core/config_coord.py
@@ -33,7 +33,7 @@ class ConfigCoordinate:
pot_f : Potential
Final state potential
W : float
- Electron-phonon coupling matrix element (eV). Must be calculated
+ Electron-phonon coupling matrix element (eV/(amu^0.5·Å)). Must be calculated
for the specific defect transition; there is no meaningful default.
degeneracy : int
Degeneracy factor
@@ -82,7 +82,7 @@ def __init__(
name : str, default=""
Identifier for this configuration coordinate
W : float, default=0.0
- Electron-phonon coupling matrix element (eV). Since the capture
+ Electron-phonon coupling matrix element (eV/(amu^0.5·Å)). Since the capture
coefficient scales as W², the default of 0.0 yields identically
zero capture coefficients — a value calculated for the specific
defect transition must be supplied for physical results.
diff --git a/src/carriercapture/core/schrodinger.py b/src/carriercapture/core/schrodinger.py
index 7c1a114..01bf8bc 100644
--- a/src/carriercapture/core/schrodinger.py
+++ b/src/carriercapture/core/schrodinger.py
@@ -56,10 +56,10 @@ def build_hamiltonian_1d(
Examples
--------
- >>> def harmonic(Q):
- ... return 0.5 * 0.02 * Q**2
+ >>> from carriercapture.core.potential import Potential
+ >>> pot = Potential.from_harmonic(hw=0.02, Q0=0.0, E0=0.0)
>>> Q = np.linspace(-10, 10, 1000)
- >>> H = build_hamiltonian_1d(harmonic, Q)
+ >>> H = build_hamiltonian_1d(pot, Q)
>>> H.shape
(1000, 1000)
"""
@@ -70,38 +70,10 @@ def build_hamiltonian_1d(
if not np.allclose(np.diff(Q), h):
raise ValueError("Q grid must be uniformly spaced for finite difference method")
- # Dimensional Schrödinger equation: [-ℏ²/(2m) d²/dQ² + V(Q)]ψ = Eψ
- # With m=1 amu, Q in amu^0.5·Å, V in eV, E in eV
- #
- # We need ℏ²/(2m) in units of eV·(amu^0.5·Å)²
- # ℏc = 0.19732697e-6 eV·m = 0.19732697e-6 * 1e10 eV·Å = 1973.2697 eV·Å
- # ℏ = ℏc/c, but c cancels in ℏ²/m, so use ℏc directly
- # Actually: ℏ²/(2m) = (ℏc)²/(2mc²) where mc² is in eV
- #
- # For m = 1 amu, mc² = AMU = 931.494e6 eV
- # ℏc in eV·Å = HBAR_C * 1e10
- # So: ℏ²/(2m) = (HBAR_C * 1e10)² / (2 * AMU) [units: eV·Å²/amu]
- #
- # But Q is in amu^0.5·Å, so dQ has units amu^0.5·Å
- # d²/dQ² has units 1/(amu^0.5·Å)² = 1/(amu·Å²)
- # So ℏ²/(2m) * d²/dQ² has units: [eV·Å²/amu] * [1/(amu·Å²)] = eV/amu²
- #
- # Wait, that doesn't work. Let me reconsider...
- #
- # Actually, in configuration coordinate space with Q in amu^0.5·Å:
- # The kinetic energy operator is: T = -ℏ²/(2μ) d²/dQ²
- # where μ = 1 amu is the effective mass
- #
- # ℏ² = (ℏc)²/c² but we're working with ℏc = 0.197e-6 eV·m
- # Let's use: ℏ² = (ℏc * 1e10 Å/m)² = (1973.27 eV·Å)²
- # And: μ = 1 amu = 931.494e6 eV/c²
- #
- # So: ℏ²/(2μ) = (1973.27)² / (2 * 931.494e6) eV·Å²·c²/eV
- # = (1973.27)² / (2 * 931.494e6) Ų·c²
- #
- # Hmm, this is getting messy. Let me just use the factor from Julia:
+ # Mass-weighted coordinates (Q in amu^0.5·Å, unit effective mass):
+ # T = -(ℏ²/2) d²/dQ², with ℏ² = (ℏc)²/(AMU·c²) = (HBAR_C·1e10)²/AMU
+ # expressed in eV·amu·Å², so T comes out in eV.
factor = (1.0 / AMU) * (HBAR_C * 1e10) ** 2
- # This has units: [c²/eV] * [eV·Å]² = Ų·c²
# Kinetic energy coefficient: ℏ²/(2m h²)
kinetic_coeff = factor / (2.0 * h**2)
@@ -145,8 +117,10 @@ def normalize_wavefunctions(
Notes
-----
- Uses trapezoidal rule for numerical integration:
+ Uses the rectangle rule for numerical integration:
∫|ψ|² dQ ≈ h * Σ|ψ|²
+ (equivalent to the trapezoid rule here, since the hard-wall boundary
+ conditions make ψ vanish at the grid edges)
Examples
--------
@@ -156,8 +130,7 @@ def normalize_wavefunctions(
>>> np.allclose(np.sum(wf_norm**2, axis=1) * h, 1.0)
True
"""
- # Integrate |ψ|² using trapezoidal rule (simplified as uniform grid)
- # norm² = ∫|ψ|² dQ = h * Σ|ψ|²
+ # norm² = ∫|ψ|² dQ ≈ h * Σ|ψ|² (rectangle rule, uniform grid)
norms_sq = grid_spacing * np.sum(wavefunctions**2, axis=1)
norms = np.sqrt(norms_sq)
@@ -229,15 +202,14 @@ def solve_schrodinger_1d(
Examples
--------
- Harmonic oscillator:
+ Harmonic oscillator with ℏω = 20 meV (E_n = ℏω(n + 1/2)):
- >>> def harmonic(Q):
- ... hw = 0.02 # ℏω = 20 meV
- ... return 0.5 * hw * Q**2
+ >>> from carriercapture.core.potential import Potential
+ >>> pot = Potential.from_harmonic(hw=0.02, Q0=0.0, E0=0.0)
>>> Q = np.linspace(-20, 20, 5000)
- >>> eigenvalues, eigenvectors = solve_schrodinger_1d(harmonic, Q, nev=10)
- >>> eigenvalues[:3] # First 3 eigenvalues
- array([0.01, 0.03, 0.05]) # E_n = ℏω(n + 1/2)
+ >>> eigenvalues, eigenvectors = solve_schrodinger_1d(pot, Q, nev=10)
+ >>> np.round(eigenvalues[:3], 6)
+ array([0.01, 0.03, 0.05])
See Also
--------
diff --git a/src/carriercapture/core/transfer_coord.py b/src/carriercapture/core/transfer_coord.py
index 5c12c92..fe30ac7 100644
--- a/src/carriercapture/core/transfer_coord.py
+++ b/src/carriercapture/core/transfer_coord.py
@@ -52,7 +52,7 @@ class TransferCoordinate:
>>> pot_1 = Potential.from_harmonic(hw=0.02, Q0=0.0, E0=0.0)
>>> pot_2 = Potential.from_harmonic(hw=0.02, Q0=5.0, E0=0.1)
>>> tc = TransferCoordinate(pot_1, pot_2, name="hole_transfer")
- >>> tc.get_coupling()
+ >>> tc.get_coupling(H_ab=0.01) # H_ab from an electronic-structure calculation
>>> tc.get_reorganization_energy()
>>> tc.get_transfer_rate(temperature=np.linspace(100, 500, 50))
"""
@@ -62,6 +62,7 @@ def __init__(
pot_1: Potential,
pot_2: Potential,
name: str = "",
+ coupling: Optional[float] = None,
):
"""
Initialize a TransferCoordinate.
@@ -74,6 +75,10 @@ def __init__(
Second diabatic state potential
name : str, default=""
Identifier for this transfer coordinate
+ coupling : float, optional
+ Electronic coupling H_ab (eV) from an electronic-structure
+ calculation. Required (here or via get_coupling) before
+ computing transfer rates.
"""
self.name = name
self.pot_1 = pot_1
@@ -82,21 +87,29 @@ def __init__(
# Computed quantities (initially None)
self.Q_cross: Optional[float] = None
self.E_cross: Optional[float] = None
- self.coupling: Optional[float] = None
+ self.coupling: Optional[float] = coupling
self.reorganization_energy: Optional[float] = None
self.activation_energy: Optional[float] = None
self.transfer_rate: Optional[NDArray[np.float64]] = None
self.temperature: Optional[NDArray[np.float64]] = None
- def get_coupling(self, Q_cross: Optional[float] = None) -> float:
+ def get_coupling(
+ self,
+ H_ab: Optional[float] = None,
+ Q_cross: Optional[float] = None,
+ ) -> float:
"""
- Calculate electronic coupling between diabatic states.
+ Set the electronic coupling and locate the diabatic crossing.
- The coupling Hab is half the energy splitting between adiabatic
- states at the intersection point of the diabatic surfaces.
+ The coupling H_ab cannot be derived from the diabatic energies alone —
+ it must come from an electronic-structure calculation, e.g. half the
+ adiabatic splitting at the diabatic crossing.
Parameters
----------
+ H_ab : float, optional
+ Electronic coupling (eV). If None, uses the value passed to the
+ constructor. Raises if neither is set.
Q_cross : float, optional
Intersection point (amu^0.5·Å). If None, will find crossing
automatically using find_crossing().
@@ -104,7 +117,7 @@ def get_coupling(self, Q_cross: Optional[float] = None) -> float:
Returns
-------
coupling : float
- Electronic coupling Hab (eV)
+ Electronic coupling H_ab (eV)
Raises
------
@@ -112,6 +125,8 @@ def get_coupling(self, Q_cross: Optional[float] = None) -> float:
If potentials don't have fit functions
ValueError
If no crossing point found
+ ValueError
+ If no H_ab is supplied (constructor or argument)
Notes
-----
@@ -123,10 +138,10 @@ def get_coupling(self, Q_cross: Optional[float] = None) -> float:
Examples
--------
- >>> tc.get_coupling() # Automatic crossing detection
+ >>> tc.get_coupling(H_ab=0.015) # Automatic crossing detection
+ 0.015
+ >>> tc.get_coupling(H_ab=0.015, Q_cross=2.5) # Specify crossing point
0.015
- >>> tc.get_coupling(Q_cross=2.5) # Specify crossing point
- 0.018
"""
if self.pot_1.fit_func is None or self.pot_2.fit_func is None:
raise ValueError("Both potentials must be fitted before calculating coupling")
@@ -145,34 +160,23 @@ def get_coupling(self, Q_cross: Optional[float] = None) -> float:
E2 = self.pot_2(Q_cross)
E_cross = 0.5 * (E1 + E2)
- # Calculate adiabatic states
- # E+ = 0.5 * (E1 + E2) + sqrt(0.25 * (E1 - E2)^2 + Hab^2)
- # E- = 0.5 * (E1 + E2) - sqrt(0.25 * (E1 - E2)^2 + Hab^2)
-
- # At diabatic crossing, E1 ≈ E2, so:
- # E+ - E- = 2 * sqrt(Hab^2) = 2 * Hab
-
- # For a more general approach (not exactly at crossing):
- # We use the splitting to estimate Hab
- delta_E = abs(E1 - E2)
-
- # If exactly at crossing (delta_E ≈ 0), coupling is half the splitting
- # Otherwise, we assume a small coupling relative to delta_E
- if delta_E < 1e-6:
- # At crossing, assume minimal splitting (numerical precision limit)
- # This is a limitation - we can't measure Hab < ~1e-6 eV this way
- coupling = 1e-6 # Placeholder
- else:
- # Away from crossing, we can't determine Hab from energies alone
- # This method only works at the diabatic crossing
- # For now, use a rough estimate
- coupling = 0.5 * delta_E
+ if H_ab is None:
+ H_ab = self.coupling
+ if H_ab is None:
+ raise ValueError(
+ "Electronic coupling H_ab cannot be determined from the diabatic "
+ "energies alone. Supply it from an electronic-structure calculation "
+ "(e.g. half the adiabatic splitting at the diabatic crossing, "
+ "2|H_ab| = E+ - E-), via the constructor or get_coupling(H_ab=...)."
+ )
+ if H_ab <= 0:
+ raise ValueError(f"H_ab must be positive, got {H_ab}")
self.Q_cross = Q_cross
self.E_cross = E_cross
- self.coupling = coupling
+ self.coupling = H_ab
- return coupling
+ return H_ab
def get_reorganization_energy(self) -> float:
"""
@@ -478,6 +482,7 @@ def from_dict(cls, data: dict) -> "TransferCoordinate":
pot_1=pot_1,
pot_2=pot_2,
name=data.get("name", ""),
+ coupling=data.get("coupling"),
)
# Restore computed values
@@ -485,8 +490,6 @@ def from_dict(cls, data: dict) -> "TransferCoordinate":
tc.Q_cross = data["Q_cross"]
if "E_cross" in data:
tc.E_cross = data["E_cross"]
- if "coupling" in data:
- tc.coupling = data["coupling"]
if "reorganization_energy" in data:
tc.reorganization_energy = data["reorganization_energy"]
if "activation_energy" in data:
diff --git a/src/carriercapture/io/doped_interface.py b/src/carriercapture/io/doped_interface.py
index dab7d98..1bce787 100644
--- a/src/carriercapture/io/doped_interface.py
+++ b/src/carriercapture/io/doped_interface.py
@@ -99,7 +99,11 @@ def load_defect_entry(file_path: Union[str, Path]) -> Any:
def get_available_charge_states(defect_entry: Any) -> List[int]:
"""
- Get list of available charge states from DefectEntry.
+ Get the charge state of a single DefectEntry as a one-element list.
+
+ A doped DefectEntry represents one charge state; to work with several
+ charge states, load multiple DefectEntry files (or use
+ DefectThermodynamics).
Parameters
----------
@@ -109,13 +113,12 @@ def get_available_charge_states(defect_entry: Any) -> List[int]:
Returns
-------
List[int]
- Available charge states
+ Single-element list with this entry's charge state
Examples
--------
- >>> charges = get_available_charge_states(defect)
- >>> print(charges)
- [-2, -1, 0, +1, +2]
+ >>> get_available_charge_states(defect)
+ [-1]
"""
_check_doped_available()
@@ -895,7 +898,7 @@ def estimate_phonon_frequency(
displacement) should be used. This estimate is useful for initial
harmonic potential approximations.
"""
- from .._constants import HBAR
+ from .._constants import AMU, HBAR, HBAR_C
Q_data = np.asarray(Q_data)
E_data = np.asarray(E_data)
@@ -967,38 +970,10 @@ def estimate_phonon_frequency(
# Ensure curvature is positive
curvature = abs(curvature)
- # Convert curvature to phonon frequency
- # hw = hbar * sqrt(k / m_eff)
- # With Q in amu^0.5*A and E in eV:
- # k has units eV / (amu * A²)
- # m_eff = 1 amu (mass-weighted coordinate)
- #
- # hw (eV) = hbar (eV*s) * sqrt(k (eV/(amu*A²)) / m (amu))
- # Need to convert units properly
-
- # hbar in eV*fs, need angular frequency in rad/fs
- # k in eV/(amu*A²), convert A² to m² and amu to kg for SI
- # Actually simpler: use natural units
-
- # k in eV/(amu*A²) -> convert to eV/(eV/c² * A²) = c²/A²
- # omega² = k/m = k (eV/(amu*A²)) / (1 amu)
- # omega = sqrt(k) * (1/A) * sqrt(eV/amu)
-
- # Conversion: 1 amu = 931.5 MeV/c², 1 A = 1e-10 m
- # sqrt(eV/amu) = sqrt(1e-6 MeV / 931.5 MeV/c²) = sqrt(1.074e-9) c = 3.28e-5 c
- # In frequency: sqrt(eV/(amu*A²)) * A = sqrt(eV/amu) / A
-
- # Simpler approach using known conversion:
- # For harmonic oscillator: E_n = hw * (n + 0.5)
- # hw = hbar * omega = hbar * sqrt(k/m)
- #
- # With k in eV/A² and m in amu:
- # hw (eV) = 0.004136 * sqrt(k (eV/A²) / m (amu))
- # But our k is in eV/(amu*A²), so m_eff = 1:
- # hw (eV) = 0.004136 * sqrt(k (eV/(amu*A²)))
-
- conversion_factor = 0.004135665 # sqrt(hbar² / amu) in eV*A
- hw = conversion_factor * np.sqrt(curvature)
+ # Mass-weighted coordinates (Q in amu^0.5*A, E = 0.5*k*Q^2):
+ # hw = hbar*omega = sqrt(hbar^2 * k) with hbar^2 = (hbar*c)^2/AMU
+ # in eV*amu*A^2, so hw (eV) = (hbar*c / sqrt(AMU)) * sqrt(k)
+ hw = (HBAR_C * 1e10 / np.sqrt(AMU)) * np.sqrt(curvature)
# Angular frequency: omega = hw / hbar
omega = hw / HBAR # rad/s
@@ -1108,9 +1083,6 @@ def calculate_Q0_crossing(
if method == "crossing":
try:
Q0, E_crossing = find_crossing(pot_initial, pot_final)
- # Calculate barriers
- E_at_Q0_initial = pot_initial(Q0) if callable(pot_initial) else E_crossing
- E_at_Q0_final = pot_final(Q0) if callable(pot_final) else E_crossing
barrier_initial = E_crossing - E0_initial
barrier_final = E_crossing - E0_final
except (ValueError, RuntimeError):
@@ -1221,7 +1193,7 @@ def create_ccd_from_defect_entries(
nev_final : int, default=60
Number of eigenvalues to compute for final potential
W : float, optional
- Electron-phonon coupling matrix element (eV).
+ Electron-phonon coupling matrix element (eV/(amu^0.5·Å)).
If None, must be set later before calculating capture coefficient.
degeneracy : int, default=1
Degeneracy factor for the capture process
diff --git a/src/carriercapture/io/readers.py b/src/carriercapture/io/readers.py
index be36af3..4b92f52 100644
--- a/src/carriercapture/io/readers.py
+++ b/src/carriercapture/io/readers.py
@@ -14,6 +14,22 @@
import numpy as np
from numpy.typing import NDArray
+# Extension -> file format, shared by readers, writers, and the CLI
+FORMAT_MAP = {
+ '.json': 'json',
+ '.yaml': 'yaml',
+ '.yml': 'yaml',
+ '.npz': 'npz',
+ '.dat': 'dat',
+ '.txt': 'dat',
+ '.csv': 'csv',
+}
+
+
+def detect_format(filepath: Union[str, Path], default: str = 'json') -> str:
+ """Detect file format from the extension, falling back to `default`."""
+ return FORMAT_MAP.get(Path(filepath).suffix.lower(), default)
+
def read_potential_data(
filepath: Union[str, Path],
@@ -243,19 +259,8 @@ def load_potential_from_file(
if not filepath.exists():
raise FileNotFoundError(f"File not found: {filepath}")
- # Auto-detect format from extension
if file_format is None:
- ext = filepath.suffix.lower()
- format_map = {
- '.json': 'json',
- '.yaml': 'yaml',
- '.yml': 'yaml',
- '.npz': 'npz',
- '.dat': 'dat',
- '.txt': 'dat',
- '.csv': 'csv',
- }
- file_format = format_map.get(ext, 'dat')
+ file_format = detect_format(filepath, default='dat')
# Load based on format
if file_format == 'json':
diff --git a/src/carriercapture/io/writers.py b/src/carriercapture/io/writers.py
index 37080fb..601adfe 100644
--- a/src/carriercapture/io/writers.py
+++ b/src/carriercapture/io/writers.py
@@ -267,19 +267,9 @@ def save_potential(
"""
filepath = Path(filepath)
- # Auto-detect format from extension
if file_format is None:
- ext = filepath.suffix.lower()
- format_map = {
- '.json': 'json',
- '.yaml': 'yaml',
- '.yml': 'yaml',
- '.npz': 'npz',
- '.dat': 'dat',
- '.txt': 'dat',
- '.csv': 'csv',
- }
- file_format = format_map.get(ext, 'json')
+ from .readers import detect_format
+ file_format = detect_format(filepath)
# Save based on format
if file_format == 'json':
diff --git a/src/carriercapture/visualization/__init__.py b/src/carriercapture/visualization/__init__.py
index 7505d50..787ca99 100644
--- a/src/carriercapture/visualization/__init__.py
+++ b/src/carriercapture/visualization/__init__.py
@@ -6,6 +6,7 @@
plot_eigenvalue_spectrum,
plot_configuration_coordinate,
plot_overlap_matrix,
+ plot_scan_heatmap,
)
from .themes import (
@@ -30,6 +31,7 @@
"plot_eigenvalue_spectrum",
"plot_configuration_coordinate",
"plot_overlap_matrix",
+ "plot_scan_heatmap",
# Themes
"COLORS",
"POTENTIAL_COLORS",
diff --git a/src/carriercapture/visualization/interactive.py b/src/carriercapture/visualization/interactive.py
index 841f79c..19bb91a 100644
--- a/src/carriercapture/visualization/interactive.py
+++ b/src/carriercapture/visualization/interactive.py
@@ -381,6 +381,8 @@ def create_scan_tab(theme: Dict[str, Any]) -> html.Div:
dcc.Input(id="scan-hw-i", type="number", value=0.008, step=0.001, style=theme["input"]),
html.Label("ℏω_f (eV):", style=theme["text"]),
dcc.Input(id="scan-hw-f", type="number", value=0.008, step=0.001, style=theme["input"]),
+ html.Label("W (e-ph coupling, eV/(amu^0.5·Å)):", style=theme["text"]),
+ dcc.Input(id="scan-W", type="number", value=0.05, step=0.001, style=theme["input"]),
html.Label("Temperature (K):", style=theme["text"]),
dcc.Input(id="scan-temp", type="number", value=300, step=10, style=theme["input"]),
@@ -550,7 +552,7 @@ def create_capture_tab(theme: Dict[str, Any]) -> html.Div:
# Calculation parameters
html.H3("Parameters", style=theme["subheader"]),
- html.Label("W (e-ph coupling, eV):", style=theme["text"]),
+ html.Label("W (e-ph coupling, eV/(amu^0.5·Å)):", style=theme["text"]),
dcc.Input(id="capture-w", type="number", value=0.068, step=0.001, style=theme["input"]),
html.Label("Q₀ (crossing point):", style=theme["text"]),
dcc.Input(id="capture-q0", type="number", value=10.0, step=0.1, style=theme["input"]),
@@ -760,13 +762,14 @@ def register_scan_callbacks(app: dash.Dash) -> None:
State("scan-de-points", "value"),
State("scan-hw-i", "value"),
State("scan-hw-f", "value"),
+ State("scan-W", "value"),
State("scan-temp", "value"),
State("scan-results-store", "data")],
prevent_initial_call=True,
)
def handle_scan_operations(upload_contents, run_clicks, plot_type, plot_options,
filename, dq_min, dq_max, dq_points, de_min, de_max, de_points,
- hw_i, hw_f, temp, current_results):
+ hw_i, hw_f, W, temp, current_results):
"""Handle parameter scan operations."""
triggered_id = ctx.triggered_id
@@ -799,6 +802,7 @@ def handle_scan_operations(upload_contents, run_clicks, plot_type, plot_options,
dE_range=(de_min, de_max, de_points),
hbar_omega_i=hw_i,
hbar_omega_f=hw_f,
+ W=W,
temperature=temp,
)
@@ -1136,42 +1140,9 @@ def create_potential_figure(pot: Potential, display_options: List[str], wf_scale
def create_scan_figure(results: ScanResult, plot_type: str, log_scale: bool) -> go.Figure:
"""Create figure for scan results."""
- Z = results.capture_coefficients.copy()
-
- if log_scale:
- Z = np.log10(Z + 1e-30)
- colorbar_title = "log₁₀(C) [cm³/s]"
- else:
- colorbar_title = "C [cm³/s]"
-
- fig = go.Figure()
-
- if plot_type == "heatmap":
- fig.add_trace(go.Heatmap(
- x=results.dE_grid,
- y=results.dQ_grid,
- z=Z,
- colorscale="Viridis",
- colorbar=dict(title=colorbar_title),
- ))
- else: # contour
- fig.add_trace(go.Contour(
- x=results.dE_grid,
- y=results.dQ_grid,
- z=Z,
- colorscale="Viridis",
- colorbar=dict(title=colorbar_title),
- contours=dict(showlabels=True),
- ))
-
- fig.update_layout(
- title="Parameter Scan: Capture Coefficient",
- xaxis_title="ΔE (eV)",
- yaxis_title="ΔQ (amu^0.5·Å)",
- template="plotly_white",
- )
+ from .static import plot_scan_heatmap
- return fig
+ return plot_scan_heatmap(results, plot_type=plot_type, log_scale=log_scale)
def create_comparison_figure(potentials: List[Potential]) -> go.Figure:
diff --git a/src/carriercapture/visualization/static.py b/src/carriercapture/visualization/static.py
index adfbd4f..cae8717 100644
--- a/src/carriercapture/visualization/static.py
+++ b/src/carriercapture/visualization/static.py
@@ -14,6 +14,15 @@
import plotly.graph_objects as go
from plotly.subplots import make_subplots
+# Shared layout defaults for all static plots
+_LAYOUT_DEFAULTS = dict(
+ template="plotly_white",
+ font=dict(size=14),
+ hovermode="closest",
+ width=900,
+ height=600,
+)
+
def plot_potential(
potential,
@@ -181,11 +190,7 @@ def plot_potential(
title=title,
xaxis_title="Q (amu0.5·Å)",
yaxis_title="Energy (eV)",
- template="plotly_white",
- font=dict(size=14),
- hovermode="closest",
- width=900,
- height=600,
+ **_LAYOUT_DEFAULTS,
)
return fig
@@ -257,11 +262,7 @@ def plot_capture_coefficient(
title=title,
xaxis_title="1000/T (K-1)",
yaxis_title="log₁₀(C) [cm³/s]",
- template="plotly_white",
- font=dict(size=14),
- hovermode="closest",
- width=900,
- height=600,
+ **_LAYOUT_DEFAULTS,
)
# Add temperature labels on top axis if requested
@@ -383,10 +384,7 @@ def plot_eigenvalue_spectrum(
title=title,
xaxis_title="Quantum Number n",
yaxis_title="Energy (eV)",
- template="plotly_white",
- font=dict(size=14),
- width=900,
- height=600,
+ **_LAYOUT_DEFAULTS,
xaxis=dict(dtick=5),
)
@@ -469,10 +467,7 @@ def plot_configuration_coordinate(
title=title,
xaxis_title="Q (amu0.5·Å)",
yaxis_title="Energy (eV)",
- template="plotly_white",
- font=dict(size=14),
- width=900,
- height=600,
+ **_LAYOUT_DEFAULTS,
)
return fig
@@ -538,10 +533,72 @@ def plot_overlap_matrix(
title=title,
xaxis_title="Final State j",
yaxis_title="Initial State i",
- template="plotly_white",
- font=dict(size=14),
- width=800,
- height=700,
+ **{**_LAYOUT_DEFAULTS, "width": 800, "height": 700},
+ )
+
+ return fig
+
+
+def plot_scan_heatmap(
+ results,
+ plot_type: str = "heatmap",
+ log_scale: bool = True,
+ title: str = "Parameter Scan: Capture Coefficient",
+ colorscale: str = "Viridis",
+ **layout_kwargs,
+) -> go.Figure:
+ """
+ Plot a parameter-scan capture coefficient map.
+
+ Parameters
+ ----------
+ results : ScanResult
+ Scan results from ParameterScanner
+ plot_type : str, default="heatmap"
+ "heatmap", "contour", or "both"
+ log_scale : bool, default=True
+ Plot log10(C) instead of C
+ title : str
+ Plot title
+ colorscale : str, default="Viridis"
+ Plotly colorscale name
+ **layout_kwargs
+ Layout overrides (e.g. width, height)
+
+ Returns
+ -------
+ fig : go.Figure
+ Plotly figure
+ """
+ if plot_type not in ("heatmap", "contour", "both"):
+ raise ValueError(f"Unknown plot_type: {plot_type}. Use 'heatmap', 'contour', or 'both'")
+
+ Z = np.asarray(results.capture_coefficients, dtype=float)
+ if log_scale:
+ Z = np.log10(Z + 1e-30)
+ colorbar_title = "log₁₀(C) [cm³/s]"
+ else:
+ colorbar_title = "C [cm³/s]"
+
+ trace_kwargs = dict(
+ x=results.dE_grid,
+ y=results.dQ_grid,
+ z=Z,
+ colorscale=colorscale,
+ colorbar=dict(title=colorbar_title),
+ )
+
+ fig = go.Figure()
+ if plot_type in ("heatmap", "both"):
+ fig.add_trace(go.Heatmap(**trace_kwargs))
+ if plot_type in ("contour", "both"):
+ fig.add_trace(go.Contour(contours=dict(showlabels=True), **trace_kwargs))
+
+ fig.update_layout(
+ title=title,
+ xaxis_title="ΔE (eV)",
+ yaxis_title="ΔQ (amu0.5·Å)",
+ **{**_LAYOUT_DEFAULTS, "width": 800, "height": 700, **layout_kwargs},
)
return fig
@@ -553,4 +610,5 @@ def plot_overlap_matrix(
"plot_eigenvalue_spectrum",
"plot_configuration_coordinate",
"plot_overlap_matrix",
+ "plot_scan_heatmap",
]
diff --git a/tests/test_doped_integration.py b/tests/test_doped_integration.py
index cb17525..bfde5c9 100644
--- a/tests/test_doped_integration.py
+++ b/tests/test_doped_integration.py
@@ -347,6 +347,20 @@ def test_estimate_phonon_frequency_curvature_method(self):
# Curvature should be close to k
assert np.isclose(result['curvature'], k, rtol=0.1)
+ @pytest.mark.skipif(not DOPED_INTEGRATION_AVAILABLE, reason="doped package not installed")
+ def test_estimate_phonon_frequency_roundtrip(self):
+ """A harmonic PES with known hw must round-trip exactly."""
+ from carriercapture._constants import AMU, HBAR_C
+
+ hw_true = 0.008 # eV
+ k = AMU * (hw_true / (HBAR_C * 1e10)) ** 2 # E = 0.5*k*Q^2
+ Q_data = np.linspace(-5, 5, 201)
+ E_data = 0.5 * k * Q_data**2
+
+ result = estimate_phonon_frequency(Q_data, E_data, method="curvature")
+
+ assert np.isclose(result['hw'], hw_true, rtol=1e-6)
+
@pytest.mark.skipif(not DOPED_INTEGRATION_AVAILABLE, reason="doped package not installed")
def test_estimate_phonon_frequency_harmonic_fit_method(self):
"""Test phonon frequency estimation using harmonic_fit method."""
diff --git a/tests/test_parameter_scan.py b/tests/test_parameter_scan.py
index f3c28da..d909dd1 100644
--- a/tests/test_parameter_scan.py
+++ b/tests/test_parameter_scan.py
@@ -118,12 +118,19 @@ def simple_params(self):
dE_range=(0, 0.5, 3),
hbar_omega_i=0.008,
hbar_omega_f=0.008,
+ W=0.05,
temperature=300.0,
nev_initial=60, # Need enough for partition function convergence at 300K
nev_final=30, # Scaled proportionally
Q_grid_points=500, # Reduced for speed
)
+ def test_scanner_requires_W(self):
+ """ParameterScanner rejects parameters without a coupling W."""
+ params = ScanParameters(dQ_range=(0, 5, 3), dE_range=(0, 0.5, 3))
+ with pytest.raises(ValueError, match="W"):
+ ParameterScanner(params)
+
def test_create_scanner(self, simple_params):
"""Test creating ParameterScanner."""
scanner = ParameterScanner(simple_params, verbose=False)
@@ -165,19 +172,6 @@ def test_create_harmonic_potentials(self, simple_params):
assert len(pot_i.eigenvalues) == 60 # nev_initial
assert len(pot_f.eigenvalues) == 30 # nev_final
- def test_calculate_W_coupling(self, simple_params):
- """Test W coupling calculation."""
- scanner = ParameterScanner(simple_params, verbose=False)
-
- W = scanner._calculate_W_coupling(
- hbar_omega_f=0.008,
- dQ=2.0,
- dE=0.5
- )
-
- assert isinstance(W, float)
- assert W > 0
-
def test_calculate_single_point(self, simple_params):
"""Test single point calculation."""
scanner = ParameterScanner(simple_params, verbose=False)
diff --git a/tests/test_sommerfeld.py b/tests/test_sommerfeld.py
new file mode 100644
index 0000000..19f24ab
--- /dev/null
+++ b/tests/test_sommerfeld.py
@@ -0,0 +1,77 @@
+"""
+Tests for the Sommerfeld (Coulomb enhancement) factor.
+"""
+
+import pytest
+import numpy as np
+
+from carriercapture import sommerfeld_parameter
+from carriercapture._constants import K_B, RYDBERG_EV
+
+
+class TestSommerfeldParameter:
+ """Test sommerfeld_parameter behavior and limits."""
+
+ def test_neutral_is_unity(self):
+ """Z=0 gives exactly 1 for both methods, scalar and array."""
+ for method in ("Analytic", "Integrate"):
+ assert sommerfeld_parameter(300.0, Z=0, m_eff=0.2, eps0=10, method=method) == 1.0
+ T = np.linspace(100, 500, 5)
+ s = sommerfeld_parameter(T, Z=0, m_eff=0.2, eps0=10, method=method)
+ np.testing.assert_array_equal(s, np.ones(5))
+
+ def test_attractive_enhances_and_decreases_with_T(self):
+ """Attractive center: s > 1 and monotonically decreasing with T."""
+ T = np.linspace(50, 800, 20)
+ for method in ("Analytic", "Integrate"):
+ s = sommerfeld_parameter(T, Z=-1, m_eff=0.2, eps0=10, method=method)
+ assert np.all(s > 1)
+ assert np.all(np.diff(s) < 0)
+
+ def test_repulsive_suppresses_and_increases_with_T(self):
+ """Repulsive center: 0 < s < 1 and increasing with T."""
+ T = np.linspace(50, 800, 20)
+ for method in ("Analytic", "Integrate"):
+ s = sommerfeld_parameter(T, Z=+1, m_eff=0.2, eps0=10, method=method)
+ assert np.all(s > 0)
+ assert np.all(s < 1)
+ assert np.all(np.diff(s) > 0)
+
+ def test_analytic_attractive_regression(self):
+ """Hand-derived analytic value: s = 4*sqrt(theta/pi)."""
+ T, Z, m_eff, eps0 = 300.0, -1, 0.2, 10.0
+ E_R = m_eff * RYDBERG_EV / eps0**2
+ theta = np.pi**2 * Z**2 * E_R / (K_B * T)
+ expected = 4.0 * np.sqrt(theta / np.pi)
+ s = sommerfeld_parameter(T, Z=Z, m_eff=m_eff, eps0=eps0, method="Analytic")
+ assert s == pytest.approx(expected, rel=1e-12)
+ assert s == pytest.approx(7.27, rel=1e-2)
+
+ def test_analytic_matches_integrate_at_low_T(self):
+ """The Pässler analytic form is the low-T limit of the integral."""
+ s_a = sommerfeld_parameter(50.0, Z=-1, m_eff=0.2, eps0=10, method="Analytic")
+ s_i = sommerfeld_parameter(50.0, Z=-1, m_eff=0.2, eps0=10, method="Integrate")
+ assert s_a == pytest.approx(s_i, rel=1e-3)
+ # The exact enhancement eta/(1 - exp(-eta)) >= eta, so the thermal
+ # average always sits at or above the analytic (low-T) limit, with
+ # the gap growing at high T.
+ s_a_hot = sommerfeld_parameter(800.0, Z=-1, m_eff=0.2, eps0=10, method="Analytic")
+ s_i_hot = sommerfeld_parameter(800.0, Z=-1, m_eff=0.2, eps0=10, method="Integrate")
+ assert s_i_hot > s_a_hot
+ assert (s_i_hot - s_a_hot) / s_a_hot > (s_i - s_a) / s_a
+
+ def test_shapes_and_types(self):
+ """Scalar in -> float out; array in -> same-shape array out."""
+ s = sommerfeld_parameter(300.0, Z=-1, m_eff=0.2, eps0=10)
+ assert np.isscalar(s) or s.ndim == 0
+ T = np.linspace(100, 500, 7)
+ s_arr = sommerfeld_parameter(T, Z=-1, m_eff=0.2, eps0=10)
+ assert s_arr.shape == T.shape
+
+ def test_unknown_method_raises(self):
+ with pytest.raises(ValueError, match="method"):
+ sommerfeld_parameter(300.0, Z=-1, m_eff=0.2, eps0=10, method="bogus")
+
+
+if __name__ == "__main__":
+ pytest.main([__file__, "-v"])
diff --git a/tests/test_transfer_coord.py b/tests/test_transfer_coord.py
index 36e976b..e9074e5 100644
--- a/tests/test_transfer_coord.py
+++ b/tests/test_transfer_coord.py
@@ -2,17 +2,11 @@
Tests for TransferCoordinate class.
Validates Marcus theory calculations for charge transfer rates and mobility.
-
-NOTE: TransferCoordinate is a Phase 3 feature that is not fully implemented yet.
-These tests are skipped to allow CI to pass while development continues.
"""
import pytest
import numpy as np
-# Skip entire module - Phase 3 transfer coordinate features not complete
-pytestmark = pytest.mark.skip(reason="Phase 3 transfer coordinate features not fully implemented")
-
from carriercapture.core import Potential, TransferCoordinate
@@ -49,11 +43,11 @@ def test_calculate_coupling_basic(self):
tc = TransferCoordinate(pot_1, pot_2)
- # Calculate coupling (auto-detect crossing)
- coupling = tc.get_coupling()
+ # Set coupling (auto-detect crossing)
+ coupling = tc.get_coupling(H_ab=0.01)
- assert coupling is not None
- assert coupling > 0
+ assert coupling == 0.01
+ assert tc.coupling == 0.01
assert tc.Q_cross is not None
assert tc.E_cross is not None
@@ -69,9 +63,9 @@ def test_calculate_coupling_specified_crossing(self):
# Specify crossing point manually
Q_cross_manual = 5.0
- coupling = tc.get_coupling(Q_cross=Q_cross_manual)
+ coupling = tc.get_coupling(H_ab=0.01, Q_cross=Q_cross_manual)
- assert coupling > 0
+ assert coupling == 0.01
assert tc.Q_cross == Q_cross_manual
# Energy should be evaluated at specified point
@@ -87,8 +81,37 @@ def test_coupling_without_fit_raises(self):
tc = TransferCoordinate(pot_1, pot_2)
with pytest.raises(ValueError, match="must be fitted"):
+ tc.get_coupling(H_ab=0.01)
+
+ def test_coupling_without_hab_raises(self):
+ """Test that get_coupling requires a user-supplied H_ab."""
+ pot_1 = Potential.from_harmonic(hw=0.03, Q0=0.0, E0=1.0)
+ pot_2 = Potential.from_harmonic(hw=0.03, Q0=8.0, E0=0.0)
+
+ tc = TransferCoordinate(pot_1, pot_2)
+
+ with pytest.raises(ValueError, match="H_ab"):
tc.get_coupling()
+ def test_coupling_from_constructor(self):
+ """Test that constructor-supplied coupling is used."""
+ pot_1 = Potential.from_harmonic(hw=0.03, Q0=0.0, E0=1.0)
+ pot_2 = Potential.from_harmonic(hw=0.03, Q0=8.0, E0=0.0)
+
+ tc = TransferCoordinate(pot_1, pot_2, coupling=0.02)
+
+ assert tc.get_coupling() == 0.02
+
+ def test_coupling_nonpositive_raises(self):
+ """Test that non-positive H_ab is rejected."""
+ pot_1 = Potential.from_harmonic(hw=0.03, Q0=0.0, E0=1.0)
+ pot_2 = Potential.from_harmonic(hw=0.03, Q0=8.0, E0=0.0)
+
+ tc = TransferCoordinate(pot_1, pot_2)
+
+ with pytest.raises(ValueError, match="positive"):
+ tc.get_coupling(H_ab=-0.01)
+
class TestReorganizationEnergy:
"""Test reorganization energy calculation."""
@@ -106,13 +129,12 @@ def test_reorganization_energy_symmetric(self):
lambda_reorg = tc.get_reorganization_energy()
- # For harmonic potentials: λ = 0.5 * k * dQ^2 = 0.5 * m * ω^2 * dQ^2
- # But our harmonic uses: E = a * (Q - Q0)^2 where a = (amu/2) * (hw/(hbar_c*1e10))^2
- # So λ = a * dQ^2
+ # For our harmonic form E = a*(Q - Q0)^2 with a = (AMU/2)*(hw/(HBAR_C*1e10))^2,
+ # the reorganization energy is exactly λ = a * dQ^2
+ from carriercapture._constants import AMU, HBAR_C
- # Let's just check it's positive and reasonable
- assert lambda_reorg > 0
- assert lambda_reorg < 1.0 # Should be less than 1 eV for typical parameters
+ a = 0.5 * AMU * (hw / (HBAR_C * 1e10)) ** 2
+ assert lambda_reorg == pytest.approx(a * dQ**2, rel=1e-6)
def test_reorganization_energy_asymmetric(self):
"""Test reorganization energy for asymmetric case."""
@@ -229,7 +251,7 @@ def test_transfer_rate_basic(self):
tc = TransferCoordinate(pot_1, pot_2)
# Calculate prerequisites
- tc.get_coupling()
+ tc.get_coupling(H_ab=0.01)
tc.get_reorganization_energy()
# Calculate transfer rate
@@ -248,7 +270,7 @@ def test_transfer_rate_temperature_dependence(self):
pot_2 = Potential.from_harmonic(hw=0.02, Q0=10.0, E0=0.0)
tc = TransferCoordinate(pot_1, pot_2)
- tc.get_coupling()
+ tc.get_coupling(H_ab=0.01)
tc.get_reorganization_energy()
temperature = np.array([100.0, 200.0, 300.0, 400.0, 500.0])
@@ -269,7 +291,7 @@ def test_transfer_rate_downhill(self):
pot_2 = Potential.from_harmonic(hw=0.02, Q0=8.0, E0=0.0)
tc = TransferCoordinate(pot_1, pot_2)
- tc.get_coupling()
+ tc.get_coupling(H_ab=0.01)
tc.get_reorganization_energy()
temperature = np.array([300.0])
@@ -297,7 +319,7 @@ def test_transfer_rate_without_prerequisites_raises(self):
tc.get_transfer_rate(temperature=temperature)
# With coupling but without reorganization energy
- tc.get_coupling()
+ tc.get_coupling(H_ab=0.01)
with pytest.raises(ValueError, match="reorganization energy"):
tc.get_transfer_rate(temperature=temperature)
@@ -311,7 +333,7 @@ def test_mobility_basic(self):
pot_2 = Potential.from_harmonic(hw=0.02, Q0=8.0, E0=0.0)
tc = TransferCoordinate(pot_1, pot_2)
- tc.get_coupling()
+ tc.get_coupling(H_ab=0.01)
tc.get_reorganization_energy()
temperature = np.linspace(100, 500, 20)
@@ -331,7 +353,7 @@ def test_mobility_temperature_dependence(self):
pot_2 = Potential.from_harmonic(hw=0.02, Q0=10.0, E0=0.0)
tc = TransferCoordinate(pot_1, pot_2)
- tc.get_coupling()
+ tc.get_coupling(H_ab=0.01)
tc.get_reorganization_energy()
temperature = np.array([100.0, 200.0, 300.0, 400.0, 500.0])
@@ -349,7 +371,7 @@ def test_mobility_distance_dependence(self):
pot_2 = Potential.from_harmonic(hw=0.02, Q0=8.0, E0=0.0)
tc = TransferCoordinate(pot_1, pot_2)
- tc.get_coupling()
+ tc.get_coupling(H_ab=0.01)
tc.get_reorganization_energy()
temperature = np.array([300.0])
@@ -384,7 +406,7 @@ def test_to_dict_from_dict(self):
tc = TransferCoordinate(pot_1, pot_2, name="test_transfer")
# Calculate some properties
- tc.get_coupling()
+ tc.get_coupling(H_ab=0.01)
tc.get_reorganization_energy()
tc.get_activation_energy(delta_G=0.0)
tc.get_transfer_rate(temperature=np.linspace(100, 500, 20))
@@ -427,7 +449,7 @@ def test_simple_charge_transfer(self):
tc = TransferCoordinate(pot_1, pot_2, name="hole_transfer")
# Calculate properties in order
- coupling = tc.get_coupling()
+ coupling = tc.get_coupling(H_ab=0.01)
lambda_reorg = tc.get_reorganization_energy()
barrier = tc.get_activation_energy(delta_G=dE)
@@ -450,8 +472,11 @@ def test_simple_charge_transfer(self):
assert np.all(np.isfinite(mobility))
# Physical reasonableness
+ from carriercapture._constants import AMU, HBAR_C
+
+ a = 0.5 * AMU * (hw / (HBAR_C * 1e10)) ** 2
assert coupling < 0.5 # Typical coupling is << 1 eV
- assert lambda_reorg < 2.0 # Typical reorganization < few eV
+ assert lambda_reorg == pytest.approx(a * dQ**2, rel=1e-6) # λ = a·dQ² exactly
assert barrier < lambda_reorg # Barrier should be < λ
assert np.max(rate) < 1e15 # Transfer rate < phonon frequency
@@ -462,21 +487,21 @@ def test_symmetric_vs_asymmetric_transfer(self):
pot_2_sym = Potential.from_harmonic(hw=0.02, Q0=8.0, E0=0.0)
tc_sym = TransferCoordinate(pot_1_sym, pot_2_sym)
- tc_sym.get_coupling()
+ tc_sym.get_coupling(H_ab=0.01)
tc_sym.get_reorganization_energy()
temperature = np.array([300.0])
rate_sym = tc_sym.get_transfer_rate(temperature=temperature, delta_G=0.0)
- # Asymmetric case (downhill)
+ # Asymmetric case (downhill: initial state 0.3 eV above final, ΔG = -0.3)
pot_1_asym = Potential.from_harmonic(hw=0.02, Q0=0.0, E0=0.3)
pot_2_asym = Potential.from_harmonic(hw=0.02, Q0=8.0, E0=0.0)
tc_asym = TransferCoordinate(pot_1_asym, pot_2_asym)
- tc_asym.get_coupling()
+ tc_asym.get_coupling(H_ab=0.01)
tc_asym.get_reorganization_energy()
- rate_asym = tc_asym.get_transfer_rate(temperature=temperature, delta_G=0.3)
+ rate_asym = tc_asym.get_transfer_rate(temperature=temperature, delta_G=-0.3)
# Downhill should be faster
assert rate_asym[0] > rate_sym[0]