Reference¶

Author:	Pauli Virtanen <pauli@ltl.tkk.fi>
Organization:	Low Temperature Laboratory, Helsinki University of Technology
Date:	2005-2006

A solver for Keldysh-Usadel 1D circuit equations.

To find out how this library can be used, you should peek at

solver.Geometry: How to specify a geometry
currents: Simple interface to calculating currents
selfconsistentiteration: A self-consistent iteration
solver: Low-level solver interface

Geometry¶

class usadel1.Geometry(nwire, nnode, x=None)¶

Geometry of the setup.

Parameters:	nwire : int Number of wires in the setup nnode : int Number of nodes in the setup x : array of floats, optional Array of x-positions at which quantities are specified. These must be in the range [0, 1].

Examples

A three-probe structure:

           3
    0 S----*----S 1
        0  |  1
           |
           |2
           |
           N
           2

>>> g = Geometry(nwire=3, nnode=4)
>>> g.t_type = [NODE_CLEAN_S_TERMINAL,
...             NODE_CLEAN_S_TERMINAL,
...             NODE_CLEAN_NODE]
>>> g.w_type = WIRE_TYPE_N
>>> g.w_ends[0,:] = [0, 3]
>>> g.w_ends[1,:] = [1, 3]
>>> g.w_ends[2,:] = [2, 3]
>>> g.w_length = [0.5, 0.5, 2]
>>> g.t_delta = [100, 100, 0]
>>> g.t_phase = array([-.5, .5, 0]) * pi/2
>>> g.t_t = 1
>>> g.t_mu = 0

coupling_lambda¶: array(nwire, nx) of $\lambda$ in wires

equal(other, compare_delta=True, compare_kinetic=True, compare_phases=True)¶

Compare two geometries.

Parameters

other : Geometry: The geometry to compare this one to.
compare_kinetic : bool, optional: Whether to compare kinetic quantities (temperatures, potentials)
compare_delta : bool, optional: Whether to compare $|\Delta|$ between the Geometries.
compare_phases : bool, optional: Whether to compare ${\rm arg}\Delta$ between the Geometries.

get_idstr()¶: Get a descriptive string for this geometry.

get_node_wires()¶: Collect lists of wires connected to each node, and label nodes by which end of the wires they correspond to.

nnode¶: how many nodes in geometry

nwire¶: how many wires in geometry

omega_D¶: array(nwire, nx) of Debye frequency in wires

t_delta¶: array(nnode) of energy gaps of nodes

t_inelastic¶: array(nnode) of $\Gamma$ -s of nodes

t_mu¶: array(nnode) of node potentials

t_phase¶: array(nnode) of phases of nodes

t_spinflip¶: FIXME: this is not yet functional array(nnode) of $\Gamma_{sf}$ -s of node

t_t¶: array(nnode) of node temperatures

t_type¶

array(nnode) of types of nodes

Possible choices are: NODE_CLEAN_N_TERMINAL (normal terminal; clean interface), NODE_CLEAN_S_TERMINAL (superconducting terminal; clean interface), NODE_CLEAN_NODE (node connecting several wires; clean interface), NODE_FREE_INTERFACE (vacuum interface).

w_conductance¶: array(nwire) of conductance*area of wires

w_delta¶: array(nwire, nx) of energy gaps in wires

w_ends¶: array(nwire,2) that maps wire ends -> nodes

w_inelastic¶: array(nwire) of $\Gamma$ -s of wires

w_length¶: array(nwire) of wire lengths

w_phase¶: array(nwire, nx) of phase in wires

w_phase_jump¶: array(nwire) phase jump between the ends of the wire, due to vector potential parallel to the wire.

w_spinflip¶: array(nwire) of $\Gamma_{sf}$ -s of wires

w_type¶

array(nwire) of wire types

Possible choices are: WIRE_TYPE_N (normal wire) and WIRE_TYPE_S (superconducting wire).

x¶

array(nx) of discretization points.

NOTE: These must be in the range [0,1] !

