getelec.transmission_solver
Numerical solvers for quantum transmission probabilities.
Provides the TransmissionSolver interface and its numerical
implementations:
Noumerov/NoumerovFast-- direct integration of the 1D Schrodinger equation, with no approximation beyond discretisation, whose local truncation error is O(h^6).NoumerovReference-- the same integration written out plainly, one energy at a time, for checking the solvers above and debugging a result.NeuralSolver-- a trained network standing in for the Noumerov result, worth it for barriers with several parameters. Train one withgetelec.training.
The closed-form and semiclassical solutions (WKB, the Airy-function solution)
live in getelec.transmission_solutions. Every solver takes an array of
energies and returns an array of transmission probabilities, so they are
interchangeable inside an emitter.
Abstract base class for all transmission solvers.
Subclasses must implement calculate_transmission(). Implementing
calculate_transmission_batch() as well is optional but lets sweeps and
multi-band calculations be evaluated in a single pass.
Transmission probability at each energy.
Parameters
potential : potential_barrier.Barrier Barrier defining the spatial potential profile. energies : np.ndarray 1D array of energies, in eV.
Returns
np.ndarray
Transmission probabilities in [0, 1], same shape as energies.
Transmission for several (potential, energies) pairs at once.
The base implementation just loops. Solvers that can genuinely batch override this. Callers should prefer it over writing their own loop, so they benefit automatically when a solver does support batching.
Parameters
potentials : sequence of potential_barrier.Barrier energy_arrays : sequence of np.ndarray
Returns
list of np.ndarray
Noumerov integration that keeps the whole wavefunction.
Same physics and same settings as NoumerovFast, and it returns the
same transmission coefficients -- but it retains psi(x) at every grid point,
so the wavefunction is available afterwards through calculate_psi().
Use it when the wavefunction is the answer: a charge density for a Poisson
solve, a probability current, or simply looking at how the wave decays
through the barrier. Use NoumerovFast when only the transmission
matters.
The cost is memory and speed. The array is n_energy x n_grid complex,
16 bytes an element, so 800 energies on a 14000-point grid is 180 MB; and
the loop order cannot be inverted the way the endpoint kernel does, so it is
several times slower per energy. Ask for the energies you actually need.
Examples
>>> solver = Noumerov() # doctest: +SKIP
>>> x, psi = solver.calculate_psi(barrier, energies) # doctest: +SKIP
>>> density = np.abs(psi) ** 2 # doctest: +SKIP
The wavefunction on the spatial grid.
Parameters
potential : potential_barrier.Barrier energies : array_like Energies in eV. Each gets its own row. x_points : array_like, optional Positions at which to return psi. Defaults to the solver's own grid. Supplying a coarser set does not make the integration cheaper -- the recurrence has to step through every grid point regardless -- but it does cut the memory that comes back.
Returns
x, psi : np.ndarray
x has shape (n_x,), ordered from vacuum into the metal as the
integration runs. psi has shape (n_energy, n_x), complex, and is
normalised to unit outgoing flux in the vacuum.
Notes
The seed is imposed at the vacuum end and the recurrence runs inward, so
psi grows by many orders of magnitude across a thick barrier. That is
the physical behaviour of the growing solution, not an instability, but
it does mean |psi| inside the metal can be enormous; divide by the
incident amplitude if you want a normalised scattering state.
Probability current at every grid point, for each energy.
For a stationary scattering state this must be independent of position, so its variation across the grid is a direct check on the integration -- more informative than any residual, because it tests a conservation law the scheme does not enforce.
Returns
x, current : np.ndarray
Noumerov integration of the 1D Schrodinger equation, endpoints only.
Keeps only psi at the last two grid points, which is all the transmission
coefficient needs. That is what makes it fast and what bounds its memory to
O(n_grid) instead of O(n_energy x n_grid). If you need the wavefunction
itself -- a charge density, a probability current, the decaying tail inside
the barrier -- use Noumerov, which keeps all of it.
Integrates right-to-left from deep vacuum into the metal, seeding an outgoing plane wave on the vacuum side and projecting the result onto incoming and outgoing waves inside the metal. The transmission probability follows from the flux ratio.
Noumerov's method has a local truncation error of O(h^6), which makes it one of the most accurate methods for this equation, and it needs only a single three-term recurrence per grid point, which makes it fast.
Parameters
x_metal : float, default -0.01
Left-hand end of the integration domain, inside the metal (nm). The
potential is identically zero there and the wave is matched at the last
two grid points, so this only has to be negative: measured over the
planar, sharp-tip and triangular barriers at 0.2-12 V/nm, the current
density from x_metal=-0.01 and from -1.0 differ by 4e-10.
x_vac_plus : float, default 3.0
Distance past the barrier at which to start the integration (nm). The
seed is only exact for a flat potential, so this has to be large enough
that the residual gradient does not matter. With the default WKB seed,
3 nm changes the current density by 0.0024%, the Nottingham heat by
0.0034% and the transmission itself by 0.025% against a 20 nm domain,
measured over the same barriers and fields. Check a particular case
with calculate_convergence_report() rather than assuming.
h : float, default 1e-3
Grid spacing (nm). The local truncation error is O(h^6).
max_barrier_width : float, default 3.0
Expected barrier width (nm); together with x_vac_plus this sets the
right-hand end of the domain.
Notes
The spatial grid and the potential profile are cached and rebuilt only when the geometry or the barrier parameters actually change, so repeated calls at a fixed field cost almost nothing in setup.
Preset tuned for speed, with the accuracy trade-offs made explicit.
Three changes relative to the defaults, each justified by a convergence study:
x_vac_plus=2.0, one nanometre less vacuum than the default. Worth 7% of the grid, at a current density 0.0072% from a 20 nm domain instead of 0.0024%, and a transmission 0.066% from it instead of 0.025%.h=2e-3, half the grid points. Against h = 1.25e-4 this costs 0.015% in current density and 0.035% in the transmission itself, measured for the planar and sharp-tip barriers over work functions of 2.5 to 6 eV, fields of 0.5 to 12 V/nm and 300 to 3000 K. The triangular barrier, the one with a jump at the surface, costs 0.01% on the aligned grid (see_grid()).energy_nodes=48. Transmission is solved on a reduced energy set and log-interpolated onto the full grid.
The short metal side, which this preset used to carry alone, is now the default for every Noumerov solve.
All three together, against a converged reference (h = 1.25e-4 on a 20 nm domain) at 1 to 12 V/nm and 300 to 1500 K: current density within 5e-5, energy distributions within 2e-3, pointwise transmission within 7e-4. The default settings reach 2e-5, 2e-4 and 2e-4 on the same cases, so the preset costs a factor of a few in accuracy for a factor of three in time -- and both are a long way inside the 1% that separates transmission algorithms from each other.
Use the default constructor when you want the transmission itself to be converged pointwise across the whole band.
The reference implementation, NoumerovReference.
The same integration written out plainly, one energy at a time, to
check a result rather than to produce many: a metal current density
takes about half a second. Keyword arguments (x_start, x_end,
h, seed) go to NoumerovReference.
Outermost classical turning point over the given barriers and energies.
The Schottky-Nordheim barrier extends to x2 ~ (E_F + phi - E) / F, which for a low field and a deep energy is tens of nm -- far beyond the 3 nm default. Integrating only part of a barrier does not produce a slightly wrong answer, it produces a meaningless one, so the domain is measured rather than assumed.
Returns
float
Required barrier width in nm, or None if it cannot be determined.
Transmission probability at each energy for a single barrier.
Solve several barriers in one parallel launch.
Used by field sweeps and by the semiconductor emitter, which needs four energy grids against the same barrier. Batching amortises the kernel launch and gives the thread pool enough work to saturate.
Natural log of the transmission probability.
Preferred over log(calculate_transmission(...)) whenever the answer
may be very small: this reaches ln T of about -1200, where taking the
log of the returned probability bottoms out near -700 and then returns
-inf. Useful for deep tunnelling and for building interpolation tables.
Returns
np.ndarray ln T, with -inf for channels that carry no flux.
Refine the grid and report how much the answer moves.
The default h is not a guarantee of accuracy at every field and
energy range. This reruns the calculation at h, h/2, h/4 and
returns the relative change, so a claimed digit count can be checked
rather than assumed.
Returns
dict Maps each step size to the largest relative deviation from the finest grid.
Noumerov integration written out plainly: the reference implementation.
Solves each energy on its own with calculate_noumerov_reference(),
from full arrays, on a fixed grid. Nothing is cached, batched, interpolated
or reordered, and the kernel is compiled without fastmath. That makes it
slow -- about half a second for a metal current density, where
Noumerov takes a few milliseconds -- and easy to follow, which is
what it is for: checking the fast solvers and debugging a result. Select it
with reference=True in the one-line API or with
Noumerov.reference().
For the planar barrier (E_F = 7.5 eV, phi = 4.5 eV, 300 K, 3-7 V/nm) the
current density and Nottingham heat agree with Noumerov to 3e-5
with either seed, and the transmission to 5e-4 with the plane-wave seed and
3e-4 with the WKB seed.
The grid is np.arange(x_end, x_start - h, -h), fixed rather than sized
to the barrier, and two things follow from that:
- It has to enclose the barrier. Where
x_endis still inside it, V(x_end) >= E, the seed is not a travelling wave and the transmission comes out as 0. The solver warns when that happens. - The defaults put a node on the surface. A continuous barrier does not
notice, but
~getelec.potential_barrier.TriangularPotentialjumps there, which costs about 1% in current density.x_end=20.0005puts the surface halfway between two nodes and brings it within 2e-5 of the exact solution.
Energies at or below zero, the bottom of the band, have no travelling state in the metal, so no incident flux: they return 0 without being solved.
Parameters
x_start : float, default -1.0
End of the integration inside the metal (nm).
x_end : float, default 20.0
Start of the integration in the vacuum (nm).
h : float, default 1e-3
Grid spacing (nm).
seed : {'plane', 'wkb'}, default 'plane'
The outgoing wave imposed at x_end. See
calculate_noumerov_reference().
Examples
>>> getelec.current_density(field=5.0, reference=True) # doctest: +SKIP
>>> solver = NoumerovReference(x_end=20.0005, seed="wkb") # doctest: +SKIP
>>> T, x, V, psi, k_metal = calculate_noumerov_reference(barrier, 7.5) # doctest: +SKIP
Transmission probability at each energy, solved one energy at a time.
Energies at or below zero return 0 without being solved. Warns when
x_end lies inside the barrier for any of the others.
Inherited Members
Transmission at one energy by Noumerov integration, with everything behind it.
The computation of NoumerovReference, written to be read line by
line. The grid runs from x_end in the vacuum to x_start in the
metal; an outgoing wave is seeded at its first two points; the Schrodinger
equation is integrated into the metal by
getelec._kernels.run_noumerov_integration(); and psi at the last two
points is matched to A exp(ikx) + B exp(-ikx). The transmission is the
outgoing flux of the seed divided by the incident flux, k |A|^2.
Parameters
potential : potential_barrier.Barrier
The barrier, carrying its own parameters.
electron_energy : float
Energy in eV, measured from the bottom of the band, where V = 0.
x_start : float, default -1.0
End of the integration inside the metal (nm).
x_end : float, default 20.0
Start of the integration in the vacuum (nm). It has to lie past the
barrier: where V(x_end) >= E the seed is not a travelling wave and the
transmission comes out as 0.
h : float, default 1e-3
Grid spacing (nm).
seed : {'plane', 'wkb'}, default 'plane'
The outgoing wave imposed at x_end: a plane wave with the local
wavevector, or the WKB form k^(-1/2) exp(i int k dx) that
NoumerovFast uses by default.
Returns
T_current : float
Transmission probability.
x_points : np.ndarray
The grid, from the vacuum into the metal (nm).
V : np.ndarray
The potential on the grid (eV).
psi : np.ndarray
The wavefunction on the grid, complex: a unit-amplitude plane wave, or
a unit-flux WKB wave, at x_end.
k_metal : float
Wavevector inside the metal (1/nm).
Notes
Only energies above the bottom of the band have a travelling state in the
metal. Below it this returns nan, and at exactly zero the matching matrix
is singular and numpy.linalg.solve() raises;
NoumerovReference returns 0 there without calling it.
Transmission from a trained neural network.
The network never predicts D. It predicts the residual
R = ln D_noumerov - ln D_kemble
against a semiclassical reference for the same barrier. The reference
already carries the exponential -- hundreds of e-folds -- so the network
only learns a correction of order unity, and the error budget is explicit:
relative error in D is exactly absolute error in R.
When it is worth it. For the planar Schottky barrier the network is
several times faster than Noumerov.fast(), but that calculation was
already cheap: the model there is the worked example, not the reason. Its cost does not grow with the number
of barrier parameters, while exact solves or a table of them do, so it pays
off for barriers with several -- tip radius and field enhancement, as in the
shipped curved-tip model. The shipped Schottky model is there as a running
example of the workflow; getelec.training trains one for your own
barrier.
Parameters
model : None, str, Path or NeuralModel, optional
None picks the shipped model for the barrier's class (see
SHIPPED_MODELS). A path or a NeuralModel uses that
model, for barriers of the class it was trained on.
fallback : {'exact', 'wkb', 'error'}, default 'exact'
What to do for a barrier with no model, or outside the trained domain.
'exact' solves it with Noumerov, keeping the answer right at the
cost of speed. 'wkb' uses the numerical semiclassical result.
'error' raises, which is what a pipeline that must not silently
change accuracy wants.
exact_solver : optional
Used by the 'exact' fallback. Defaults to Noumerov.fast().
energy_nodes : int, default 48
Evaluate on this many energies and spline ln D onto the rest. Both
the reference and the forward pass cost time proportional to the number
of energies, and ln D is smooth in energy, so this makes the cost
flat in grid size and changes the answer by well under 1%. Set to 0 to
evaluate at every energy.
Examples
>>> solver = NeuralSolver() # doctest: +SKIP
>>> barrier = SmallRadiiPotential(7.5, 4.5, 5.0, radius=50.0, gamma=100.0)
>>> solver.calculate_transmission(barrier, energies) # doctest: +SKIP
Whether a model exists for this barrier and its parameters were trained on.
Energies whose barrier height falls inside the trained range.
Transmission probability at each energy.
Inherited Members
A trained network and everything needed to use it, stored as plain arrays.
save() writes an .npz that loads with allow_pickle=False: the
weights and biases, the names of the feature function, the semiclassical
reference and the barrier class, and the domain the model was trained on.
No pickle means no scikit-learn at load time and no coupling to the version
that trained it.
Parameters
weights, biases : list of np.ndarray
Layer by layer. Hidden layers use tanh, the output is linear.
features : str
Key into FEATURES.
reference : str
Key into REFERENCES.
barrier : str
Class name of the barrier the model was trained on.
domain : dict
name -> (lo, hi) for barrier_height, total_height, field
and each entry of parameters.
parameters : sequence of str
Extra barrier attributes the features take, in order -- for example
("radius", "gamma").
metadata : dict, optional
Provenance: how the model was trained and how well it validated.
Forward pass: tanh hidden layers, linear output.
Runs in single precision, which cuts its cost by about a third. The rounding this adds, ~1e-6 on an output of order unity, is three orders of magnitude below the model's own error; the stored weights stay in double precision.
The model shipped for a barrier class, or None if there is none.
Read through importlib.resources, so it works from a source
checkout, an installed wheel and a frozen application alike. Loaded once
and shared.
Make a feature function available to models that record name.
function(barrier_height, total_height, field, *parameters) must return
an (n, k) array, where parameters are the model's extra barrier
attributes in order.
Feature vector for a planar Schottky-Nordheim barrier.
The first three entries are the scaled physical parameters. The rest are
physics-derived: f = k_e F / h**2 is the scaled barrier field, equal to 1
exactly at the barrier top, so tanh(log f) and tanh(f - 1) give the
network a coordinate that locates the top regardless of where the other
parameters put it. The products let it represent the leading cross terms
without having to build them from scratch.
Parameters
barrier_height : array_like
h = E_F + phi - E, in eV.
total_height : array_like
W = E_F + phi, in eV.
field : array_like
Field, V/nm.
Returns
np.ndarray, shape (n, 8)
Feature vector for the curved-tip barrier: the planar features plus seven terms in tip radius and field enhancement.
The radius enters as its inverse, 20/R: the curvature corrections are a
series in x/R, so the transmission is smooth -- to leading order linear
-- in 1/R, and the planar limit R -> infinity is the finite point 0.
It is scaled to [-1, 1] over R = 20-1000 nm, the radii over which
~getelec.potential_barrier.SmallRadiiPotential is valid; gamma is
taken in its logarithm, scaled to [-1, 1] over 1-200.
Parameters
barrier_height, total_height, field : array_like
As in get_schottky_features().
radius : array_like
Tip radius, nm.
gamma : array_like
Field enhancement factor.
Returns
np.ndarray, shape (n, 15)