wavepacket.solver

This module contains classes to solve the equations of motion set up by the classes in the wavepacket.expression module.

Classes

`ChebychevSolver`	Solver that expands the time evolution operator into Chebychev polynomials.
`RelaxationSolver`	Solver for propagation in imaginary time.
`OdeSolver`	A solver that uses scipy's ODE integrators as backend.
`SolverBase`	Abstract base class for all differential equation solvers.

Functions

diagonalize(→ Iterator[tuple[float, ...)

Calculates the eigenstates and -values of an operator.

Package Contents

class wavepacket.solver.ChebychevSolver(expression: wavepacket.expression.ExpressionBase, dt: float, spectrum: tuple[float, float])

Bases: wavepacket.solver.solverbase.SolverBase

Solver that expands the time evolution operator into Chebychev polynomials.

This solver is fast and accurate, but works only for time-independent closed systems, and requires explicit bounds for the spectrum of the Hamiltonian or Liouvillian, which is more setup work.

There is a dedicated tutorial for this solver, see Using Chebychev solvers.

Parameters:

expression: wp.expression.Expression: The Schrödinger equation or Liouvillian that describes the right-hand side of the differential equation.
dt: float: The time step that is propagated in one go.
spectrum: tuple[float, float]: Lower and upper bound of the spectrum of the Hamiltonian or Liouvillian. If the bound is generous, the solver is less efficient. If the spectrum extends beyond the bounds, the solution may diverge.

Attributes:

alpha: float, readonly: The Kosloff alpha number “spectral range * dt / 2”. Ideally, it should be around 40. Smaller values reduce efficiency, much larger values can cause problems with overflows.
order: int, readonly: The order of the expansion.

Notes

The name of this class deviates from the nowadays official transcription, which is “Chebyshev”. The original references by Tal-Ezer and Kosloff, however, use the outdated spelling, so we tried to be consistent with these references.

alpha: Final[float]

order: Final[int] = 0

step(state: wavepacket.grid.State, t: float) → wavepacket.grid.State

Evolves the given state for one elementary time step.

Given a wave function or density operator at time t, this function should return the propagated wave function or density operator at the new time t+dt, where dt is the elementary time step.

Parameters:

statewp.grid.State: The state to be evolved in time
tfloat: The time at which the time evolution starts

Returns:

wp.grid.State: The state at the new time t+dt.

class wavepacket.solver.RelaxationSolver(hamiltonian: wavepacket.operator.OperatorBase, dt: float, spectrum: tuple[float, float])

Bases: wavepacket.solver.solverbase.SolverBase

Solver for propagation in imaginary time.

This class is similar to the wp.solver.ChebychevSolver, but in imaginary time. It relaxes an initial wave function to the ground state of a system or a density operator to the density operator of a system at finite temperature.

There is a dedicated tutorial for this solver, see Ground-state relaxation and finite temperatures.

Parameters:

hamiltonian: wp.operator.OperatorBase: The Hamiltonian operator that describes the right-hand side of the differential equation.
dt: float: The time step that is propagated in one go.
spectrum: tuple[float, float]: Lower and upper bound of the spectrum of the Hamiltonian or Liouvillian. If the bound is generous, the solver is less efficient. If the spectrum extends beyond the bounds, the solution may diverge.

Attributes:

alpha: float, readonly: The Kosloff alpha number “spectral range * dt / 2”. Ideally, it should be around 40. Smaller values reduce efficiency, much larger values can cause problems with overflows.
order: int, readonly: The order of the expansion.

alpha: Final[float]

order: Final[int] = 0

step(state: wavepacket.grid.State, t: float) → wavepacket.grid.State

Evolves the given state for one elementary time step.

Given a wave function or density operator at time t, this function should return the propagated wave function or density operator at the new time t+dt, where dt is the elementary time step.

Parameters:

statewp.grid.State: The state to be evolved in time
tfloat: The time at which the time evolution starts

Returns:

wp.grid.State: The state at the new time t+dt.

class wavepacket.solver.OdeSolver(expr: wavepacket.expression.ExpressionBase, dt: float, **kwargs)

Bases: wavepacket.solver.solverbase.SolverBase

A solver that uses scipy’s ODE integrators as backend.

This solver merely wraps scipy.integrate.solve_ivp() in a Wavepacket solver. ODESolvers tend to be workhorses that can be applied to most systems without a second thought at the expense of performance. See [1] for more details on the convergence of solvers.

Parameters:

exprwp.expression.ExpressionBase: The right-hand side of the differential equation that encapsulates the quantum system.
dtfloat: The elementary time step for the time evolution.
**kwargs: Additional keyword parameters can be supplied that are directly forwarded to scipy. The most important parameters are “method” (defaults to “RK45”) and “rtol”, “atol” for the relative and absolute error tolerances (booth default to 1e-6).

References

[1]

https://sourceforge.net/p/wavepacket/cpp/blog/2021/04/convergence-2

step(state: wavepacket.grid.State, t: float) → wavepacket.grid.State

Evolves the given state for one elementary time step.

Given a wave function or density operator at time t, this function should return the propagated wave function or density operator at the new time t+dt, where dt is the elementary time step.

Parameters:

statewp.grid.State: The state to be evolved in time
tfloat: The time at which the time evolution starts

Returns:

wp.grid.State: The state at the new time t+dt.

class wavepacket.solver.SolverBase(dt: float)

Bases: abc.ABC

Abstract base class for all differential equation solvers.

Each solver must take in the constructor the size of an elementary time step, and it must implement a method step() to evolve a state by one elementary time step forward in time.

Usually, a solver does not need to know what step it evolves or how a differential equation looks in detail. The form of the differential equation is encapsulated in a wavepacket.expression.ExpressionBase instance that is usually supplied to the specific implementations.

A solver may, however, implicitly assume properties of the differential equation. It may also take additional arguments, these details are documented in the specific implementations.

Parameters:

dtfloat: The size of an elementary time step.

Raises:

wp.InvalidValueError: If the timestep is not positive.

dt: Final[float]

abstractmethod step(state: wavepacket.grid.State, t: float) → wavepacket.grid.State

Evolves the given state for one elementary time step.

Given a wave function or density operator at time t, this function should return the propagated wave function or density operator at the new time t+dt, where dt is the elementary time step.

Parameters:

statewp.grid.State: The state to be evolved in time
tfloat: The time at which the time evolution starts

Returns:

wp.grid.State: The state at the new time t+dt.

propagate(state0: wavepacket.grid.State, t0: float, num_steps: int, include_first: bool = True) → Iterator[tuple[float, wavepacket.grid.State]]

Generator function that yields the propagated wave function at multiple time steps.

This function allows you to propagate in one go with a for loop. It repeatedly calls step() and returns the wave function and current time.

Parameters:

state0wp.grid.State: The initial state to be propagated in time.
t0float: The initial time at which the state is given
num_stepsint: For how many elementary time steps the state should be propagated.
include_firstbool: If true, the function will start by yielding the initial state.

Yields:

Tuple[float, wp.grid.State]: A tuple consisting of the time and the result of the propagation at that time. The time starts optionally with t0 and progresses in units of dt.

Raises:

wp.InvalidValueError: If num_steps is negative.

Examples

>>> solver = ...
>>> psi0 = ...
>>> for time, psi in solver.propagate(psi0, t0, 5):
>>>    print(f't = {time}, trace = {wp.grid.trace(psi)}')

wavepacket.solver.diagonalize(op: wavepacket.operator.OperatorBase, t: float | None = None) → Iterator[tuple[float, wavepacket.typing.ComplexData]]

Calculates the eigenstates and -values of an operator.

See Calculating eigenstates and -energies of operators for a discussion of this topic.

This function is a wrapper around numpy.linalg.eigh that calculates a matrix representation of the operator transforms the calculated eigenstates into a wp.grid.State for easier consumption, and provides a generator for looping instead of a matrix with all eigenvalues in one go. This function diagonalizes a full, dense operator matrix, so it requires

Typically, you solve the eigenproblem for time-independent operators, but you can also calculate instantaneous states and energies by specifying a time value.

Parameters:

op: wp.operator.OperatorBase: The operator whose eigenstates and -values are calculated.
t: float | None = None: The time at which the operator is evaluated. Required only for time-dependent operators.

Yields:

Tuples consisting of the eigenenergy and the eigenstate of the operator.
The output is sorted by the eigenvalues.

Raises:

wp.InvalidValueError: If no time was supplied for a time-dependent operator

Examples

Iterate over the eigenvalues and -vectors

>>> hamiltonian = ...
>>> for energy, state in wp.solver.diagonalize(hamiltonian):
>>>     print(f'E = {energy}, trace norm = {wp.grid.trace(state)}')