Note also that this does not affect the actual mesh chosen by all solvers – however, it always affects the discretization used for spectral solutions and kinetic coefficients.

Solver¶

class usadel1.Solver¶

Low-level interface to the Usadel solver.

set_geometry(g, preserve=False)¶

Set the geometry to use.

Parameters:	g : Geometry The geometry to use preserve : bool, optional Whether to preserve currently set values for $\Delta$ and the kinetic coefficients.

set_kinetic(coefs)¶

Set the coefficients for the kinetic equations.

Parameters:	coefs : KineticCoefficients Kinetic coefficients for the equations.

sp_solve(E, x, continued=False)¶

Solve spectral equations.

Parameters:	E : array of floats Energies to solve the equations at x : array of floats, optional Positions to return quantities at. Defaults to the same as those in Geometry. continued : bool, optional Whether to use old solution as a starting point.
Returns:	sol : SpResult Solution to the equations.

kin_solve(E, x, continued=False)¶

Solve kinetic equations.

Parameters:	E : array of floats Energies to solve the equations at x : array of floats, optional Positions to return quantities at. Defaults to the same as those in Geometry. continued : bool, optional Whether to use old solution as a starting point.
Returns:	sol : KinResult Solution to the equations.

Results¶

class usadel1.SpResult(E=None, x=None, r=None)¶

The result from a spectral calculation.

The parameters a and b correspond to the representation

$\hat{G}^R = \frac{1}{1 - a b} \begin{pmatrix} 1 + a b & 2 a \\ -2 b & -(1 + a b) \end{pmatrix}$

of the retarded Green function.

Attributes:

E : array of floats, shape (ne,): Energies the solutions are evaluated at
x : array of floats, shape (nx,): Positions the solutions are evaluated at; scaled to [0, 1]
a : array of floats, shape (ne, nwire, nx): Riccati parameter a
b : array of floats, shape (ne, nwire, nx): Riccati parameter b
da : array of floats, shape (ne, nwire, nx): Derivative of a
db : array of floats, shape (ne, nwire, nx): Derivative of b
ne : int: Number of energy points
nx : int: Number of x-points
nwire : int: Number of wires

class usadel1.KinResult(E=None, x=None, r=None)¶

The result from a kinetic calculation.

Attributes:

E : array of floats, shape (ne,): Energies the solutions are evaluated at
x : array of floats, shape (nx,): Positions the solutions are evaluated at; scaled to [0, 1]
fL : array of floats, shape (ne, nwire, nx): The distribution function $f_L$
fT : array of floats, shape (ne, nwire, nx): The distribution function $f_T$
dfL : array of floats, shape (ne, nwire, nx): Derivative of fL
dfT : array of floats, shape (ne, nwire, nx): Derivative of fT
jL : array of floats, shape (ne, nwire, nx): The spectral current $j_L$
jT : array of floats, shape (ne, nwire, nx): The spectral current $j_T$
ne : int: Number of energy points
nx : int: Number of x-points
nwire : int: Number of wires

class usadel1.KineticCoefficients(g=None, r=None, to_evaluate=None)¶

Coefficients in the kinetic equations.

Attributes:

x : array of floats, shape (nx,): x-positions the coefficients are evaluated at; scaled to range [0, 1]
E : array of floats, shape (ne,): Energies the coefficients are evaluated at
DL : array, shape (ne, nwire, nx): Coefficient $D_L$
DT : array, shape (ne, nwire, nx): Coefficient $D_T$
TT : array, shape (ne, nwire, nx): Coefficient ${\cal T}$
rjE : array, shape (ne, nwire, nx): Coefficient $\Re j_E$
ijE : array, shape (ne, nwire, nx): Coefficient $\Im j_E$
dDL : array, shape (ne, nwire, nx): Derivative of DL
dDT : array, shape (ne, nwire, nx): Derivative of DT
dTT : array, shape (ne, nwire, nx): Derivative of TT
drjE : array, shape (ne, nwire, nx): Derivative of rjE
dijE : array, shape (ne, nwire, nx): Derivative of ijE
cTL : array, shape (ne, nwire, nx): prefactor of $f_L$ in the sink term for $j_T$
cTT : array, shape (ne, nwire, nx): prefactor of $f_T$ in the sink term for $j_T$

Currents¶

class usadel1.CurrentSolver(geometry, E=None, chunksize=None, output_function=None, automatic_energy=False, maxE=None, ne=None)¶

Solver for currents flowing in a given geometry.

Parameters:

geometry : Geometry: The geometry to solve for
E : array of floats, optional: Energy points to evaluate all quantities at. If omitted, the default is to use a sensibly chosen grid.
maxE : float, optional: If E is omitted: the maximum energy for the automatic energy grid.
ne : int, optional: If E is omitted: the number of energy points to use.
automatic_energy : bool, optional: Whether to choose an energy grid automatically. (Yes, if E is omitted.)
chunksize : int: How often to display the progress of calculation.
output_function : func(message): Function to use to print any output. If omitted, print to stderr

Attributes:

geometry : Geometry: Current geometry
spectral : SpResult: Solution to the spectral equation
kinetic : KinResult: Solution to the kinetic equations
coefficient : KineticCoefficient: Coefficients in the kinetic equations
solver : Solver: The low-level solver.
G : array of floats, shape (ne, nw, 2, 2, nnode): The spectral conductance/thermoelectric matrix. See eg. [1].

References

[1]	Applied Physics A, 89, 625-637 (2007)

Solving equations

solve()¶: Solve both spectral and kinetic equations.

solve_spectral()¶

Solve the spectral equations.

Store the results to self.spectral and self.kinetic.

solve_kinetic(E=None, ne=None)¶

Solve the kinetic equations and store the result to self.kinetic.

Parameters:	E : array of floats, optional The energy discretization to use. None indicates that use either the same as for spectral, or, if self.automatic_energy is `True`, pick a sensible choice. ne : int, optional How many points to use in energy discretization, if choosing the energy discretization automatically.

solve_spectral_if_needed(calculate_G=True)¶

Solve for spectral data if there is none yet.

Parameters:	calculate_G : bool, optional Whether to calculate the spectral conductance/thermoelectric matrix.

solve_spectral_and_save_if_needed(filename, calculate_G=True, **kw)¶

Solve for spectral data if there is none yet, and save the result to a file.

Parameters:

filename : str: Name of the file to solve data to.
calculate_G : bool, optional: Whether to calculate the spectral conductance/thermoelectric matrix.

Other parameters as in save.

calculate_G(superconductors_in_equilibrium=False, only_terminals=None)¶

Compute the spectral conductance matrix.

Warning

Spectral quantities (such as G) are not conserved in superconducting wires. Be aware that the G is evaluated at x[0] in each wire.

Parameters:	superconductors_in_equilibrium : bool, optional Assume superconductors are at equilibrium, and skip calculating the conductance matrix entries for them. only_terminals : list of int, optional Terminals for which to calculate the conductance matrix entries. If omitted, calculate for all terminals.

approximate_G(use_jS=True, use_T=True)¶

Find an approximation for the G-matrix, up to first order in $j_S$ and ${\cal T}$ .

Stores the result to self.G.

Parameters:	use_jS : bool, optional Use the spectral supercurrent in the approximation. use_T : bool, optional Use the coefficient ${\cal T}$ in the approximation.

set_solvers(*args, **kwargs)¶

Set solver types and tolerances.

Parameters:	sp_solver : int, optional Spectral solver to use. Choices are `SP_SOLVER_COLNEW` (default) and `SP_SOLVER_TWPBVP`. sp_tol : float, optional Tolerance to use in the spectral solver. kin_solver : int, optional Kineticsolver to use. Choices are `KIN_SOLVER_COLNEW` and `KIN_SOLVER_TWPBVP` (default), and `KIN_SOLVER_BLOCK`. kin_tol : float, optional Tolerance to use in the kinetic solver.

Computing currents

get_currents(E=None, ix=None, w_jL=None, w_jT=None, ne=None)¶

Calculate currents, using fixed-grid discretization.

Parameters:

ix : int, optional: Index of position to evaluate currents at in each wire. If omitted, current is evaluated at all positions.
w_jT : list of int, optional: For which wires to compute charge current for. If omitted, compute for all wires.
w_jL : list of int, optional: For which wires to compute energy current for. If omitted, compute for all wires.
ne : int, optional: How many points to put in energy discretization, if chosen automatically.

Returns:

Ic : array of floats, shape (nwire, nx’): Charge current in each wire. Contains nan in entries that were not calculated. If ix was given, nx’ == 1, otherwise nx’ == nx.
IE : array of floats, shape (nwire, nx’): Energy current in each wire.

get_currents_lazy(w_jL=None, w_jT=None, epsabs=None, epsrel=None, ix=None)¶

Calculate currents, using an adaptive integrator.

Will solve kinetic equations at energy points where the solutions are needed.

Parameters:

epsabs : float, optional: Absolute tolerance for the currents.
epsrel : float, optional: Relative tolerance for the currents.

Also takes the same parameters as get_currents.

Returns:

Returns similar output as get_currents.

get_currents_from_G(*args, **kw)¶

Calculate currents, from the spectral conductance matrix.

Parameters:	w_jT : list of int, optional For which wires to compute charge current for. If omitted, compute for all wires. w_jL : list of int, optional For which wires to compute energy current for. If omitted, compute for all wires.
Returns:	Ic : array of floats, shape (nwire,) Charge current in each wire. Contains nan in entries that were not calculated. IE : array of floats, shape (nwire,) Energy current in each wire.

get_currents_from_G_with_f(f_func, w_jT=None, w_jL=None, ix=None)¶

Calculate currents from the G-matrix, for given distribution functions at the terminals.

Parameters:

f_func : func(E, mu, T): Function that takes an energy grid of shape (ne, 1) and potentials and temperatures [shape (1, nnode)] and returns (fL, fT) where each entry is an array of shape (ne, nnode) describing the distribution function at each energy in each terminal.
w_jT : list of int, optional: For which wires to compute charge current for. If omitted, compute for all wires.
w_jL : list of int, optional: For which wires to compute energy current for. If omitted, compute for all wires.

Returns:

Ic : array of floats, shape (nwire,): Charge current in each wire. Contains nan in entries that were not calculated.
IE : array of floats, shape (nwire,): Energy current in each wire.

Computing linear response

get_linear_response_from_G(*args, **kw)¶

Compute the linear response in currents to changes in potentials and temperatures in the terminals.

This function assumes that the distribution function in each terminal is a Fermi function.

Parameters:

w_jT : list of int, optional: For which wires to evaluate the charge current related entries.
w_jL : list of int, optional: For which wires to evaluate the energy current related entries.

Returns:

G_VV : array, shape (nwire, nnode): Conductance matrix; contains entries $dI_{c,i}/dV_j$ . Has nan at entries that were not calculated.
G_VT : array, shape (nwire, nnode): Thermoelectric matrix; contains entries $dI_{c,i}/dT_j$ .
G_TV : array, shape (nwire, nnode): Thermoelectric matrix; contains entries $dI_{E,i}/dV_j$ .
G_TT : array, shape (nwire, nnode): Thermal conductance matrix; contains entries $dI_{E,i}/dT_j$ .

get_linear_response_from_G_with_f(df_func, w_jT=None, w_jL=None)¶

Compute the linear response in currents to changes in distribution functions in the terminals.

Parameters:

df_func : func(E, mu, T): Function that takes an energy grid of shape (ne, 1) and potentials and temperatures [shape (1, nnode)] and returns (dfL_dT, dfL_dV, dfT_dT, dfT_dV) where each entry is an array of shape (ne, nnode) describing the change in the distribution function with regard to perturbation in each parameter.
w_jT : list of int, optional: For which wires to evaluate the charge current related entries.
w_jL : list of int, optional: For which wires to evaluate the energy current related entries.

Returns:

G_VV : array, shape (nwire, nnode): Conductance matrix; contains entries $dI_{c,i}/dV_j$ . Has nan at entries that were not calculated.
G_VT : array, shape (nwire, nnode): Thermoelectric matrix; contains entries $dI_{c,i}/dT_j$ .
G_TV : array, shape (nwire, nnode): Thermoelectric matrix; contains entries $dI_{E,i}/dV_j$ .
G_TT : array, shape (nwire, nnode): Thermal conductance matrix; contains entries $dI_{E,i}/dT_j$ .

Saving and loading data:

classmethod load(src, path=None, **kw)¶

Load data from an open HDF5 file.

Parameters:	src: a tables.Group where to load data from, or a file name joined with node path. Other arguments (except `geometry`) are passed on to __init__

save(to, path=None, save_coefficient=True, save_spectral=True, save_kinetic=True)¶

Save date to an open HDF5 file.

Parameters:	to: str or tables.File A file name or a parent HDF5 node to save to path : str, optional HDF5 name to save under save_coefficient : bool, optional Whether to save kinetic coefficients save_spectral : bool, optional Whether to save spectral data save_kinetic : bool, optional Whether to save kinetic data.

Notes

Geometry is saved in all cases.

classmethod resume(src, geometry, path=None, compare_E=False, **kw)¶: Load data from a data file if parameters match or return a new solver if not.

Optimizing parameters¶

usadel1.optimize_parameters_for(x0, zero_func, adjust_func, **kwargs)¶

Find a configuration where the specified quantity vanishes, by adjusting the specified variables.

Parameters:

x0 : array: Initial guess for all n parameters
zero_func : func(): Function whose zero to look for; should return a vector of the same size as x0.
adjust_func : function(x): Function to call to set current parameters

Other keyword arguments accepted by scipy.optimize.fsolve can be given.

Notes

If AlreadyConvergedException is raised by zero_func, the optimization terminates successfully with the current value.

exception usadel1.AlreadyConvergedException¶: An exception that the user can raise during optimize_parameters_for to indicate that the optimum has been already reached.

Self-consistent iteration¶

usadel1.self_consistent_realtime_iteration(currents, max_iterations=100, output_func=<built-in method write of file object at 0x2aee570cb210>, iterator=<class 'usadel1.nonlinearsolver.FixedPointSolver'>)¶

Self-consistent real-time iteration for the order parameter.

Adjusts the parameters in the given Geometry object until convergence is met.

Parameters:

solver : CurrentSolver: A solver containing the geometry in which to solve for $\Delta$
max_iterations : int, optional: Maximum number of iterations to make
output_func : func(message), optional: Function to use for printing messages. If omitted, prints to stderr.
iterator: Fixed-point iteration accellerator. Suitable choices are FixedPointSolver (Broyden solver, default) and DummyFixedPointSolver (simple fixed-point iteration).

Yields:

k : int: Number of current iteration
d : fixed-point solver object: Has the methods residual_norm() that returns a max-norm estimating the error in $\Delta$ , and relative_residual_norm() that returns the norm scaled by the initial norm.
v : float: Amount by which current conservation is violated.

Notes

In real-time iteration, some care must be taken to choose the energy grid so that it is sufficiently dense in the energy range where the peak in the $F^K$ -function will be.

The energy cutoff is eliminated by using the formula:

$\Delta = [\int_0^{\epsilon_0}d\epsilon F^K_0(\epsilon)]^{-1} \left( \Delta_0 \int_0^{\epsilon_0}d\epsilon F^K(\epsilon) + \int_{\epsilon_0}^{\omega_c}d\epsilon [\Delta_0 F^K(\epsilon) - \Delta F^K_0(\epsilon)] \right)$

and in the latter term taking the bulk expressions for $F^K$ and $F^K_0$ , the latter being the bulk function corresponding to $\Delta_0$ (the bulk gap). Here, $\epsilon_0$ is the energy limit up to which numerical solutions are calculated for $F^K$ .

Examples

>>> solver = CurrentSolver(geometry)
>>> it = self_consistent_realtime_iteration(solver)
>>> for k, d, v in it:
>>>     if d.residual_norm() < 0.1 and v < 1e-4:
>>>         break
>>> else:
>>>     raise RuntimeError("Iteration did not converge")

usadel1.self_consistent_matsubara_iteration(geometry, max_iterations=100, max_ne=300, output_func=<built-in method write of file object at 0x2aee570cb210>, iterator=<class 'usadel1.nonlinearsolver.FixedPointSolver'>, E_max=None, force_integral=False)¶

Self-consistent Matsubara iteration for the order parameter.

Adjusts the parameters in the given Geometry object until convergence is met. Applicable only to equilibrium situations.

Parameters:

geometry : Geometry

The geometry in which to solve for $\Delta$

max_iterations : int, optional

Maximum number of iterations to make

max_ne : int, optional

Number of Matsubara frequencies to sum over.

Note

If the number of frequencies given is smaller than that required to sum to frequencies over $\Delta$ , an energy integral is performed instead of a discrete sum.

output_func : func(message), optional

Function to use for printing messages. If omitted, prints to stderr.

iterator

Fixed-point iteration accellerator. Suitable choices are FixedPointSolver (Broyden solver, default) and DummyFixedPointSolver (simple fixed-point iteration).

E_max : float

Cutoff energy to use (smaller than Debye, but larger than energy scales of the structure).

If None, determined automatically. The automatic cutoff satisfies: - It is no larger than omega_D. - The corresponding length scale is smaller than 0.05 of shortest wire

or 0.1 of coherence length

force_integral : bool

Force integration, even if max_ne is large enough for summing up to E_max.

Yields:

k : int: Number of current iteration
d : fixed-point solver object: Has the methods residual_norm() that returns a max-norm estimating the error in $\Delta$ , and relative_residual_norm() that returns the norm scaled by the initial norm.
v : float: Amount by which current conservation is violated. (NB: not currently implemented in Matsubara iteration, always zero.)

Notes

The energy cutoff is eliminated by using the formula:

$\Delta = [\sum_{\omega_k < \epsilon_0} F_0(\omega_k)]^{-1} \left( \Delta_0 \sum_{\omega_k < \epsilon_0} F(\omega_k) + \sum_{\epsilon_0 < \omega_k < \omega_c} [\Delta_0 F(\omega_k) - \Delta F_0(\omega_k)] \right)$

and in the latter term taking the bulk expressions for $F$ and $F_0$ , the latter being the bulk function corresponding to $\Delta_0$ (the bulk gap). Here, $\epsilon_0$ is the energy limit up to which numerical solutions are calculated for F.

Examples

>>> it = self_consistent_matsubara_iteration(geometry)
>>> for k, d, v in it:
>>>     if d.residual_norm() < 0.1 and v < 1e-4:
>>>         break
>>> else:
>>>     raise RuntimeError("Iteration did not converge")

HDF5 file layout¶

The HDF5 files produced by CurrentSolver.save() have the layout:

/geometry/w_type
/geometry/t_type
...
/spectral/a
/spectral/b
...
/kinetic/fL
/kinetic/fT
...
/coefficient/DL
/coefficient/DT
...

i.e. the attributes of the Geometry, SpResult, KinResult, and KineticCoefficient objects are simply dumped to the file to an appropriate place.

You can easily load the HDF5 files by using the supplied h5load.m script (probably requires a recent version of Matlab).

Table Of Contents

Previous topic

Next topic

Reference¶

Geometry¶

Solver¶

Results¶

Currents¶

Optimizing parameters¶

Self-consistent iteration¶

HDF5 file layout¶

Navigation

Table Of Contents

Previous topic

Next topic

Quick search

Reference¶

Geometry¶

Solver¶

Results¶

Currents¶

Optimizing parameters¶

Self-consistent iteration¶

HDF5 file layout¶

Navigation